Homework

The material in the homework is part of the curriculum.

• The homework is compulsory.
• These will count for 5% of the final grade whether you do the homework or not.
• The deadline for handing in is one week after supervision.
• Each exercise set will only be evaluated as "approved/not approved", no further classification or feedback.
• You are allowed to work together with others, but you will hand in your homework individually.
• For each exercise set it is specified what to hand in. Do not hand in a detailed report or any Matlab code.

Exercise classes will be held in Nullrommet (unless announced otherwise). If you don't have access to Matteland, you need to send your full name, NTNU username and the card number on your student identity card to Sølve as soon as possible.

17.01., 20.01. 23.01. Boundary value problems Implementation of the numerical solution of simple boundary value problems with various boundary conditions. Here are some examples of implementations of other problems in Matlab that you may use as inspiration and reference. Numerical solution of a convection diffusion boundary value problem in 1D ($\nu u''-u'=f, u(0)=\alpha, u(1)=\beta$) convdiff.m. Numerical solution of the pendulum problem ch. 3.3 of the note pendulum.m. You find here an example of how to compute error plots giving numerical evidence that the method is implemented correctly and achieves the expected order.
20.01., 25.01., 01.02. 03.02. Parabolic equations Implement a finite differences discretization of the heat equation, with various boundary conditions: use the Euler method, the backward-Euler method and the Crank-Nicholson method. As a reference solution, you may use the method you are studying with a very small time step and the same spatial step when measuring the error in time, and a very small spatial step and the same time step when measuring the error in space. If your numerical method is NOT unconditionally stable then the time step must be chosen according to the stability restriction and the size of h, so for a tiny h you'll have to use an even smaller k. This becomes a challenging numerical experiment, but still doable, unless the exact solution of the problem is available. You should consider the piecewise constant error function as in the previous exercise, for a fixed time. That is, run your method until time T and find the error from the difference between you numerical solution and your reference solution in the spatial nodes.
25.01., 01.02., 08.02. 10.02. Elliptic equations Implement a discretization of the Laplace equation on the square.
08.02., 10.02., 15.02. 17.02. Advection equation and hyperbolic problems Simple experiments with the linear advection equation and the Burgers equation. Semi-discretization of the KdV equation: Kdv.m Kdvfunc.m Kdvinit.m.
From 15th of February until the 29th of March supervision of the project

Feedback on the homework

1) Boundary value problems

Overall, most of you have solved these tasks satisfactory. The most common errors are listed below.

• For the convergence plots, use logarithmic scales on both axes. In Matlab: If using “hold on” to plot more than one graph in the same figure, call this after calling loglog for the first time. Otherwise, you will not get logarithmic axes.
• If your numerical solution is not symmetric around x = 0.5, you should check your matrix A and vector F. If you are new to Matlab, note that indices start at 1.
• Remember to scale the norms correctly when calculating the function norm from the vector norm for the error (see p. 16 in the lecture notes). In task 2 all norms should give an error plot with slope 2.
• To visualize the correct order in a convergence plot, it is illustrative to include a straight line with the expected slope to see if this is parallel to the slope of the error.
• Note that in task 3, if you use f = sin(π x), you will get f_0 = 0, and thus case 1 and case 2 will give identical schemes. Hence the order of convergence for case 1 will be 2 in this specific case.
• Please, only hand in one file (pdf, doc, odt, png, jpeg. Not .fig and not .zip).

2) Parabolic equations

Most of you seem to have had a few or many difficulties with this exercise. Below are some common errors and hints on how you should solve each of the tasks.

• You should find the error in a function norm as in the previous exercise. That is, use a modified version of the vector norm defined on the spatial nodes for a given time. Do not integrate over time (and hence do not multiply your error vector with the square root of the time step). (See section 5.2 in the lecture notes.)
• If you compare the numerical solution with the exact solution, you must use a small spatial step on your method when studying the convergence in time, and a small time step when studying the convergence in space. You will however have problems with verifying the convergence in time for Forward Euler. This is due to the stability criterium r ≤ 1/2 for this method (see section 5.6 in the lecture notes). The solution to this problem is to compare with a solution obtained using the same method, with the same spatial step, but with a very fine time step.
• Remember that your matrix-vector system is a system of equations. In Task 3, you have one or two (depending on whether you choose to consider Neumann conditions on one or both boundaries) unknowns, and you must thus modify your system to include one or two more equations. The equations determining the values at the boundary should be as when you implemented Neumann conditions in last weeks exercise, while the rest of the equations should be as in task 1 here.
• In the last task, when implementing periodic boundary conditions, you should use your code from Task 1 and make small adjustments to your matrix A while getting an extra element in your u-vector. To better understand this problem, you could think of the domain as a cylinder. Note that you do not know the value of u on the boundaries x = 0 and x = 1 in this case, but you know u(0,t) = u(1,t). When you discretise to find an approximation to the second derivative in x=0, you can use the neighbouring values larger than 0 and smaller than 1.
• As in the previous exercise, you should include a reference line with your convergence plots, to verify the order of your methods.

3) Elliptic equations

• The intended way to solve the problem was to construct a linear system of equations as described in the note on page 66. From the non-homogeneous boundary at y=1 you will get some non-zero elements in the vector on the right hand side.
• When you want to analyze the convergence in the spatial directions seperately by comparing to the exact solution, you need to ensure that the fixed stepsize is small compared to the varying step size throughout the analysis. Otherwise, the main source of error will come from the spatial direction of the fixed step size. Alternatively, you could use as a reference solution a solution obtained by your method using a small step size in the direction you are considering.
• Many of you have only showed the rate of convergence in one direction, or for h = k. You should verify the correct order of convergence in both directions. To do this, you would need to implement the method to also work for h ≠ k.
• As before, remember to include a reference line close to your error curves in the convergence plots, and use the loglog() function in Matlab for plotting these plots.
• Some of you had handed in 2D line or contour plots of the numerical solution, with no labels on the u values. You should either include labels on u or plot the solution in three dimensions, e.g. by using the surf function in Matlab.

4) Advection equation and hyperbolic problems

• Exercise 1: The oscillation most of you observed to the right, is not due to the stability condition of the explicit scheme (k < h^2/(2*nu)). In fact, if you implement the corresponding implicit scheme (which is unconditionally stable) you can still observe this effect. The phenomena is due to under resolution in the spatial direction where the gradient/derivative of the solution is large.
• Exercise 2: Remember that the domain in the spatial direction is unbounded. So when you want to solve the problem on some finite interval, you should consider the corresponding domain of dependence as described in the note. Either employ periodic boundary conditions or consider an initial solution and domain in space and time such that you can use uniform Dirichlet boundary conditions.
• If you do not clearly observe the unstability of (7.12), your maximum time is too low. No matter the discretisation parameters, (7.12) will diverge with increasing time.
• You should observe second order convergence both in space and time for both the Lax-Wendroff and Leap-Frog schemes (see p. 85 of the note).

Getting started with MATLAB

We expect that you have basic programming skills in Matlab and that you have implemented some numerical algorithms before. If not, you should read up and practice on your own. You will find tons of tutorials and introductory notes on the web. Below are some resources to get you started.