Numerical Differential Equations

winter semester 2014-15

Syllabus

1. ordinary differential equations (ODEs)
2. elliptic partial differential equations (PDEs)
3. evolutionary PDEs (parabolic and hyperbolic of first order)

1. one-step and linear multi-step schemes for initial ODEs problems
2. finite difference method
3. finite element method

Lecture notes

References

Text books

1. Deuflhard, Peter, Bornemann, Folkmar, Scientific Computing with Ordinary Differential Equations, Series: Texts in Applied Mathematics, Vol. 42, Springer-Verlag, New York, 2002. (theory of ODEs, ODE schemes, Boundary Value Problems in 1D) One can download a pdf file from MIMUW computers (valid Dec 2014): Springer link
2. David F. Griffiths, Desmond J. Higham, Numerical Methods for Ordinary Differential Equations, Springer-Verlag, 1st Edition, London, 2010. An elementary textbook on ODE schemes. One can download a pdf file from MIMUW computers (valid Dec 2014): Springer link
3. Claes Johnson, Numerical solution of partial differential equations by the finite element method, Cambridge University Press, Cambridge, 1987.
4. Randall J. LeVeque, Finite difference methods for ordinary and partial differential equations, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007, Steady-state and time-dependent problems. (numerical schemes for ODES, finite difference methods for elliptic and parabolic PDEs)
5. Alfio Quarteroni, Riccardo Sacco, and Fausto Saleri, Numerical mathematics, Texts in Applied Mathematics, vol. 37, Springer-Verlag, New York, 2000. (numerical schemes for ODEs and some PDES - hyperbolic nad parabolic) One can download a pdf file from MIMUW computers (valid Dec 2014): Springer Link
6. John C. Strikwerda, Finite difference schemes and partial differential equations, second ed., Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2004. (FD schemes for PDEs - all types)

Monographs or advanced text books

1. Dietrich Braess, Finite elements, third ed., Cambridge University Press, Cambridge, 2007, Theory, fast solvers, and applications in elasticity theory, Translated from the German by Larry L. Schumaker. (advanced text book)
2. Susanne C. Brenner and L. Ridgway Scott, The mathematical theory of finite element methods, third ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008.
3. J. C. Butcher, Numerical methods for ordinary differential equations, second ed., John Wiley and Sons Ltd., Chichester, 2008.
4. P. G. Ciarlet and J.-L. Lions (eds.), Handbook of numerical analysis. Vol. II, Handbook of Numerical Analysis, II, North-Holland, Amsterdam, 1991, Finite element methods. Part 1.
5. Philippe G. Ciarlet, The finite element method for elliptic problems, Classics in Applied Mathematics, vol. 40, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002, Reprint of the 1978 original [North-Holland, Amsterdam].
6. E. Hairer, S. P. Norsett, and G. Wanner, Solving ordinary differential equations. I, second ed., Springer Series in Computational Mathematics, vol. 8, Springer-Verlag, Berlin, 1993, Nonstiff problems.
7. E. Hairer and G. Wanner, Solving ordinary differential equations. II, second ed., Springer Series in Computational Mathematics, vol. 14, Springer-Verlag, Berlin, 1996, Stiff and differential-algebraic problems.
8. Alfio Quarteroni and Alberto Valli, Numerical approximation of partial differential equations, Springer Series in Computational Mathematics, vol. 23, Springer-Verlag, Berlin, 1994. (FD schemes and FE for PDEs) One can download a pdf file from MIMUW computers (valid Dec 2014): Springer Link

LAB

Labs

• Lab 1 Introduction to octave. Euler's schemes
  • 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?
    m-file with explicit Euler scheme and a script with some tests:
    exEuler.m - explicit Euler scheme (in many dimensions)
    testexEul.m - basic tests of explicit Euler scheme (in 1D and 2 D)
• Lab 2 - Function lsode() in octave and simple ODE schemes cont. In multistep ODE schemes take x_1,x_2 etc as exact solutions if not take x_k k=1,..,p- p-1 computed by using some explicite 1-step method of the same order e.g. Taylor scheme for midpoint one. midpoint - midpoint unstable for u'=-u u(0)=1 for long time
  • 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.
  
  OrderExEuler.m - tests of local truncation error order of explicite Euler scheme
  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)
  taylor.m - explicit Taylor scheme (in many dimensions)
  testtaylor.m - basic tests of Taylor scheme (in 2 D)
