How to solve Possion equation have a parameter h for different values of h using loop

General Discussion
1

Dear, I want to solve the Poisson equation involving h for different values of h given in Pram (in the following code). Then, I used a loop to pick each value from Pram as a value for h and then solve the problem. I want to compute AvgNu at different values of h. Kindly see my code and guide me to achieve the required results.

real Mesh=60;
real ra=0.1;
border C(t=0,1){x=0.5+racos(2pit);y=0.2sqrt(3)+rasin(2pit);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real Pram=[1,2,3];
for(int i;i<=3;i++)
{
real h=Pram[i];
fespace Vh(Th,P2);
func f= xy;
Vh U,V;
problem Poisson(U,V,solver=LU)=
int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))

  • hint2d(Th)( fV )

on(C,U=1) ;
Poisson;
real AvgNu;
macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'*[N.x,N.y]);
cout<<"Avg = "<<AvgNu[i]<<endl;
}

1 Like
2

They is no problem

use the following writing:

problem Poisson(U,V,solver=LU)=
int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))
- int2d(Th)( h* f * V )
on(C,U=1) ;

in Poisson problem

3

Dear Professor @frederichecht thanks for your response. I think you don’t understand my point. I want to solve the Poisson equation for different values of h using the loop. For example,
for i=1, the code should solve the Poisson problem for h = 1 and compute the value for AvgNu,
for i=2, the code should solve the Poisson problem for h = 2 and compute the value for AvgNu, and
for i=3, the code should solve the Poisson problem for h = 3 and compute the value for AvgNu.

4

No , for( i = 0; i<3; ++i) read the doc please.

5

That means we can’t solve using the loop.

6

the corrected code with only factorisation (init=0) because the matrix is constante here.

// M_Usman Muhammad Usman
real Mesh=60;
real ra=0.1;
border C(t=0,1){x=0.5+ra*cos(2*pi*t);y=0.2*sqrt(3)+ra*sin(2*pi*t);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real[int]  Pram=[1,2,3],AvgNu(3);
real h;
fespace Vh(Th,P2);
func f= x*y;
Vh U,V;
int i=0;
problem Poisson(U,V,solver=LU,init=i)=
    int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))
    - int2d(Th)( h*f*V )
    + on(C,U=1) ;
    verbosity=0; 
for( i=0;i<3;i++)
{
  h=Pram[i];

  Poisson;
  plot(U,wait=1);
  macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'*[N.x,N.y]);
cout<<i << "Avg = "<<AvgNu[i]<<endl;
}

the output :

ff++-13/bin/FreeFem++ 'Usman.edp' -glut '/usr/local/ff++-13/bin/ffglut' ;exit;
-- FreeFem++ v4.13 (Mer 20 sep 2023 11:48:02 CEST - git v4.13-141-g86b66a0e)
   file : Usman.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : // M_Usman Muhammad Usman
    2 : real Mesh=60;
    3 : real ra=0.1;
    4 : border C(t=0,1){x=0.5+ra*cos(2*pi*t);y=0.2*sqrt(3)+ra*sin(2*pi*t);}
    5 : mesh Th=buildmesh(C(Mesh));
    6 : plot(Th,wait=1);
    7 : real[int]  Pram=[1,2,3],AvgNu(3);
    8 : real h;
    9 : fespace Vh(Th,P2);
   10 : func f= x*y;
   11 : Vh U,V;
   12 : int i=0;
   13 : problem Poisson(U,V,solver=LU,init=i)=
   14 :     int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))
   15 :     - int2d(Th)( h*f*V )
   16 :     + on(C,U=1) ;
   17 :     verbosity=0; 
   18 : for( i=0;i<3;i++)
   19 : {
   20 :   h=Pram[i];
   21 : 
   22 :   Poisson;
   23 :   plot(U,wait=1);
   24 :   macro grad(u) [dx(u),dy(u)] )  //
   25 : AvgNu[i]=-int1d(Th,C)((grad(U)  [dx(U),dy(U)])'*[N.x,N.y]);
   26 : cout<<i << "Avg = "<<AvgNu[i]<<endl;
   27 : } sizestack + 1024 =2168  ( 1144 )

  --  mesh:  Nb of Triangles =    638, Nb of Vertices 350
0Avg = 0.00545748
1Avg = 0.010915
2Avg = 0.0163724
7

Thank you, dear Professor @frederichecht . For nonlinear PDEs, I got the same value for AvgNu. Maybe I did something wrong in the following code. Kindly have a look and guide me to overcome this error.

real Mesh=60;
real ra=0.1;
border C(t=0,1){x=0.5+racos(2pit);y=0.2sqrt(3)+rasin(2pit);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real a=1;
real err=1;
real[int] Pram=[1,2,3],AvgNu(3);
real b;
fespace Vh(Th,P2);
func f= xy;
Vh U,U0,V;
U0=0;
int i=0;
problem Poisson(U,V,solver=LU,init=i)=
int2d(Th)(dx(U)dx(V) + dy(U)dy(V))
-int2d(Th)(abfU0UV)
+on(C,U=1) ;
verbosity=0;
for( i=0;i<3;i++)
{
b=Pram[i];
for(int ij=0; (ij)<=10; ij++)
while (err>1e-4)
{
Poisson;
err=int2d(Th)(U-U0)^2;
U0 = U;
endl;
}
macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'[N.x,N.y]);
cout<<i << "Avg = "<<AvgNu[i]<<endl;
}

8

For non linear equation first you have to choose which algorithm you use to solve your problem
fixed point , Newton, …

so you put a loop in a loop…

9

I solved the nonlinear problem as:
for(int ij=0; (ij)<=10; ij++)
while (err>1e-4)
{
Poisson;
err=int2d(Th)(U-U0)^2;
U0 = U;
endl;
}
and then put it in the loop to compute the solution for various values of b. The got the solution successfully but the solution is same for each value of b.

10

Normal, your problem do not change in the iteration process, and the err must set before big value.

11

I already set err=1 here, then I used it in the iterative process.

12

Dear Professor @frederichecht , I also defined err=1 just after the value b and before starting the second loop. But, I got the same values for AvgNu. Where I have to define the value err so that I get the correct results for AvgNu. I am waiting for your positive response.

13

Dear Professor @frederichecht , what do you meant by big value.

14

bigger than 1e-4 to enter , in the loop. so 1 for exemple…

15

Dear Professor @frederichecht , I followed your instructions (code given below) but found the same results for AvgNu.

real Mesh=100;
real ra=0.1;
border C(t=0,1){x=0.5+racos(2pit);y=0.2sqrt(3)+rasin(2pit);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real a=1;
real[int] Pram=[1,20,30],AvgNu(3),E=[1,1,1];
real b;
fespace Vh(Th,P2);
func f= xy;
Vh U,U0,V;
U0=0;
int i=0;
problem Poisson(U,V,solver=LU,init=i)=
int2d(Th)(dx(U)dx(V) + dy(U)dy(V))
-int2d(Th)(abfU0UV)
+on(C,U=1) ;
verbosity=0;
for( i=0;i<3;i++)
{
b=Pram[i];
real err=1;
for(int ij=0; (ij)<=10; ij++)
while (err>1e-7)
{
cout<<err<<endl;
Poisson;
err=int2d(Th)(U-U0)^2;
U0 = U;
endl;
}
macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'[N.x,N.y]);
cout<<i << "Avg = "<<AvgNu[i]<<endl;
}

16

Dear Professor @frederichecht , kindly see this and guide me to fix the query.

17

I have no idea of your No linear problem, because you only put linear problem.
So it is difficult to help you.