We use \(utrue=u(t,x,y)=r(t)\cos(2\, x)\sin(4\, y)\) to test BDF(1), or implicit Euler scheme, for the heat operator. The problem is
$$ \begin{array}{rcll} \partial_t u-\Delta u +g(u)& =& f& \mbox{on } \Omega=(0,1)^2,\quad t>0\\ u(t,x,y)&=& utrue &\mbox{on } \partial\Omega,\quad t>0\\ u(t=0)&=&utrue(t=0)& \mbox{on } \Omega \end{array} $$
Using \(F(t,u)=\Delta u(t) -g(u)+f(t)\) and \(F^n=\Delta u^n-g(u^n)+f^n\) BDF1 scheme for \( \partial_t u-F (u) = 0\) is
Given \(u^0\) compute $$ \begin{array}{rcl} u^{n+1}-\Delta t\,F^{n+1}& =& u^n,\quad n=0,1,2,... \end{array} $$
and the variational formulation looks for \(u^{n+1}\in H^1(\Omega)\) such that \(u=utrue\) on \(\partial\Omega\) and for all \(v\in H^1_0(\Omega)\)
Given \(u^0\) $$ \begin{array}{rcl} \frac{1}{\Delta t}(u^{n+1},v)+(\nabla u^{n+1},\nabla v)-(g(u^{n+1}),v)& =&\frac{1}{\Delta t} (u^n,v)+ (f^{n+1},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=5\) the results are
\(\Delta t=h\) | \(L^\infty(L^2)\) | \(L^\infty(H^1)\) |
---|---|---|
0.1 | 1.56634 | 7.0906 |
0.05 | 0.75159 | 3.23678 |
0.025 | 0.367014 | 1.55258 |
download this example: BDF1_Heat.edp or return to 2D examples
verbosity=0; // // The BDF1 scheme is used to solve the heat equation // // Aplication of BDF1 to u'(t)-[Lap(u)+g(u)+fh(t)]=0 gives // // BDF(1)---> u_{n+1}-dt*[Lap(u_{n+1})+g(u_{n+1})+fh_{n+1}]=u_{n} // // BDF(1) is just backward Euler // // Variational formulation (in the case boundary integrals are cero) // // (u_{n+1},v)+dt*(Grad(u_{n+1}),Grad(v))-dt*(g(u_{n+1})+fh_{n+1},v)=(u_{n},v) // // Use of exact solution u=r(t)*Phi(x,y) with Phi(x,y)=cos(p*x)*sin(q*y) // // satisfies u_t-[Lap(u)+(p^2+q^2)*u+r'(t)*Phi(x,y)]=0 // // because u_t=r'(t)*Phi(x,y) // Lap(u)=r(t)*Phi(x,y)*(-p^2-q^2)=-(p^2+q^2)*u // // using fh=r'(t)*Phi(x,y) // // Variational formulation is // // (u_{n+1},v)+dt*(Grad(u_{n+1}),Grad(v))-dt*((p^2+q^2)*u_{n+1},v)-dt*(fh,v)=(u_{n},v) // // // solution parameters // int p=2,q=4; real t=0; // // Domain square // int factor=4; int nn=10*factor; mesh Th=square(nn,nn); real dt=1.0/nn;// time step fespace Xh(Th,P2); Xh uh,vh; real Tfin=5, t0=0.; real r,drdt,errL2sq=0.0,errH1sq=0.0,errLiL2=0.0,errL2H1=0.0,locL2=0.0,locH1=0.0; int i; Xh Phi=cos(p*x)*sin(q*y); t=0;// initial time r=t^3-50;// r(t) function Xh uhold0=r*Phi;// finite element function to keep previous value (set at time =0) Xh fh; Xh bdata; // // Main problem definition // problem parabolicbdf1(uh,vh) =int2d(Th)(uh*vh/dt-(p^2+q^2)*uh*vh+ (dx(uh)*dx(vh)+dy(uh)*dy(vh)) ) - int2d(Th)(fh*vh+uhold0*vh/dt) +on(1,2,3,4,uh=bdata); // //Time loop // for (t=t0+dt;t<=Tfin;t+=dt) { //cout << "Time : " << t << endl; // // solving one step of the problem // r=t^3-50;// r(t) function drdt=3*t^2;// r'(t) function fh=drdt*Phi; bdata=r*Phi; parabolicbdf1; // Xh uexh=r*Phi; Xh dxuexh=-p*r*sin(p*x)*sin(q*y); Xh dyuexh=q*r*cos(p*x)*cos(q*y); Xh err=uh-uexh; plot(cmm="Time level "+t,uh,dim=3,fill=1,wait=0,value=1); //Error computing locL2=int2d(Th)((uh-uexh)^2); locH1=locL2+int2d(Th)((dx(uh)-dxuexh)^2)+int2d(Th)((dy(uh)-dyuexh)^2); errL2sq=max(errL2sq,locL2); errH1sq=errH1sq+dt*locH1; // // update the time level // uhold0=uh; } //Errors in time errLiL2 = sqrt(errL2sq); errL2H1 = sqrt(errH1sq); cout << "Using BDF(1) and (p,q) = " << p<<" "<<q<<endl; cout << "Errors at time T = " << Tfin<<" using dt=dh = " << dt << endl; cout << "Errors : Linfty(L2) " << errLiL2 << " L2(H1) " << errL2H1 << endl;