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 Local Truncation 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 local LTE 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.

Octave scripts with solution of some problems from our labs