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.
  midpoint.m - implementation of midpoint scheme
  testmidpoint.m - testing midpoint scheme - order of convergence and instability for long time interval (for dx/dt=-x with x(0)=1)

• Explicit Adams-Bashforth scheme. 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 dla 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)|\leq h or \leq x_n*h
• Repeat tests for explicit 3step Adams-Bashford scheme:: x(n+1)=x(n)+(h/12)*[23f(n)-16f(n-1)+5f(n-2)] z 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 dla f=2*sin(x) i b=\pi (we know the solution).
  2. -y''+y=-2+(1-x*x) z y(-1)=0; y(1)=0. (we know the solution).
  3. -y''+exp(-x) y =cos(x) z y(-1)=0; y(1)=0.
  4. -y''+ay'+ y =cos(x) z y(-1)=0; y(1)=0 for a=0,1,10,100,-1,-10,-100.
  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...

• FEM linear 1D; regular grid; Dirichlet homogeneous bnd conditions Find FEM P_1 spproximation of -u'' +cu=f on [a,b] on a regular grid (x_k=a+k*h h=(b-a)/N) for u=\sin(\pi*x) and a=0;b=\pi using the trapezoidal rule for rhs for N=10,20,40,80,160 . Compute the L^2 , discrete \infty and H^1 norms of u_h-I_h u. Check the order by the halving method. Compare with FDM method.
• FEM linear 1D; simple irregular grid; Dirichlet homogeneous bnd conditions Find FEM P_1 spproximation of -u'' +cu=f u(a)=0=u(b) on [a,b] on [0,2\pi] for u=x^3*\sin(x) on a simple irregular grid with the 2h mesh on [0,\pi] and h mesh on [\pi,2\pi] . Compute the L^2 , discrete \infty and H^1 norms of u_h-I_h u . Check the order by the halving step method starting with h0=2\pi/40 .
• FEM linear 1D; simple irregular grid; Dirichlet homogeneous bnd conditions Find FEM P_1 spproximation of -u'' +cu=f u(a)=0=u(b) on [0,2\pi] for u=x^3*\sin(x) on the following irregular grid: take grid getting denser closer to the right end : take \{x_{n+1}\}\cup\{b\} for x_{n+1}=x_n+h_n\leq b with h_n=(b-a)/(sqrt(n)*N) and x_0=a , Compute the L^2 , discrete \infty and H^1 norms of u_h-I_h u . Check the order by doubling N and computing the ratio of the error for N and 2N .
• FEM linear 1D; regular grid; mixed homogeneous bnd conditions} Find FEM P_1 spproximation of -u''+cu=f u(0)=0; u'(b)=0 on equaidistant mesh. Compute the rhs using trapezoidal composite quadrature rule. on a regular grid for u=\sin(x) and a=0;b=1 using the trapezoidal rule for rhs for N=10,20,40,80,160 . Compute the L^2 , discrete \infty and H^1 norms of u_h-I_h u . Check the order by the halving step method. Compare with FDM methods. Note that we must take larger submatrices i.e. our FEM space is spanned by (\phi_i)_{i=1}^N (\phi_N is this extra function)
• (extra problem)FEM linear 1D; general elliptic operator regular grid; Dirichlet homogeneous bnd conditions Find FEM P_1 approximation of -(a(t)u')'+c(t)u=f(t) u(le)=u(re)=0 on equaidistant mesh. Compute the stifness matrix and the rhs vector using trapezoidal composite quadrature rule. Test the convergence order for the known solution u=\sin(t) , a=1+t^2 , c=0 on [0,\pi] . Repeat the problem taking c(t)=1+t . Weak formulation is find u \in V_h \int_le^re a(t)u'v' + c(t)uv dx=\int_le^re fvdx forall v \in V_h. i.e. we get the system Au=F with A=(\int_le^re a(t)\phi_i'\phi_j' + c(t)\phi_i \phi_j)_{i,j}
• (extra problem)FEM linear 1D; regular grid; Dirichlet homogeneous bnd conditions; higher order quadrature rule for RHS Find FEM P_1 spproximation of -u'' +cu=f on [a,b] on a regular grid for u=\sin(\pi*x) and a=0;b=\pi using the Simpson rule for rhs for N=10,20,40,80,160 . Compute the L^2 , discrete \infty and H^1 norms of u_h-I_h u . Check the order by the halving method. Compare with the results where the RHS was computed by the trapezoidal rule.

Octave scripts with solution of some problems from our labs