Numerical Differential Equations

• Read the code in nrrbasic.m
• Create any matrices 3x4 A i 3x5 B - and then a matrix 3x8 C whose first 3 columns are A and next B. Extract from C a block D from C(1,1) to C(3,3). Flip colums of D. Write D to a file. (binary or ASCII) - change D(1,1) to -100 in the ascii-file and read the new matrix to octave'a. Compute norms of D e.g. 1st, 2nd, Frobenius ones etc.
• Compute discrete max norm of (sin(x))^2 on [0,1] (without using loops).
• Write Euler's schemes for y'=ay y(0)=1, a=1,10,100,-1,-10,-100 on [0,T] - write graphs for T=1,10,100, compute errors: ae(T;h)=|y(T)-y_h(T)| and ce(T;h)=ae(T;h)/y(T) and ae(h;h)=|y(h)-y_h(h)| i ce(h,h)=ae(h;h)/y(h). Print ae(t;h) i ce(t;h) for t\in[0,T] on the display.
• Repeat the problem for y''=-y z y(0)=0;y'(0)=1 - draw the computed orbits (y,y'). Are they contained in a circle or periodical?
• Do the same but for y'=1+y; y(0)=0 on [0,T) for T=0.5,1,1.5. Compute error t=0.5,1.5 for h=0.5,1/4,1/8 etc.

