Poisson-Boltzmann Equation

[nassboy]

I've been working on finding a closed form solution or an approximate solution for a variant of the Poisson-Boltzmann equation.

y''[x]=A*(1-e^(-y[x])

I'm unable to solve it using methods such as separation of variables, but I have been able to generate a power series using recursive derivatives. I've pasted the mathematica code below...I've been unable to find the pattern to generate the series...can anybody help?

In[2]:= (*Define Recursive Derivative Function *)

In[3]:= Derivative[N_][Y][0]:=D[Y''[x],{x,N-2}] /. x->0

In[4]:= (*Define the Problem*)

In[5]:= Y''[x]=A*(1-E^(-Y[x]))

Out[5]= A (1-E^-Y[x])

In[6]:= Y''[0]=Y''[x] /. x->0

Out[6]= A (1-E^-Y[0])

In[7]:= Y'[0]=0

Out[7]= 0

In[8]:= (*Check the Series*)

In[9]:= Y[0]+Integrate[Integrate[Series[Y''[x],{x,0,10}],x],x]

Out[9]= Y[0]+1/2 (A-A E^-Y[0]) x^2+1/24 A^2 E^(-2 Y[0]) (-1+E^Y[0]) x^4-1/720 (A^3 E^(-3 Y[0]) (4-7 E^Y[0]+3 E^(2 Y[0]))) x^6+(A^4 E^(-4 Y[0]) (-34+82 E^Y[0]-63 E^(2 Y[0])+15 E^(3 Y[0])) x^8)/40320-((A^5 E^(-5 Y[0]) (496-1510 E^Y[0]+1638 E^(2 Y[0])-729 E^(3 Y[0])+105 E^(4 Y[0]))) x^10)/3628800+(A^6 E^(-6 Y[0]) (-11056+40540 E^Y[0]-56568 E^(2 Y[0])+36684 E^(3 Y[0])-10545 E^(4 Y[0])+945 E^(5 Y[0])) x^12)/479001600+O[x]^13

 

[coomast]

I've been working on finding a closed form solution or an approximate solution for a variant of the Poisson-Boltzmann equation.

y''[x]=A*(1-e^(-y[x])

...

Hello nassboy,

I am not able to help you in finding a solution using mathematica but I can help you in finding a solution analytically. That is to say, the final integral you need to solve numerically or by another method. OK, so first substitute the following variables:
$$u=y \qquad v=y' \qquad w=\frac{dv}{du}=\frac{y''}{y'}$$
or:
$$y=u \qquad y'=v \qquad y''=v\frac{dv}{du}$$

This comes from the fact that the original equation does not have x explicitely used. This gives now:

$$v\frac{dv}{du}=A*\left(1-e^{-u}\right)$$

And solved for v:

$$v=\sqrt{K_1+2A*\left(u+e^{-u}\right)}$$

Or using the opposite of the substitution:

$$\frac{dy}{dx}=\sqrt{K_1+2A*\left(y+e^{-y}\right)}$$

From which:

$$x+K_2=\int \frac{dy}{\sqrt{K_1+2A*\left(y+e^{-y}\right)}}$$

In which $K_1$ and $K_2$ are the integration constants. This final integral is probably not solvable in an ordinary way, so other methods are necessary, perhaps it is the starting point for the series you are looking for.

Hope this helps,

coomast

 

[nassboy]

Thanks for looking at my problem...

The reason I went to a series approach in the first place was that I didn't know how to tackle that ugly integral. It has very steep curvature so quadrature methods don't work well...

 

[coomast]

Thanks for looking at my problem...

The reason I went to a series approach in the first place was that I didn't know how to tackle that ugly integral. It has very steep curvature so quadrature methods don't work well...

Mmm, if really necessary I can look up the quadrature methods for finding integrals in order to come up with a formula for finding it numerically with the greatest degree of algebraic accuracy -if I recall well- it was called. This can take a day or two because this has been some time. Maybe I can find a different -and for this case- a better way of solving it. Do you want me to look for something?

coomast

 

[nassboy]

Nah,

I'm more interested in a closed form approximate solution...

some function Y[x]=, that is close to the real solution

 

[coomast]

Nah,

I'm more interested in a closed form approximate solution...

some function Y[x]=, that is close to the real solution

I will have a look at it tomorrow. I'm not promising anything but I will give it a try.

coomast

 

[nassboy]

thanks!

 

[coomast]

hello nassboy,

I had a look at the integral again and came up with a closed form solution. However this is not an easy solution and difficult to work with, although, if you are going to use it in a program it might be the solution you are looking for. But first some relations I used in deriving this solution.
$$\frac{1}{\sqrt{1+x}}=1-\frac{1}{2}x+\frac{1 \cdot 3}{2 \cdot 4}x^2- \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6}x^3+...= \sum_{i=0}^{\infty}\frac{(2i)!\cdot (-1)^i}{\left(2^i\cdot i!\right)^2}x^i$$

