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.

• Write a function with an implementation of the explicit Euler scheme working for any 1st order ODEs. The INPUT parameters - a handle to the vector field, the ends of the interval, an initial value, no of discrete time points plus one. OUTPUT the vectors of discrete times and the approximations of the discrete solution.
  Test the scheme for y'=ay y(0)=1 on [0,T] with a=-1,1,10,-10,100,-100 and T=1,10,100 and h=0.1,0.01,1e-3. Plot the graphs of the approximation and exact solutions and compute errors: ae(T;h)=|y(T)-y_h(T)| and ce(T;h)=ae(T;h)/y(T). Print ae(t;h) i ce(t;h) for t\in[0,T] on the display.
• Repeat the problem for the linear penduulum model: 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.

• Plot the solutions of y'=ay y(0)=1 a=1,-1 using Taylor or midpoint schemes (e.g. N=20 T=4). Test the convergence order od Taylor and midpoint scheme for this problem (one can use the scripts from this page...)
• 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?
• Test the convergence order for Taylor(2nd order), Heun

  y(n+1)=y(n)+ 0.5h*(f(n)+f(t(n+1),y(n)+hf(n))

  and modified Eulera schemes

  y(n+1)=y(n)+ h*f(t(n)+0.5h,y(n)+0.5hf(n))

  on the IVP:
  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 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 Banach 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 or preset no of iterations e.g. one or two or three.
• 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.

• Some older problems e.g. test the order of Loca lTruncation Error (LTE) of some schems: forward Euler, cakward Euler, midpoint Heun, Taylor two step Adams-Bashforth etc For IVP: y'=ay a=1,-1 with solution y(t)=exp(a*t) and in 2D : y''=-y with y(t)=cos(t).
  How to do it e.g. compute loca lLTE for a given t=1 for forward Euler for x'=f(x,t): compute e(h,t):=||x(t+h)-x(t)-h*f(x(t),t)|| then e(h/2,t) and compute e(h,t)/e(h/2,t) - if it is 2^p - it means that for this case e(h,t) is like h^p. We can also compute e(h)=max_k e(h,t(k)) taking t(k)=t0+h*k N+1 (N*h=T-t0) equidistant points on [t0,T].
• Solve some simple problems using lsode() - plot the graphs of the solution or trajectories in 2D cases e.g.

  y'=ay y(0)=1; (penduulum equations) y''=-y (linear) or y''=-sin(y) with y(0)=1 y'(0)=0.

• 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 a predictor-corrector scheme based on Heun method (predictor) and trapezoid method (corrector) i.e. for a known x(n) compute y(n) by the Heun method

  y(n+1)=x(n)+ 0.5h*(f(n)+f(t(n+1),x(n)+hf(n))

  and then substitute y(n) into (replacing x(n+1) in the rhs of the trapezoid scheme i.e. we get

  x(n+1)=x(n)+0.5 h(f(n)+f(t(n+1),y(n+1)).

  Then test the convergence order for y'=ay a=1,-1, y(0)=1 on [0,1] and [0,10] - starting with N=20.
• One an also change the p-c scheme correcting twice i.e. computing y(n) by Heun method, then z(k+1)=x(n)+0.5 h(f(n)+f(t(n+1),z(k)) k=0,1 and z(0)=y(n+1) obtaining z(1) which is taken a new approximation of x(n+1).
  We can actually correct many times e.g. p-times i.e. taking k =0,..,p-1 or as many times as long ||z(k+1)-x(n)-0.5 h(f(n)+f(t(n+1),z(k))||>TOL taking e.g. TOL=Ch^2 for some constant C. This is equivalent to solving the equation in one step of the implicit trapezoid method by a version of the Banach iterative method. (in trapezoid method x(n+1) is a fix point of the function G(y)=x(n)+h*f(t(n+1),y) and our iteration is z(k+1)=G(z(k))).
  Implement the p-c scheme which use Heun as the predictor and correct twice (corrector - trapezoid method) and test the convergence order as before.

-u''+d(x)u' c(x)u=f with Dirichlet bnd cond : u(a)=alpha u(b)=beta and for mixed bnf conditions: u'(a)=alpha u(b)=beta

• Implement shooting method for

  y''-d(x)y'- c(x)y=f(x), y(a)=ya; y(b)=yb, a, b, ya,yb given, d,c,f given functions on [a,b]

  (2 shoots with y'(a) equal to s_1=0 and s_2=1 i.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; y(0)=1; y(b)=1 for b=1,4,10,15,20.

  Draw graphs. Compute errors |y(b;s)-yb|.
• Solve using the shooting method:
  1. Linear 2nd order BVP with sin() RHS

     -y''+y=f z y(0)=0; y(b)=0 for f=2*sin(x) and b=\pi (the solution can be computed).

  2. Linear 2nd order BVP with polynomial quadratic RHS

     -y''+y=-2+(1-x*x); y(-1)=0; y(1)=0. (the solution can be computed).

  3. Linear 2nd order BVP with nonconstant coeff.

     -y''+exp(-x) y =cos(x); y(-1)=0; y(1)=0.

  4. Linear 2nd order BVP with cos() RHS and a large coeff. and the 1st order term (large advection a)

     -y''+ay'+ y =cos(x); y(-1)=0; y(1)=0 for a=0,1,10,100,-1,-10,-100.

  5. Nonconstant coeff.

     -y''+(1+x*x)*y=sin(x)+ (1+x*x)*sin(x) with y(0)=0 and y(1)=sin(1) (the solution can be guessed i.e. y(x)=sin(x)).

  Plot the graphs. Compute |y(b;s)-yb|.
• Implement the shooting method for

  y''=F(t,y,y'), y(a)=ya; y(b)=yb.

  Use lsode() and fsolve() or fzero() (or implement your own nonlinear solver e.g. the bisection method) Test for F(t,y,y')=sin(y), y(0)=1, y(1)=2 and (2) y''=F(x,y,y')=y^2 x(0)=1 x(1)=-4 (two solutions s1 approx -6 and s2 \approx -22!)

-y''+d(x)y'+c(x)y=f(x) y(a)=al y(b)=be

y''=F(x,y,y') y(a)=al y(b)=be

y''=F(x,y,y') y(a)=al y(b)=be

• 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 compute 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 (mixed bnd conditions)

  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 [0,1] , u=sin(x) on the boundary for the known solution u(x)=\sin(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=0,...,N but for h < 1/N and introduce discrete zero bnd conditions for k=0,N. One can use the old scripts e.g. taking the old mesh but be=sol(1+0.5h) so in fact we solve -u''=f on (0,1+h/2) u(0)=0 u(1+h/2)=sin(1+h/2)
• Collatz approximation - test on - u''=f in [0,1] , u=0 on the boundary for the known solution u^*(x)=sin(x) or (1-x*x)sin(x), pick the mesh such that x_0=0 but x_N < 1 as in the previous problem THE FDM scheme:

  u_0=0, (-u_{k-1}+2u_k-u_{k+1})/(h*h)=f(x_k) k=1,..,N-1.

  The last equation - Collatz approximation:

  l1(1)=0=u^*(1) with l1 linear interpolant s.t. l1(x_k)=u_k k=N-1,N.

  We have

  l1(t)=u_{k-1}(t-x_k)/(-h)+u_k(t-x_{k-1})/h.

• Test the fdm scheme for

  - u''=f in [-1,1] , u(-1)=0 u'(1)=a

  for the known solution u^*(x)=(x+1)\sin(x) - standard equidistant mesh, standard 3point stencil inside interval, but increase the order of the scheme at the right end. We add an extra so called ghost point x_{N+1]=1+h and approximate u'(b) by the central difference getting the equation (u_{N+1}-u_{N-1})/(2h)=a. Test the convergence error in the discrete max norm (and L^2). THis scheme has sense if we can extend the solution outside the interval i.e. if the solution is regular enough up to the boundary.
• Test the fdm scheme for

  - u''+cu =f in [a,b] , u(a)=al u'(b)=be

  for the known solution e.g. u^*(x)=(x+1)\sin(x) - standard equidistant mesh, standard 3point stencil inside interval, but increase the order of the scheme at the right end. We add an extra Taylor expansion term i.e [u(b)-u(b-h)]/h=u'(b)-0.5hu''(b)+O(h^2) i.e. we approximate

  u'(b)=be by [u(b)-u(b-h)]/h+0.5hu''(b)=be

  but

  -u''(b)=f(b)-c*u(b)

  (assuming that the diff. equation is also valid at x=b) thus our last equation: by

  [u(N)-u(N-1)]/h+0.5hc*u(N)=be + 0.5*h*f(b).

  Test the convergence error in the discrete max norm (and L^2). This scheme has sense if the solution is regular enough up to the boundary.
• 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.
• 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 discrete maximum norm (matrix) of the matrix obtained in FDM for

  -Laplacian u=f in [0,1]^2 , u=0 on the boundary

  for difffrent N (or h=1/N).
  This norm gives us the stability constant for the FDM scheme in the discrete maximum norm.
• 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 + \vec{b}\nabla u=f in [-1,1]^2 ,

  u=0 on the boundary for \vec{b} given constant vector e.g. \vec{b}=[1;1] or [1;0] etc for the known solution u(x)=\sin(\pi*x)\sin(\pi*y) in L^2 and max discrete norms. The gradient of u we can approximate by the respective central differences i.e. \nabla u(x,y) \approx [u(x+h,y)-u(x-h,y),u(x,y+h)-u(x,y-h)]/(2h)

- Laplacian u=2pi^2*sin(x)sin(y)
u=0 on boundary

\int_a^b u'v' + c*u*v dx =\int_a^b f v dx \forall v \in V_h

d/dt u -au''+cu=f in (a,b)
u(t,a)=al(t), u(t,b)=be(t) boundary conditions u(0,x)=u0(x) an initial value

Octave scripts with solution of some problems from our labs

ODE schemes - simple implementations and tests

-y''+d(x)y'+c(x)y=f(x) y(a)=al y(b)=be

y''=F(x,y,y') y(a)=al y(b)=be

y''=F(x,y,y') y(a)=al y(b)=be

FDM schemes - simple implementations and tests

- Laplacian u=2pi^2*sin(x)sin(y)
u=0 on boundary

FDM schemes - simple implementations and tests

-au''+cu=f

A=(\int_a^b \phi_k'\phi_l')_{k,l=0}^N and B=(\int_a^b \phi_k\phi_l)_{k,l=0}^N

Solvers for paraboli problems

d/dt u -au''+cu=f in (a,b)
u(a)=al, u(b)=be (al be constant but can be easily changed to al(t) be(t)) u(0,x)=u0(x)