Numerical Differential Equations

winter semester 2015-16

Exam questions

1. One step schemes. Examples and convergence theory
2. Multistep schemes . Examples including Adams schemes. Convergence theory.
3. Stiffness. Definition. Examples of ODEs. Idea of adaptive step control.
4. Finite difference method for elliptic equation in 1D and 2D. 1D Example: - u''+u=f u(a)=u(b)=0 and detailed convergence analysis of this problem in a discrete maximum norm
5. Finite difference method for elliptic equations. 2D Example: -Laplacian u=f , u=0 on boundary. Abstract convergence anaylisi: the order of local truncation error and stability. Discrete convergence - definition and Lax theorem.
6. Idea of Finite Element Method. Analysis of convergence of linear finite element in 1D for -u''=f, u(a)=u(b)=0.
7. Elements of abstract FEM theory: Lax Milgram theorem. Cea Lemma. Application to 2d elliptic boundary problem with homogeneous Dirichlet boundary element.
8. FEM and different boundary conditions for 2nd order elliptic boundary problem: Dirichelt, Neumann, Robin, mixed. (it is not discussed in Lect Notes. I discussed it during a lecture but I will not ask this question unless somebody choose it)
9. Elements of abstract FEM theory: continuous FEM spaces, affine families of FEM spaces, shape regularity etc
10. Schemes for parabolic differential equations: FDM and FEM in 1/2D (as in lecture notes or as was presented in the lecture)
11. Mixed formulation of elliptic 2nd order problem and its mixed FEM discretization. Stokes problem in a mixed formulation. (it is not discussed in Lect Notes. I will present it during the last lecture but I will not ask this question unless somebody choose it)

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?
  • Solve the same problem 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.
  
  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 (Nov 3, 2015) 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.
  • 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
  
  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, local truncation error order etc
  testHeun.m - tests for Heun scheme (explicit Runge-Kutta scheme of order 2)
• Lab 4 and 5 (Nov 24, 2015 and Dec ?, 2015)- Finite difference method (FDM) for -u''+bu'+ cu=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, 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} .
  • 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\}
  
  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 does not work at all - why?)
  linshoot15.m -m-file with a function linshoot15() solving the linear boundary value ODE problem: -y''+p(x)y+q(x)y=f(x); y(a)=alpha y(b)=beta 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
  fdmsolve15.m function solving -u''+c(t)u=f(t) with Dirichlet bnd cond : u(a)=alpha u(b)=beta by FDM method (order two)
  testfdm1D15.m function testing the order of convergence of -u''=sin(t) with Dirichlet bnd cond : u(a)=alpha u(b)=beta (solution u=sin(t)) by FDM method (order two) and the order of local truncation error
• Lab 6 and 7 FDM for Poisson equation in 2D (Dec 15 and 22 2015)
  • 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...
  
  test2dlapFDM15.m - octave script testing order of convergence of FDM for -Laplacian u=f on a 2D unit square; u=0 on boundary; with uniform mesh (5 point stencil) in max discrete norm; for u solution = sin(x_1)*sin(x_2) + p1, (then f=2*u; u=p1 linear poly on boundary)
• Lab 8- 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.
• Lab 9 -- FD method for parabolic equations in 1D i.e. we discretized equation u_t-u_{xx}=f by FDM with respect to x - and apply octave ODE solver (lsode()) to the resulting ODEs system
  • Run the scripts - i.e. solve u_t=u_{xx} for diccontinuous u0
  • Change the script in order to solve u_t-u_{xx}=f u(t,a)=ga(a) u(t,b)=gb(b) 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
  • 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)
  
  FD1dtest.m - octave script with code testing convergence order for FDM discretization of u_t-u_{xx}=0 u(0)=u(pi)=0 u(0,x)=sin(x) - u(t,x)=exp(-t)sin(x) is the solution (the ODE solver is lsode() - octave black box for ODEs)
  FEM1dtest.m - octave script with code testing convergence order for linear FEM discretization of u_t-u_{xx}=0 u(0)=u(pi)=0 u(0,x)=sin(x) - u(t,x)=exp(-t)sin(x) is the solution -one can use any 1D traingulation - please test a few different (the ODE solver is lsode() - octave black box for ODEs)
  FDstep.m - octave script with code testing diffusion in u_t-u_{xx}=0 u(0)=u(pi)=0 u(0,x)= indicator function of [1,2.5] and tests with same u0 but with force term f<>0 i.e. u_t-u_{xx}=f
  nrr15-parab1d.tgz - octave scripts and functions with code testing the order of convergence 3 basic schemes (ex/imlicit Euler, Crank-Nicholson) for u_t-u_{xx}=0 u(a)=ae u(b)=be u(0,x)= u0(x) with FDM or FEM space discretizations; please read README file in the archive

Octave scripts with solution of some problems from our labs