I compute the solution for several spatial mesh sizes and measure the L2 error with respect to the exact solution. My goal is to study the order of convergence, but the computed rates seem inconsistent or lower than expected.
verbosity=0;
macro dn(u) (N.xdx(u)+N.ydy(u)) // def the normal derivative
// parameter
real dt, h, T= 1;
int nref =4;
real [int] L2error ( nref ); // initialize the L2 error array
real [int] Dx( nref ); // initialize the Space discretization array
real [int] DT( nref ); // initialize the Time discretization array
real L2err;
real eorr;
macro Grad(u)[dx(u),dy(u)]//
macro uex(t)(exp(-t)(x^2 - y^2)) //
macro g(t) (exp(-3(t))*(x^2 - y^2)^3) //
//Mesh
for ( int n=0; n< nref ;n++) {
int M =52^(n) ;
h =1./ M;
Dx[n]=h;
dt=1./256;//0.4h;//0.4h^2;//
DT[n]= dt;
// border bottom(t=0,1){x = 0 + t2 - 1; y = -1; label = 1;}
// border right(t=0,1){x = 1; y = -1 + t2; label = 1;}
// border top(t=0,1){x = 1 - t2; y = 1; label = 1;}
// border left(t=0,1){x = -1; y = 1 - t*2; label = 1;}
mesh Th = square(M,M,[2x-1,2y-1]);//buildmesh(bottom(M) + right(M) + top(M) + left(M));
cout << "level M = " << M << endl;
//plot(Th, wait=false);
//Fespace
fespace Vh(Th,P1);
fespace Vh2(Th,P1);
//fespace Vh2(Th,[P1dc,P1dc]);
//fespace Vh2(Th,[P1,P1]);
// Function
func real df(real u){return u^3; }
func real ddf(real u){return 3*u^2; }
Vh u, uold, uoldd, auxv, vv;
Vh2 dfalpha;
Vh2 ddfalpha;
uoldd = x^2 - y^2;
uold = uoldd - dt*(x^2 - y^2);
plot(u,wait=0,cmm=“Klein-Gordan Exact Solution” ,value=true);
for (real t=0.;t<T;t+=dt){
uold=uoldd;
u=uold;
real res;
int iter=0;
for (int k=0;k<50;k++) // Newton loop
{
dfalpha=df(u);
ddfalpha= ddf(u);
varf Fv(v, phi)=int2d(Th)
((u-2uold+uoldd)phi1./(dtdt)+(dx(u)*dx(phi)+dy(u)*dy(phi))1./(2dt)
-(dx(uoldd)*dx(phi)+dy(uoldd)*dy(phi))1./(2dt)+(dx(u)*dx(phi)+dy(u)*dy(phi))1./2+(dx(uoldd)dx(phi)+dy(uoldd)dy(phi))1./2+ (u-uoldd)phi1./(2dt) + (dfalpha)phi1./4 +((uoldduoldd)u)phi1./4 +((uu)uoldd)phi1./4+((uoldduoldduoldd)phi1.4))// dfalphaphi)
+int2d(Th)(-g(t)*phi);// + on(1,2,3,4, u=0);
auxv=Fv(0,Vh); // nabla J
res= sqrt(auxv'auxv[]); // Residual
// cout << " residu = " << res << endl;
if (res< 1e-6) break;
varf dFv(v, phi) = int2d(Th)(
vphi*(1./(dt*dt))
- (dx(v)dx(phi) + dy(v)dy(phi))(1./(2dt) + 1./2)
- vphi(1./(2*dt))
- ddfalphavphi // df/du = 3u^2 * v
+(uvuoldd)phi1./2 - (1./4)uoldd^2v*phi // df/du = 3u^2 * v
);
// varf dFv(v,phi)=int2d(Th)
// (vphi1./(dtdt) + (dx(v)dx(phi)+dy(v)dy(phi))1./(2dt) +(dx(v)dx(phi)+dy(v)dy(phi))1./(2)+ vphi1./(2dt)+((vv)phi3./4)+((uoldduoldd)phi1./4)+((uuoldd)phi1./2));
matrix H=dFv(Vh,Vh, factorize=1,solver=LU); // HJ
vv=H^-1*auxv;
u-=vv;
iter++;
}// end of Newton loop
// uold=2.*u-uold;
plot(u,wait=0,cmm=“Klein-Gordan Approximate Solution”, value=true);
eorr= abs(u-uex(t+dt))^2;
L2err=sqrt(int2d(Th)(abs(u-uex(t+dt))^2));
cout << "time = " << t+dt << " iter = " << iter << " residu = " << res << " L2 error = " << L2err << " eror = "<< eorr << endl;
}//end of time loop
cout << “final time T is reached” << endl;
L2error [n]= L2err;
}//end of loop on n< nref
for ( int n =0;n< nref ;n ++)
cout << " L2error " << n << " = " << L2error [n] <<endl;
for ( int n =1; n< nref; n ++){
cout <<" Space convergence rate = " << log ( L2error [n-1]/L2error
[n])/log(Dx[n-1]/Dx[n]) <<endl;
}