• Lab 3 Experimental testing the order of ODE schemes cont. Testing starting for multistep schemes (for Adams-bashforth of 2nd order). the shooting method
  • 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 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.
  • 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.
  
  AdamsB2.m - explicit Adams-Bashforth scheme (of order 2)
  testAB2start.m - tests starting value (x^h_1) for explicit Adams-Bashforth scheme (of order 2)
  testAdamsBo2.m - tests for explicit Adams-Bashforth scheme (of order 2) like for midpoint - i.e. error order etc
  testshoot.m -shooting method for linear boundary value ODE problem: -y''+c(x)y=0; y(0)=1 y(b)=1 (b=1 - shooting works fine, b=20 - shooting doees not work at all - why?)
  linshooting.m -m-file with a function linshooting() solving the linear boundary value ODE problem: -y''+d(x)y+c(x)y=f(x); y(a)=ya y(b)=yb using shooting method
  shooting.m -m-file with a function shooting() solving the boundary value ODE problem: y''=F(x,y,y'); y(a)=ya y(b)=yb using shooting method
• Lab 4 - Finite difference method (FDM) for -u''+cu=f with Dirichlet bnd cond : u(a)=alpha u(b)=beta and for mixed bnf conditions: u'(a)=alpha u(b)=beta
  • Shooting - method- problems from the previous lab
  • 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''=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} .
  • 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\}$
  
  FDsolver.m function solving -u''+cu=f with Dirichlet bnd cond : u(a)=alpha u(b)=beta by FDM method (order two)
  testFD.m tests of FDM method for solving -u''+cu=f with Dirichlet bnd cond : u(a)=alpha u(b)=beta -prob 1 u''=u u(0)=u(b)=1 b=1,4,8,16 etc plots of graphs; problem 2: order of FDM with 3 points stencil for -u''=sin(x) u(0)=0 u(b)=sin(b); prob 3 - same as prob 2 but with u_h(b)=solution(b+h/2) - order 1.
  FDmixlft.m function solving -u''+cu=f with mixed bnd cond : u'(a)=alpha u(b)=beta by FDM method (order one)
  FDtest.m - tests of order of the both methods for a few model examples
• Lab 5 FDM for Poisson equation in 2D
  • 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 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 e.g. discontinuous e.g. $c(x,y)=1 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 the known solution $u=u1^p$ for $u1(x)=\max(0,\sin(\pi*x)\sin(\pi*y))$ and p=1,2,3,4 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.
  • 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\}$
  
  Lap2Dsq.m - octave m-file with function creating sparse matrix FD -Laplacian on a unit square on uniform mesh (5 point stencil)
  FDM2dtest.m - octave script testing order of convergence of FDM for -Laplacian u=f on a 2D unit square; u=g on boundary; with uniform mesh (5 point stencil) in max and L^2 discrete norms; for u solution = sin(pi*x_1)*sin(pi*x_2), (then f=2*(pi)^2*u; u=0 on boundary)
• Lab 6- FEM in 1d Finite element method (FDM) for $-u'' +cu=f$ on $[a,b]$ with Dirichlet bnd condition. In all problems we use linear continuous finite element method i.e conforming $P_1$ FEM.
  • 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 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/2]$ 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; Dirichlet non-homogeneous bnd conditions Find FEM $P_1$ spproximation of $-u'' +cu=f$ $u(a)=0=u(b)$ on a regular grid for $u=\sin(\pi*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 method. Compare with FDM method.
  • 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.
  • FEM linear 1D; general elliptic operator regular grid; Dirichlet homogeneous bnd conditions Find FEM $P_1$ spproximation 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$.
  • 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.
  
  FEM1Dsolver.m - - linear 1D FEM solver for -au''+bu'+cu=f with Dirichlet or Robin bc - na siatce dowolnej
  testMESDbc.m - a m-file with a function testing 1D linear FEM with Dirichelt bc
  testMESor.m - script with few calls to testMESDbc() (in m-file testMESDbc.m -which should be downloaded) - a few tests of linear 1D FEM order by the halving step method.

Octave scripts with solution of some problems from our labs