Solution to the Poisson Boltzmann equation

General Discussion
1

Hello community, I am new to using FreeFem++. I have been trying to solve the Poisson-Boltzmann equation, however the code I am using does not solve said equation. I attach the image of the equations I am trying to solve and the code I am using, I would appreciate your suggestions.

EQ

load "medit";
real a=1.5,b=0.75;
int n=300,m=150;
mesh Th=square(n,m,[a*x,b*y]);// Uniform mesh on [0,a]x[0,b]
plot(Th,wait=1,ps="PRmalla0.eps");
fespace Vh(Th,P2);
Vh u=0,v=0;
func f=0;
func Ka=10;
int i=0;
real error=0.00001, coef=0.001^(1./5.);
problem Problem1(u,v,solver=sparsesolver,init=i,eps=1.0e-6)=
int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v)+Ka*Ka*sinh(u)*v)
+on(1,u=1)
+on(2,u=1)+on(3,u=1)+on(4,u=1);
real cpu=clock();
Problem1;
plot(u,wait=2,ps="PRuinicial.eps");

plot(u,wait=1,ps="PRufinal.eps");
plot(Th,wait=1,ps="PRmallafin.eps");
Problem1; // solve the problem plot(uh);
plot(u,ps="Problem1.eps",fill=1);// to see the result
medit("Problem1",Th,u);
plot(u, wait=true, value=true, fill=true, ps="HeatExchanger.eps",dim=3);
2

Your problem is nonlinear. You must use other method like Fixed point, Newton-Raphson or so…

3

I propose this fixed point algorithm:

4

Thanks for the suggestion, I have already implemented the code and it is indeed correct. However, I have tried to export the data in an mxn matrix array, where the columns are the “x” values ​​and the rows are the “y” coordinate.

load "medit";
//Parameters
real Ka=10.;
real Zeta1=1.;
real Zeta2=1.;
real Zeta3=1.;
real Zeta4=1.;
// Mesh
border L1(t=0,1.5){x=t;y=0;}
border L2(t=0,0.75){x=1.5;y=t;}
border L3(t=0,1.5){x=1.5-t;y=0.75;}
border L4(t=0,0.75){x=0;y=0.75-t;}
int n=100;
mesh Th=buildmesh(L1(2*n)+L2(n)+L3(2*n)+L4(n));
//plot(Th,wait=1,ps="rectangulo.eps");

fespace Vh(Th,P2);
Vh u=0, v=0;

real error=0.00001, coef=0.0001^(1./5.);

// Initialize up (previous value)
Vh up;
up = u;

problem Problem1(u, v, solver=sparsesolver) =
    int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v) + Ka*Ka*(u)*v )
  + int2d(Th)( Ka*Ka*(sinh(up)-up)*v )
  + on(L1, u=Zeta1)
  + on(L2, u=Zeta2)
  + on(L3, u=Zeta3)
  + on(L4, u=Zeta4);

for(int iter = 0; iter < 20; ++iter) {
    up[] = u[];
    Problem1; // Solve the problem
    up[] -= u[];
    cout << iter << " error " << up[].linfty << endl;
    if (up[].linfty < 1e-5) break; // Convergence check
}

//plot(u, wait=1, ps="PRufinal.eps");
//plot(Th, wait=1, ps="PRmallafin.eps");
Problem1; // Solve the problem
plot(u, ps="Problem1.eps", fill=1); // To see the result
medit("Problem1", Th, u);
plot(u, wait=true, value=true, fill=true, ps="HeatExchanger.eps", dim=3);