• read help for lsode - draw solutions for y'=-0.1y y(0)=1 on [0,20]using lsode(). Do the same for y'=(cos(x))^2 z y(0)=0 on [0,T] for T=1,10,20 or y'=1+y*y, y(0)=0.
• Draw graphs of solution and its derivative and of orbit (y(t),y'(t)) for y''=-y y(0)=0 y'(0)=1. Check if (y(T)^2+(y'(T)^2==1 for large T and different h.
• Implement midpoint and Taylor schemes of order 2 for y'=ay, y(0)=1, a=1,10,100,-1,-10,-100 on [0,T] draw the graphs for T=1,10,100, compute errors for different T etc Compare the graphs and errors with the results for explicit Euler scheme
• Do the same for y''=-y, y(0)=0;y'(0)=1 - draw the orbit (y,y'). Are orbits periodical? Repeat the problem for y''=-sin(y), y(0)=0,y'(0)=1. Are computed solutions periodical? Compare with lsode() solutions.
• Test error of midpoint for y'=-y, y(0)=1 on [0,T] for growing T with h=0.1,1e-2,1e-3 etc Draw graphs.
• Test experimentally the order of local truncation error of Taylor or midpoint scheme for all IVP with known solutions: : y'=ay,y(0)=1; a=-100,-10,-1,1,10,100 for t=1,10,100 (if there is no fl overflow) do same for y''=y y(0)=0 y'(0)=1; y'=cos^2(y) y(0)=0 for t=1 i 100 itd Do same for the pendulum eq: y''=-sin(y) y(0)=0 y'(0)=1 taking the lsode() solution as the exact one.
  midpoint17.m - implementation of midpoint scheme as in the lab this year
  midpoint.m - implementation of midpoint scheme (old one)
  testmidpoint.m - testing midpoint scheme - order of convergence and instability for long time interval (for dx/dt=-x with x(0)=1) - using midpoint.m m-file
  taylor17.m - implementation of Taylor scheme as in the lab this year
  testlte17.m - tests of the order of LTE for Euler and midpoint schemes

• Explicit two-step Adams-Bashforth scheme x.
  x(n+1)=x(n)+(h/2)*[3*f(n)-f(n-1)] with f(n)=f(t(n),x(n))
  Draw graphs of solutions: y'=-y, y(0)=1 or y''=-y y(0)=0;y'(0)=1 itp (y_1 take as a solution or compute by Heun method)
• Test the order of local truncation error of explicit 2-step Adams-Bashford scheme on y'=ay, y(0)=1, : a=1,-1.
• Test the convergence order of explicit 2-step Adams-Bashford scheme on y'=ay, y(0)=1, x_1=exp(ah) : a=1,-1.,1 for T=0.1,1,10,100. (using the method of halving the step i.e. compute e(T,h)/e(T,h/2) for e(t,h)=||x(T,h)-sol(T)||; x(t,h) approximation of the soultion sol(T) computed by the scheme using the step h) Take different x_1: e.g. x_1=exp(ah), x_1 computed by Heun or explicit Euler: x_1=x_0+h*f(t_0,x_0). Do you see any difference?
• Repeat the previous problem for the Taylor(2nd order) Heun and modified Eulera schemes for y'=ay, y(0)=1, x_1=exp(ah) for a=1,-1 and t=0.1,1,10,100.
• Test the order of convergence of Heun,, implicit two-step Adams and modified Eulera schemes on y'=1+y^2, y(0)=0 for t=0.1,1,1.5, 1.54. What happens close to pi/2?
• Test the order of convergence of Heuna and modified Eulera schemes on for y''=-y'', y(0)=0 y'(0)=1 (for Adams x_a1 computed by explicit Euler) for t=0.1,1,10,100,1000. Are the computed orbits periodical? Repeat the problem on y'=-sin(y), y(0)=0, y'(0)=1.
• Write predictor corrector scheme - for explicit and implicit Euler schemes nonlinear system is solved by Banch iterations: x^{k+1}=x_n+h*f(x^k,t^k+h) (if x^{k+1}=x^k then x^k=x_{n+1}) for x^0=x_n+hf_n (predictor ex. Euler) - no of nonlinear iterations M should be a parameter - test for M=1,2,3 -It can be written with a stopping criteria like |x^{k+1}-x_n-h*f(x^k,t^k+h)|<= h or <= x_n*h
• Implement 2step Adams-Moulton scheme
  x(n+1)=x(n)+(h/12)*(5*f(n+1)+8*f(n)-f(n-1)) with f(n)=f(t(n),x(n))
  and test the order of convergence - use 1step Heun scheme for computing x(1).
• Implement 2step predictor - corrector scheme: take 2step Adams-Bashford as a predictor and 2-step Adams-Moulton as a corrector. Test the order of convergence
• Repeat tests for explicit 3step Adams-Bashford scheme:: x(n+1)=x(n)+(h/12)*[23*f(n)-16*f(n-1)+5*f(n-2)] with f(n)=f(t(n),x(n))
  x(1),x(2) may be computed by Heun scheme (order 2) - or substitute x(t0+h) x(t0+2h) with x(t) the solution of IVP.

• Compare computation times of X(T) for X the solution of dX/dt=AX with X(0)=[1;2] for the matrix [-1,-10;-1,-11] and T=10,100,1000,1000 etc once using lsode with stiff solvers and once with non-stiff (lsode_options("integration method","non-stiff")). To get real time use: (tic;lsode(..);toc)
• Implement shooting method for y''-d(x)y'- c(x)y=f(x), y(a)=ya; y(b)=yb, a, b, ya,yb given, c,f given functions on [a,b] (2 s shoots with y'(a) equal to s_1=0 and s_2=1i.e. we solve the equation with y(a)=ya; y'(a)=s_k - get two solutions for t=b : y(b;s_k) k=1,2 and hence we compute s such that the solution with y(a)=ya; y'(a)=s_k is such that y(b)=tb).
• Solve by the shooting method: -y''+y=0 z y(0)=1; y(b)=1 for b=1,4,10,15,20. Draw graphs. Compute errors |y(b;s)-yb|.
• Solve using the shootin method:
  1. -y''+y=f z y(0)=0; y(b)=0 for f=2*sin(x) and b=\pi (we know the solution).
  2. -y''+y=-2+(1-x*x); y(-1)=0; y(1)=0. (we know the solution).
  3. -y''+exp(-x) y =cos(x); y(-1)=0; y(1)=0.
  4. -y''+ay'+ y =cos(x); y(-1)=0; y(1)=0 for a=0,1,10,100,-1,-10,-100.
  5. -y''+(1+t*t)y=sin(x)+ (1+x*x)*sin(x) with y(0)=0 and y(1)=sin(1) (we know the solution y(x)=sin(x)).
  Plot the graphs. Compute |y(b;s)-yb|.
• (homework) Implement the shooting method for y''=F(t,y,y'), y(a)=ya; y(b)=yb. Use lsode() and fsolve() Test for F(t,y,y')=sin(y), y(0)=1, y(1)=2.
• Form a tridiagonal symmetric matrix with two on main diagonal and minus one on sub and superdiagonals: using a loop and the octave function diag() Form this matrix in a sparse format using the octave function sparse() without using diag() or forming any full format matrix.
• Solve the FDM problem -u''+u=0 , u(0)=u(T)=1 on [0,T] for the known solution u(x) (you should copmute it...) for T=1,5,10,15,20 for N=100 and 1000. Compare with the shooting method.
• Compute the order of local truncation error of standard FDM for -u''=\sin(x) , u(0)=u(\pi)=0
• Solve the FDM problem -u''=f , u(0)=u(\pi)=0 on [0,\pi] for the known solution u(x)=\sin(x) on the mesh of 10,20,40,80 points. Compute the discrete error in max and discrete L^2 norms.
• Compute the order of local truncation error of standard FDM for -u''=\sin(x) , u(0)=u(\pi)=0 in L^2 and max discrete norms.
• Compute the order of error of standard FDM for -u''=\sin(x) , u(0)=u(\pi)=0 in L^2 and max discrete norms by the halving the mesh size method. I.e we compute the error e_h and e_{h/2} and then check the ratio: e_h/e_{h/2} . Repeat for discrete H1 norm (discrete L2 norm of difference of u_h) i.e. ||u||_{1,h}^2=\sum_{k=0}^{N-1} h*|D_hu(x_k)|^2 with D_h a forward difference.

• Form the matrix and the right hand side vector for the FDM discretiation of the problem -u''=f on [0,1], u'(0)=\alpha;u(1)=\beta picking f , \alpha,\beta for a known solution e.g., u(x)=\sin(x+1) (We would like to avoid sitution that u^{(k)}(\alpha)=0 for any k which could artificially increase the order of convergence or of the local truncation error). The Neumann condition at the left end of the interval we approximate by the forward difference. Compute the FDM solution - plot the FDM solution - plot the error - compute the discrete norms of the error. Check experimentally the order of the error. Compute the error of the local truncation error.
• Compute the order of error of standard FDM for - u''=f in [-1,1] , u=0 on the boundary for the known solution u(x)=\sin(\pi*x) in L^2 and max discrete norms but on NO mesh which is ON THE BOUNDARY e.g. take the mesh (h*k) for k,l=-N-1,...,N+1 but for h=0.9/(N+1) and introduce discrete zero bnd conditions for k\in\{0,N+1\}
• Form a 5-diagonal symmetric matrix - a FDM discretization of 2D Laplacian on a regular mesh on a square.
• Solve the FDM problem -Laplacian u=f in [0,1]^2 , u=0 on the boundary for the known solution u(x)=\sin(\pi*x)\sin(\pi*y) on the mesh of 10,20,40,80 points in each direction. Compute the discrete error in max and discrete L^2 norms. Plot the error and the FDM solutions.
• (homework) Compute the order of local truncation error of standard FDM for -Laplacian u=f in [0,1]^2 , u=0 on the boundary for the known solution u(x)=\sin(\pi*x)\sin(\pi*y) in L^2 and max discrete norms.
• Compute the order of error of standard FDM for -Laplacian u=f in [0,1]^2 , u=0 on the boundary for the known solution u(x)=\sin(\pi*x)\sin(\pi*y) in L^2 and max discrete norms by the halving the mesh size method. I.e we compute the error e_h and e_{h/2} and then check the ratio: e_h/e_{h/2} .
• Compute the order of error of standard FDM for -Laplacian u +c(x,y)u=f in [-1,1]^2 , u=0 on the boundary for c(x) nonnegative and discontinuous e.g. c(x,y)=1 for x<0 ; c(x,y)=0 otherwise for the known solution u(x)=\sin(\pi*x)\sin(\pi*y) in L^2 and max discrete norms.
• Compute the order of error of standard FDM for -Laplacian u=f in [-1,1]^2 , u=0 on the boundary for the known solution u(x) which is of low regularity e.g. C^2 or C^3 - take u(x,y)=\sin(\pi*x)*f(y) with f(x) a function which is a piecewise cubic polynomial on [-1,0] and [0,1] such that f(-1)=f(1)=0 and f,f',f'' continuous at x=0 . (Hermite interpolation). Consider the mesh such that there mesh points on the line y=0 and without meshpoints on this line.
• consider a mixed bnd - i.e .on the lower edge put \partial_n u = g_1 and the remaining u=g - modify the code - test the order of convergence for a known solution (just take the matrix from Dirichlet problem and add N-1 rows/column (first if bottom edge nodes are forst or lst if bottom edge nodes are numbered at the end)
• Compute the order of error of standard FDM for -Laplacian u=f in [-1,1]^2 , u=0 on the boundary for the known solution u(x)=\sin(\pi*x)\sin(\pi*y) in L^2 and max discrete norms but on NO mesh which is ON THE BOUNDARY e.g. take the mesh (h*k,l*h) for k,l=-N-1,...,N+1 but for h=0.9/(N+1) and introduce discrete zero bnd conditions for k,l\in\{0,N+1\}
• Stability . Compute the matrix norms L^1_h, L^2_h max - A+cI and her inverse for some c (A FDM matrix of 2D Laplacian on a square - equidistant mesh) for N=10,20,40,80,160.- Attention! I can be done using the octave's function norm(). Formally we should put Dirichlet bnd into the matrix...

• Run the scripts - i.e. solve u_t=u_{xx} for discontinuous u0
• Change the script in order to solve u_t-u_{xx}=f u(t,a)=ga(t) u(t,b)=gb(t) u(0,x)=u0(x) for any functions f, ga,gb,u0.
• Using the function from the previous problem solve u_t-u_{xx}=0 with zero bnd and u0 being a peak e.g. u0(x)=1e10 on t0+[-1e-7,1e7] (t0 given point) and zero otherwise
• Change the Dirichelt bnd to Neumann (or mixed Neumann at a and Dirichlet at b)
• Implement explicit, implicit Eulers and Crank-Nicholson FDM schemes for u_t-u_{xx}=f u(t,a)=ga(a) u(t,b)=gb(b) u(0,x)=u0(x) Test them for different valuers of parameters taking f=0 and ga=gb=0. Take different u0 (sin(x), sin(10x), 'a peak' etc) Test the order of convergence in discrete max norm by the halving steps methos with respecrt to \tau and h. (for a known solution e.g. u(t)=exp(-t)sin(x) a=0, b=\pi)

Octave scripts with solution of some problems from our labs