$$(A+B)^i=\sum_{k=0}^{i}\frac{i!}{k!\cdot (i-k)!}A^k\cdot B^{i-k}$$

$$e^{-by}=\sum_{p=0}^{\infty}\frac{(-by)^p}{p!}$$

So the first step is to rewrite the integral as:

$$x+K_2=\frac{1}{\sqrt{K_1}}\int \frac{dy} {\sqrt{1+\frac{2A}{K_1}(y+e^{-y})} }$$

Expanding the denominator using the first relation gives:

$$x+K_2=\frac{1}{\sqrt{K_1}} \sum_{i=0}^{\infty} \frac{(2i)!\cdot (-1)^i}{\left(2^i\cdot i!\right)^2} \cdot \left(\frac{2A}{K_1}\right)^i \cdot \int (y+e^{-y})^i dy$$

This under the assumption that the switching of the summation and the integral is allowed, meaning that it must converge. I didn't check this. Using the second relation we can rewrite this new integral as:

$$\int (y+e^{-y})^i dy=\sum_{k=0}^{i} \frac{i!}{k!\cdot (i-k)!} \int {y^ke^{-(i-k)y}}dy$$

The final integral now can be found by the third relation:

$$\int {y^ke^{-(i-k)y}}dy=\sum_{p=0}^{\infty} \frac{(-1)^p (i-k)^p}{p!} \int y^{k+p}dy= \sum_{p=0}^{\infty} \frac{(-1)^p (i-k)^p}{p!}\cdot \frac{y^{k+p+1}}{k+p+1}$$

The complete solution now becomes:

$$x+K_2=\frac{1}{\sqrt{K_1}} \sum_{i=0}^{\infty} \left[ \frac{(2i)!\cdot (-1)^i}{\left(2^i\cdot i!\right)^2} \cdot \left(\frac{2A}{K_1}\right)^i \cdot \right. \\ \left. \left(\sum_{k=0}^{i} \left[\frac{i!}{k!\cdot (i-k)!} \cdot \left( \sum_{p=0}^{\infty} \frac{(-1)^p (i-k)^p}{p!}\cdot \frac{y^{k+p+1}}{k+p+1} \right)\right]\right)\right]$$

I can't garantee that this is without error. As stated before the convergence checks needed for getting approval for the switching between the integral and the sums are not done.

Hope this helps a bit,

coomast

 

[nassboy]

Wow I'm impressed! I wouldn't have come up with that in ages! I'll try to verify it tomorrow or Monday!

I have one more, perhaps, naive question. Is one able to invert a series like that to get a function for y, for example y(x)=Series?

 

[nassboy]

I tried to implement the series solution in Mathematica and did not get the results that I expected. Perhaps the series does not converge...perhaps I made a mistake.

 

[Mute]

Maybe you can use the http://en.wikipedia.org/wiki/Lagrange_inversion_theorem" to get a series expansion?

If you set your initial conditions to be $y(x_0) = y_0$ and $y'(x_0) = y'_0 \neq 0$, then you can solve for x in terms of the integral over y as coomast did (I kept the plus/minus sign from taking the square root):

$$x = x_0 \pm \int_{y_0}^{y}\frac{du}{\sqrt{2A(u + e^{-u}) + C_0}} \equiv f(y)$$

where $C_0 = -2A(y_0 + e^{-y_0}) + (1/2)(y'_0)^2$. Then we want y = g(x), which we can find with the inversion theorem. Now, hopefully the solution is analytic at the initial value (or anywhere - if it's not then the whole thing breaks down), so referring to the wikipedia article, let's choose $a = y(x_0) = y_0$ and $b = x_0$. So, then by the theorem,

$$g(x) = y_0 + \sum_{n=1}^{\infty}\frac{d^{n-1}}{dy^{n-1}}\left.\left(\frac{y-y_o}{f(y)-x_0}\right)^n\right|_{y=y_0}\frac{(x-x_0)^n}{n!}$$.

In the first order (n=1) term in the sum we have to use L'Hopital's rule as y -> y_0. Doing so, I get,

$$y = g(x) = y_0 \pm\sqrt{2A(y_0+e^{-y_0})+C_0}(x-x_0) + \sum_{n=2}^{\infty}\frac{d^{n-1}}{dy^{n-1}}\left.\left(\frac{y-y_o}{f(y)-x_0}\right)^n\right|_{y=y_0}\frac{(x-x_0)^n}{n!}$$.

The higher order terms look like they'll be ugly to evaluate, but in principle you can churn out each term in the expansion to get your series expansion. You would likely want to check your result against numerical solutions to the equation to make sure it's capturing the correct behaviour at small values of y, as using the theorem does assume that f(y) is analytic when I haven't actually checked that.