BDF 3

We use $utrue=u(t,x,y)=r(t)\cos(2\, x)\sin(4\, y)$ to test BDF(3) 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$ BDF(3) scheme for $\partial_t u-F (u) = 0$ is:

Given $u^0,\;u^1,\,u^2$ compute $$\begin{array}{rcl} u^{n+1}-\frac{6}{11}\Delta t\,F^{n+1}& =& \frac{18}{11}u^n-\frac{9}{11}u^{n-1}+\frac{2}{11}u^{n-2}\quad n=2,3,... \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,\;u^1,\,u^2$ $$\begin{array}{rcl} \frac{1}{\Delta t}(u^{n+1},v)+\frac{6}{11}(\nabla u^{n+1},\nabla v)-\frac{6}{11}(g(u^{n+1}),v)& =&\frac{1}{\Delta t} (\frac{18}{11}u^n-\frac{9}{11}u^{n-1}+\frac{2}{11}u^{n-2},v)+ \frac{6}{11}(f^{n+1},v),\quad n=2,3,... \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.10.002498280.98048
0.050.0001585810.248803
0.0259.69022e-060.0620477

verbosity=0;
//
// The BDF3 scheme is used to solve the heat equation
//
//  Aplication of BDF3 to u'(t)-[Lap(u)+g(u)+fh(t)]=0 gives
//
// BDF(3)---> u_{n+1}-(6/11)*dt*[Lap(u_{n+1})+g(u_{n+1})+fh_{n+1}]=
//                                                     (18/11)*u_n-(9/11)*u_{n-1}+(2/11)*u_{n-2}
//
//
// Variational formulation (in the case boundary integrals are cero)
//
//                                     ((18/11)*u_n-(9/11)*u_{n-1}+(2/11)*u_{n-2},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
//            -(6/11)*dt*((p^2+q^2)*u_{n+1},v)
//            -(6/11)*dt*(fh,v)=((18/11)*u_n-(9/11)*u_{n-1}+(2/11)*u_{n-2},v)
//
//
// solution parameters
//
int p=2,q=3;
real t=0;

// Domain square
//
int factor=4;
int nn=10*factor;
mesh Th=square(nn,nn);
//
fespace Xh(Th,P2);
Xh uh,vh,uhrhs,w;
real Tfin=5, t0=0.;
real dt=1.0/nn;

real r,drdt,errL2sq=0.0,errH1sq=0.0,errLiL2=0.0,errL2H1=0.0,locL2=0.0,locH1=0.0;
Xh Phi=cos(p*x)*sin(q*y);

t=t0;
r=t^3-50;// r(t) function
Xh uhold0=r*Phi;// finite element function to keep previous value (set at time t0)
t=t0+dt;
r=t^3-50;// r(t) function
Xh uhold1=r*Phi;// finite element function to keep previous value (set at time t0+dt)
t=t0+2*dt;
r=t^3-50;// r(t) function
Xh uhold2=r*Phi;// finite element function to keep previous value (set at time t0+dt)

Xh fh;
Xh bdata;

//
// Main problem definition

problem parabolicbdf3(uh,vh) =int2d(Th)(uh*vh/dt-(6.0/11.0)*(p^2+q^2)*uh*vh+
(6.0/11.0)*(dx(uh)*dx(vh)+dy(uh)*dy(vh))
)
- int2d(Th)((6.0/11.0)*fh*vh
+((18.0/11)*uhold2-(9.0/11)*uhold1+(2.0/11)*uhold0)*vh/dt)
+on(1,2,3,4,uh=bdata);

//
//Time loop
//
for (t=t0+3*dt;t<=Tfin;t+=dt)
{

//
// 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;

parabolicbdf3;

//
Xh uexh=r*Phi;
Xh dxuexh=-p*r*sin(p*x)*sin(q*y);
Xh dyuexh=q*r*cos(p*x)*cos(q*y);

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=uhold1;
uhold1=uhold2;
uhold2=uh;
}

//Errors in time
errLiL2 = sqrt(errL2sq);
errL2H1 = sqrt(errH1sq);
cout << "Using BDF(3) 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;