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).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?y(n+1)=y(n)+ 0.5h*(f(n)+f(t(n+1),y(n)+hf(n))
and modified Eulera schemesy(n+1)=y(n)+ h*f(t(n)+0.5h,y(n)+0.5hf(n))
on the IVP:y'=-sin(y), y(0)=0, y'(0)=1.
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.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))
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.y'=ay y(0)=1; (penduulum equations) y''=-y (linear) or y''=-sin(y) with y(0)=1 y'(0)=0.
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 getx(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.-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
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).-y''+y=0; y(0)=1; y(b)=1 for b=1,4,10,15,20.
Draw graphs. Compute errors |y(b;s)-yb|.-y''+y=f z y(0)=0; y(b)=0 for f=2*sin(x) and b=\pi (the solution can be computed).
-y''+y=-2+(1-x*x); y(-1)=0; y(1)=0. (the solution can be computed).
-y''+exp(-x) y =cos(x); y(-1)=0; y(1)=0.
-y''+ay'+ y =cos(x); y(-1)=0; y(1)=0 for a=0,1,10,100,-1,-10,-100.
-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)).
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
-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.-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.-u''=\sin(x) , u(0)=u(\pi)=0
in L^2 and max discrete norms.-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.-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.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 havel1(t)=u_{k-1}(t-x_k)/(-h)+u_k(t-x_{k-1})/h.
- 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.- 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 approximateu'(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.-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.-Laplacian u=f in [0,1]^2 , u=0 on the boundary
for difffrent N (or h=1/N).-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} .-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.-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
- in nodal basis we get the system (A+cB)u=F with A, B (or their submatrices...) obtained from the script below, the rhs for linear FEM we can approximate by composite trapezoidal rule i.e. F=(f_i)_i with f_i=\int_a^b f\phi_i dx \approx 0.5*f(x_i)(x_{i-1}-x_{i+1}) i=1,...,N-1 (note that in case of equidistant mesh x_k=a+k*h we get f_i=h*f(x_i) as in FDM method...). The H^1 seminorm and L^2 norm may be computed for u\inV_h as |u|_1^2=u^TAu and ||u||_0^2=u^TMu - u vector of coefficients in nodal basis.| \int_a^b \phi_k\phi_l dx= | 0 | |k-l|>1; |
|---|---|---|
| (1/3)*(|x_k-x_{k+1}|+|x_k-x_{k-1}|) | k=l ; | |
| (1/6)*|x_k-x_l| | |k-l|=1; |
| \int_a^b phi_k'phi_l'dx= | 0 | |k-l|>1; |
|---|---|---|
| (|x_k-x_{k+1}|^{-1}+|x_k-x_{k-1}|^{-1}) | k=l; | |
| -|x_k-x_l|^{-1} | |k-l|=1 |
-au''+cu=f
i.e.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
for \phi_k nodal basis - piecewise linear - any grid - to get matrix for zero Dirichlet bnd - take submatrices of A and B The returned matrices are for the nodal basis (\phi_i)_{i=0}^N i.e. includes the basis functions related to bnd poits, in case of Dirichlet bnd we have to take submatrices obtained by removing first and last rows/cols.
-au''+cu=f in (a,b)
u(a)=al, u(b)=be
u_t=u_{xx} x\in(0,pi), t\in(0,T] u(t,0)=u(t,pi)=0, u0(x)=sin(x),
for known solution u(t,x)=exp(-t)sin(x). Use lsode() for ODE solver for FDm semi-discrete problem: U'=-T U where T FDM discretization of - 2nd derivative i.e. 3diagonal matrix (N-1)x(N-1) with 2/h^2 on the main diagonal, and -1/h^2 od sub and super diagonals, here h=(b-a)/N, U=(u(1),..,u(N-1))^T vector of values at discrete points.u_t=u_{xx} u(t,a)=u(t,b)=0
with discontinuus u0 e.g.on [0,2] u0=0 on [0,1] and u0(x)=1 on [1,2]. Plot the FDM solution for t=10,1,1e-1,1e-2 for h=1/10,1/100,1/1000.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)
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
-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
- Laplacian u=2pi^2*sin(x)sin(y)
u=0 on boundary
-au''+cu=f
i.e.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
for \phi_k nodal basis - piecewise linear - any grid - to get matrix for zero Dirichlet bnd - take submatrices of A and B The returned matrices are for the nodal basis (\phi_i)_{i=0}^N i.e. includes the basis functions related to bnd poits, in case of Dirichlet bnd we have to take submatrices obtained by removing first and last rows/cols.
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)