# Midpoint method giving flickering result

**URL:** <https://community.freefem.org/t/midpoint-method-giving-flickering-result/1743>\
**Category:** General Discussion\
**Created:** [May 13, 2022, 10:56pm UTC](https://community.freefem.org/t/midpoint-method-giving-flickering-result/1743 "2022-05-13T22:56:23Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![fas41](https://avatars.discourse-cdn.com/v4/letter/f/5e9695/32.png) [@fas41](https://community.freefem.org/u/fas41)\
**Post date:** [May 13, 2022, 10:56pm UTC](https://community.freefem.org/t/midpoint-method-giving-flickering-result/1743/1 "2022-05-13T22:56:23Z")

</div>

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](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](https://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.5_x^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.5_x^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 );

 ![Screen Shot 2022-05-13 at 6.52.18 PM](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/67dd826370f923ca0632a3f53c51508d55f85458.jpeg)

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [May 16, 2022, 11:36am UTC](https://community.freefem.org/t/midpoint-method-giving-flickering-result/1743/2 "2022-05-16T11:36:10Z")

</div>

Yes it is well know that the crank-nicolson method (mid point method) is at the limit for the stability condition, and have no dump on oscillation in time.

Use a other schema like BDF 2 for example to remove occupation in time.

---

<div class="post-metadata">

**Author:** ![fas41](https://avatars.discourse-cdn.com/v4/letter/f/5e9695/32.png) [@fas41](https://community.freefem.org/u/fas41)\
**Post date:** [May 18, 2022, 1:49pm UTC](https://community.freefem.org/t/midpoint-method-giving-flickering-result/1743/3 "2022-05-18T13:49:27Z")

</div>

Actually in my theory , i have worked on backward euler and midpoint, so i wanted to show numerically that midpoint gives good results . But in case of midpoint it goes back and forth