# Solution to the Poisson Boltzmann equation
**URL:** https://community.freefem.org/t/solution-to-the-poisson-boltzmann-equation/3337
**Category:** General Discussion
**Created:** [June 21, 2024, 8:52pm UTC](https://community.freefem.org/t/solution-to-the-poisson-boltzmann-equation/3337 "2024-06-21T20:52:32Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@EJimenez](https://community.freefem.org/u/EJimenez)
#### Post date: [June 21, 2024, 8:52pm UTC](https://community.freefem.org/t/solution-to-the-poisson-boltzmann-equation/3337/1 "2024-06-21T20:52:32Z")
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.
```auto
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);
```
---
### Author:  [@julienG](https://community.freefem.org/u/julienG)
#### Post date: [June 22, 2024, 9:00am UTC](https://community.freefem.org/t/solution-to-the-poisson-boltzmann-equation/3337/2 "2024-06-22T09:00:12Z")
Your problem is nonlinear. You must use other method like Fixed point, Newton-Raphson or so…
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [June 24, 2024, 12:02pm UTC](https://community.freefem.org/t/solution-to-the-poisson-boltzmann-equation/3337/3 "2024-06-24T12:02:34Z")
I propose this fixed point algorithm:
> [@EJimenez](#):
>
> ```auto
> real error=0.00001, coef=0.001^(1./5.);
> Vh up; // previous value
> up=u;
> problem Problem1(u,v,solver=sparsesolver)=
> int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v)+Ka*Ka*(u)*v)
> + int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v)+Ka*Ka* (sinh(up)-up)*v)
> +on(1,u=1)
> +on(2,u=1)+on(3,u=1)+on(4,u=1);
>
> for(int iter=0; iter< 20;++iter)
> {
> up[]=u[];
> Problem1
> up[]-=u[];
> cout << iter << " errr " << up[].linfty < if( up[].linfty < 1e-5) break; // converge
> }
>
> ```
---
### Author:  [@EJimenez](https://community.freefem.org/u/EJimenez)
#### Post date: [June 28, 2024, 9:35pm UTC](https://community.freefem.org/t/solution-to-the-poisson-boltzmann-equation/3337/4 "2024-06-28T21:35:33Z")
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.
```auto
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);
```