I have shared the code here and screenshot of the pde:Please ignore the error calculating part in the code since I am not doing it for known exact solution:
// A step of the midpoint method can be carried out by first taking
// a backward Euler step of size dt/2, followed by a forward Euler
// step of size dt/2.
//
//constant isotherm
// Location:
//
// https://people.sc.fsu.edu/~jburkardt/freefem_src/midpoint/midpoint.edp
cout << endl;
cout << “Midpoint:” << endl;
cout << " FreeFem++ version." << endl;
cout << " Solve a time-dependent boundary value problem in a square." << endl;
cout << " Use the midpoint method to handle the time variation." << endl;
real omega= 0.5;
//real rho=1.0;
real m;
real a;
real b;
real c;
real h;
real tm;
real hold;
//real epsr=1e-10;
real totalmid;
//int nt = 100;
real t;
real tmin = 0.0;
//real tmax =0.01 ;
real tmax = 100.0;
//not sure what to choose for tmax in our problem
real told;
ofstream fout(“concentration_mid”);
//real dt = ( tmax - tmin ) / nt; we want dt to be h^2/16, h^2/256 according to thesis which I defined later after writing h.
// Express the number of nodes as (almost) a power of 2.
//for ( int nlog2 = 1; nlog2 <= 5; nlog2++ )
//skipping for loop for now to see atleast for first h
// h values given in thesis is 1/8,1/16,1/32,1/64,1/128
{
int nlog2 = 3;
//int nlog2 = 3;
//int nlog2 = 4;
//int nlog2 = 5;
//int nlog2 = 6;
//int nlog2 = 7;
// n: number of nodes in each spatial direction.
//
int n = pow ( 2.0, nlog2 );
//
// h: spacing between nodes.
//0.01
h = 1.0 / n;
int nt = 32;
real dt = (tmax-tmin) / nt;
// we want dt to be 1/32,1/64,1/128,1/256,1/512 according to thesis
//int nt = 32;
//real dt = (tmax-tmin) / nt;
//int nt = 64;
//real dt = (tmax-tmin) / nt;
//int nt = 128;
// real dt = (tmax-tmin) / nt;
//int nt = 256;
//real dt = (tmax-tmin) / nt;
// int nt = 512;
// real dt = (tmax-tmin) / nt;
cout << " H = " << h << endl;
cout << " DT = " << dt << endl;
//
// Th: the triangulation of the square.
real x0 = 0;
real x1 = 3;
real y0 = 0;
real y1 = 20;
//int n = 5;
//real m = 20;
mesh Th = square(n, n, [x0+(x1-x0)x, y0+(y1-y0)y]);
//
//mesh Th = square ( n, n );
//
// Define Vh, the finite element space. “P1” = piecewise linears.
//
fespace Vh ( Th, P1 );
//
// ch: the finite element function that will approximate the solution.
// coldh: the solution at the previous time step.
//
Vh ch;
Vh cmh;
Vh coldh;
//Vh qh;// probably i don’t need it. because everywhere i used qh=1+ch directly
//Vh qoldh;//probably i don’t need it. because everywhere i used qoldh=1+coldh directly
//
// vh: the test functions used to project the error.
//
Vh vh;
// Cexact: the exact solution.
//
// func Cexact = t^2( x^3 -1.5x^2+1.0 ) * cos ( (pi/4) * y );
// func q = constant;
//
//
// f: the right hand side.
//
func f = 0;
func fm = 0;
problem halfbackwardeulerstep ( cmh,vh ) =
int2d ( Th )
(
omega* cmh * vh
)
- int2d ( Th )
(
omega* coldh * vh
)
+ int2d ( Th )
(
0.5*dt*( (0.0*(dx(cmh))) + (0.1*(dy(cmh))) ) * vh
+ 0.5*dt*(dx(cmh) * dx(vh) + dy(cmh) * dy(vh) )
)
- int2d ( Th )
(
0.5*dt*fm* vh
)
+ on ( 3, cmh = 1 );
cout << endl;
// cout << " T ||cerroroo0 ||cerror00|| ||gradcerror00|| ||crror01||" << endl;
cout << endl;
real cerroroo0 = 0.0;
real cerror00 = 0.0;
real gradcerror00 = 0.0;
real cerror01 = 0.0;
for ( int it = 0; it <= nt; it++ )
{
if ( it == 0 )
{
t = tmin;
//ch = t^2*( x^3 -1.5x^2+1.0 ) * cos ( (pi/4) * y );
ch=0;
totalmid=0;
//qh = 1+t( x - x^2 ) * sin ( pi * y );//probably don’t need since I directly used qh=1+ch
}
else
{
told = t;
t = ( ( nt - it ) * tmin + it * tmax ) / nt;
tm = ( told + t ) / 2.0;
coldh = ch;
//qoldh = qh;// probably don’t need since I directly used qoldh=1+coldh
//or we can use this t=t+dt;
//t=t+dt;
halfbackwardeulerstep;
ch = (2.0 * cmh - coldh);
totalmid=int2d(Th)(ch);
//if(abs(t-0.1)<1e-10||abs(t-1)<1e-10||abs(t-2)<1e-10||abs(t-3)<1e-10||abs(t-4)<1e-10||abs(t-5)<1e-10||abs(t-6)<1e-10 || abs(t-7)<1e-10 || abs(t-8)<1e-10 || abs(t-9)<1e-10 || abs(t-10)<1e-10 ){
fout << t << " “<< totalmid<< endl;
plot(ch,value=true,wait=0,fill=1,ps=“mid_nonsteady_”+t+”.eps");
}
//
// Compute error.
//
a = sqrt
(
int2d ( Th )
(
( ch - (t^2*( x^3 -1.5*x^2+1.0 ) * cos ( (pi/4) * y )))^2
)
);
b = sqrt (
int2d ( Th ) ( ( dx(ch) - (t^2*(3.0*x^2-3.0*x)*cos((pi/4)*y) ))^2 )
+ int2d ( Th ) ( ( dy(ch) + (t^2*(pi/4)*(x^3 -1.5*x^2+1.0)*sin((pi/4)*y) ))^2 )
);
//is the following error refer to 0,1 norm of (c- ch) of the table in thesis?
c = sqrt ( a^2 + b^2 );
if ( it == 0 || it == nt )
{
cerror00 = cerror00 + 0.5 * a^2;
gradcerror00 = gradcerror00 + 0.5 * b^2;
cerror01 = cerror01 + 0.5 * c^2;
}
else
{
cerror00 = cerror00 + a^2;
gradcerror00 = gradcerror00 + b^2;
cerror01 = cerror01 + c^2;
}
//cerroroo0 =(tmax-tmin)*max ( cerroroo0, a );
cerroroo0 =max ( cerroroo0, a );
}
//
// Ai = sqrt ( L2 integral over X ( E^2 ) )
//
// ||E|| 0,0 = sqrt ( L2 integral over T ( L2 integral over X ( E^2 ) )
// Use trapezoidal rule:
// ||E|| 0,0 = dt * sqrt ( 1/2 A0^2 + A1^2 + A2^2 + … + An-1^2 + 1/2 An^2 )
//
// ||E|| oo,0 = max over T ( sqrt ( L2 integral over X ( E^2 ) )
// = max ( A0, A1, …, An )
//
// cerroroo0=(tmax-tmin)*cerroroo0;
cerror00 = sqrt(dt) * sqrt ( cerror00 );
gradcerror00 = sqrt(dt) * sqrt ( gradcerror00 );
cerror01 = sqrt(dt) * sqrt ( cerror01 );
//c = sqrt ( cerror00^2 + cerror01^2 );
cout << endl;
cout << "||cerror||_oo,0 = " << cerroroo0 << endl;
cout << "||cerror||_0,0 = " << cerror00 << endl;
cout << "||grad_cerror||_0,0 = " << gradcerror00 << endl;
cout << "||cerror||_0,1 = " << cerror01 << endl;
}
//
// Terminate.
//
cout << “\n”;
cout << “forward_euler:\n”;
cout << " Normal end of execution.\n";
exit ( 0 );
