We use \(u=u(t,x,y)=r(t)\sin(2\pi x)\cos(\pi y)\) to test Crank-Nicolson scheme for the heat operator. We choose the function \(f(t,x,y)\) such that the problem is
$$ \begin{array}{rcll} \partial_t u-\Delta u & =& f& \mbox{on } \Omega=(0,1)\times (-0.5,0.5),\quad t>0\\ u(t,x,y)&=& 0 &\mbox{on } \partial\Omega,\quad t>0\\ u(t=0)&=&u_0(x,y)& \mbox{on } \Omega \end{array} $$
Using \(F(t,u)=\Delta u(t) +f(t)\) and \(F^n=\Delta u^n+f^n\) Crank-Nicolson scheme for \( \partial_t u-F (u) = 0\) is
Given \(u^0\) compute $$ \begin{array}{rcl} u^{n+1}-\frac{\Delta t}{2}F^{n+1}& =& u^n+\frac{\Delta t}{2}F^{n},\quad n=0,1,2,... \end{array} $$
and the variational formulation looks for \(u^{n+1}\in H^1_0(\Omega)\) such that for all \(v\in H^1_0(\Omega)\)
Given \(u^0\) $$ \begin{array}{rcl} (u^{n+1},v)+\frac{\Delta t}{2}(\nabla u^{n+1},\nabla v)& =& (u^n,v)-\frac{\Delta t}{2}(\nabla u^{n}, \nabla v)+\frac{\Delta t}{2}(f^{n+1}+f^n,v),\quad n=0,1,2,... \end{array} $$
Using \(P_2\) finite elements and computing the \(L^\infty(0,T;L^2(\Omega))\) and \(L^\infty(0,T;H^1(\Omega))\) errors, for \(T=3\) the results are
| \(\Delta t=h\) | \(L^\infty(L^2)\) | \(L^\infty(H^1)\) |
|---|---|---|
| 0.1 | 0.0469857 | 0.615365 |
| 0.05 | 0.0105996 | 0.15044 |
| 0.025 | 0.00260145 | 0.0374534 |
download this example: Crank_Nicolson_Heat.edp or return to 2D examples
//
// Testing Crank-Nicolson for heat equation
//
// Using of keyword:problem
//
//----------------------------------------------------------------
// Definition of the domain: D={square [0,1]x[-0.5,0.5]}
//
// then u(t,x,y)=g(t)*sin(2*pi*x)cos(pi*y) is cero on boundary of D all time
//
// u_t(t,x,y)=dgdt(t)*sin(2*pi*x)cos(pi*y)
// u_{xx}(t,x,y)=-4*pi^2*u(t,x,y)
// u_{yy}(t,x,y)=-pi^2*u(t,x,y)
//
// As a consequence
//
// u_t-Lap(u)=dgdt(t)*sin(2*pi*x)*cos(pi*y)+5*pi^2*g(t)*sin(2*pi*x)*cos(pi*y)
//
//
verbosity=0;
func real g(real t){
real s= 5*cos(10*t);//-t^3-t^2+1;
return s;
}
func real dgdt(real t){
real s=-50*sin(10*t);// -3*t^2-2*t;
return s;
}
int factor=1;
int m=10*factor;// number of points
mesh th=square(m,m,[x,-.5+y]);
real h=1.0/m;// size of h
real dt=h;// size of dt
//
// Definition of the axis OX-OY just for fun
//
border OX(t=-0.5,1){x=t;y=0;}
border OY(t=-0.5,1){x=0;y=t;}
// plot initial mesh
plot(th,OX(1),OY(1),wait=0);
//
// Definition of the fespace
//
fespace Vh(th,P2);
Vh u,v;
Vh u0=g(0)*sin(2*pi*x)*cos(pi*y);// initial data
Vh fh1,fh0;// force terms at time levels t^{n+1} and t^n respectively
//
// plot initial data
//
plot(u0,OX(1),OY(1),fill=1,dim=3,wait=0);
//
// Heat equation u_t-Lap(u)=f solved by Crank-Nicolson
// u(t=0)=u0
// u=0 on boundary for all t
//
// u^n means u at time level t^n
//
// (u^{n+1}-u^n)/dt-(Lap(u^{n+1})+Lap(u^n))/2=(f^{n+1}+f^{n})/2
//
// or
// u^{n+1}-(dt/2)*Lap(u^{n+1})=u^n+(dt/2)*Lap(u^{n})
// +(dt/2)*(f^{n+1}+f^{n})
// Variational formulation
// (u,v) means integral on computational domain
//
// (u^{n+1},v)+(dt/2)*(Grad (u^{n+1}),Grad(v))=
// (u^n,v)-(dt/2)*(Grad(u^{n}),Grad(v))
// +(dt/2)*((f^{n+1}+f^{n}),v)
problem heat(u,v)=int2d(th)(
u*v+0.5*dt*(dx(u)*dx(v)+dy(u)*dy(v))
)
-int2d(th)(
u0*v-0.5*dt*(dx(u0)*dx(v)+dy(u0)*dy(v))
)
-int2d(th)(0.5*dt*(fh1+fh0)*v
)
+on(1,2,3,4,u=0);
real t=0;// initial time
real tf=3;// final time
real errL2sq=0,errH1sq=0;// to hold errors
for(int n=0;n<tf/dt;n++)
{ fh0=(dgdt(t)+5*pi^2*g(t))*sin(2*pi*x)*cos(pi*y);
t=t+dt;
fh1=(dgdt(t)+5*pi^2*g(t))*sin(2*pi*x)*cos(pi*y);
heat;
Vh uexh=g(t)*sin(2*pi*x)*cos(pi*y);
Vh dxuexh=2*pi*g(t)*cos(2*pi*x)*cos(pi*y);
Vh dyuexh=-pi*g(t)*sin(2*pi*x)*sin(pi*y);
//Error computing
real locL2=int2d(th)((u-uexh)^2);
real locH1=locL2+int2d(th)((dx(u)-dxuexh)^2)+int2d(th)((dy(u)-dyuexh)^2);
errL2sq=max(errL2sq,locL2);
errH1sq=errH1sq+dt*locH1;
plot(cmm="Approximated solution at time "+t,u,
OX(1),OY(1),fill=1,dim=3,wait=0,value=1);
//plot(cmm="Approximated solution at time "+t,errf,fill=1,dim=3,wait=0,value=1);
u0=u;
}
//Errors in time
real errLiL2 = sqrt(errL2sq);
real errL2H1 = sqrt(errH1sq);
cout << "Using Crank-Nicolson " <<endl;
cout << "Errors at time T = " << tf<<" using dt=dh = " << dt << endl;
cout << "Errors : Linfty(L2) " << errLiL2 << " L2(H1) " << errL2H1 << endl;