Poisson Boltzmann eq. with an infinite boundary condition

Question

Roi Bar-On 2020-10-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/623238-poisson-boltzmann-eq-with-an-infinite-boundary-condition

评论： Roi Bar-On 2020-11-3

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

Hi everybody!

I've written a code that solves the following PB problem

This is a 2nd order ode and the dependant variable here is psi (tilde) and the independant variable is x (tilde). The boundary conditions are:

The code includes an optimizer based on the secant method. I'm trying to guess the solution that will ''land'' on the 2nd boundary condition via the shooting method:

function [x2,p2]=optimumup(y21,y22,bc2,l)
%optimizing the choice of y2 for the shooting method, using the Newton
%method (secant) by solving the equation p(length(p(:,1))-bc2=0
[x1,p1] = BP(y21,l);
[x2,p2] = BP(y22,l)
size(p2)
for i=1:10
dp=((p2(length(p2(:,1)),1)-bc2)-(p1(length(p1(:,1)),1)-bc2))/(y22-y21);
y23=y22-(p2(length(p2(:,1)),1)-bc2)/dp;%recurrence relation
p1=p2;
[x2,p2] = BP(y23,l);
if abs(p2(length(p2(:,1)))-bc2)<0.001;
    disp(p2(length(p2(:,1)),1));
    disp(p2(length(p2(:,1)),1)-bc2);
    plot(x2,p2(:,1),x2,p2(:,2))
    return 
end
    y21=p1(1,2);
    y22=p2(1,2);
end
disp('No success')
end
function [x,p] = BP(y2,l)
%Shooting method for the BP eq
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[x,p] = ode45(@eqs,[0 l],[90/24 y2],options);
%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2))
end
function [dy] = eqs(x,y)
%PB equation
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = sinh(y(1));
end

y21 and y22 are the initial x guesses for the secant method. bc2 is the boundary condition at x goes to inifinity. l is the interval size.

I have two main problems:

I don't know how to incorperate the inifinite boundary condition.
For some reason my solution does not converge. I'm working dimensionless.

I would appreciate some help guys.

Thanks,

Roi

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2020-10-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/623238-poisson-boltzmann-eq-with-an-infinite-boundary-condition#answer_522663

在 MATLAB Online 中打开

I suspect the initial gradient required to get zero at infinity is an irrational number (or possibly a rational number with a large number of digits in its decimal expansion!). See the following for example:

xf = 25;
xspan = [0 xf];
dydx0 = -1.042119689394;
dx = 10^-11;
dydx0hi = dydx0+dx;
dydx0lo = dydx0-dx;
Y0 = [1, dydx0];
[x,Y] = ode45(@eqs,xspan,Y0);
psi = Y(:,1);
Y0hi = [1, dydx0hi];
[xhi,Yhi] = ode45(@eqshi,xspan,Y0hi);
psihi = Yhi(:,1);
Y0lo = [1, dydx0lo];
[xlo,Ylo] = ode45(@eqslo,xspan,Y0lo);
psilo = Ylo(:,1);
figure(1)
plot(x,psi,'k',xhi,psihi,'r',xlo,psilo,'b'),grid
axis([0 xf -2 2])
xlabel('x'),ylabel('\psi')
lbl1 = ['d\psi dx(x=0) = ',num2str(dydx0,12)];
lbl2 = ['d\psi dx(x=0) = ',num2str(dydx0hi,12)];
lbl3 = ['d\psi dx(x=0) = ',num2str(dydx0lo,12)];
legend(lbl1,lbl2,lbl3)
function dY = eqs(~,Y)
%PB equation
dY = [Y(2);
      sinh(Y(1))];
end
function dY = eqshi(~,Y)
%PB equation
dY = [Y(2);
      sinh(Y(1))];
end
function dY = eqslo(~,Y)
%PB equation
dY = [Y(2);
      sinh(Y(1))];
end

This produces the following graph:

Although it might look as though the black line has converged on zero, if you extend the end time to, say 30, you will see it start to diverge.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Alan Stevens 2020-10-24

In effect I used the shooting method, except that I did it manually rather than writing a routine to do it! It’s actually not difficult, and took much less time than I would have taken to write a routine. Had the result looked as though the curve was tending to zero in a stable manner, I would then have written a routine to do it more accurately. However, as you can see from the fact that small perturbations in the initial slope lead to divergences in opposite directions, it seems that the final equilibrium is likely to be unstable.

A boundary condition at infinity suggests an asymptotic approach to zero (in this case). In physical systems, such as fluid flow boundary layers, for example, it is a convenient fiction standing in for a finite distance, so can often be replaced by a (relatively) large but finite value. If that is the case for you then you might be able to make use of the result. However, your equation doesn’t look as though it represents a physical system (I could be wrong!), so you will have to decide if it is useful or not.

Roi Bar-On 2020-10-24

Physically it describes the potential between two surfaces where one has a dimensionless potential of 1 and the other one has the potential of zero. As a matter a fact I'm solving almost similar problem with the same boundary conditions so i need a robust algorithm to work here.. I am seeking to incorporate my infinite boundary condition within the initial guess in the code as I mentioned.

Alan Stevens 2020-10-24

So why do you need to set an infinite boundary condition, rather than use the distance between the surfaces?

请先登录，再进行评论。

Answer 2

David Goodmanson 2020-10-24

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/623238-poisson-boltzmann-eq-with-an-infinite-boundary-condition#answer_523308

编辑：David Goodmanson 2020-10-24

在 MATLAB Online 中打开

Hi Roi

For convenience I will use u in place of psi. An exact solution is

tanh(u/4) = A*exp(-x)    where  A = tanh(1/4)             [1]
so 
u = 4*atanh(A*exp(-x))

which meets both of the boundary conditons. For du/dx at the origin, it turns out that

du/dx = -2*sinh(1/2) = -1.042190610987495

which disagrees slightly from what Alan has, for reasons not known to me.

===========derivation=============

start with

u'' = sinh(u)

multiply by integrating factor u' and integrate

u'^2/2 = cosh(u)      (+ constant)
u'^2 = 2*cosh(u) + constant

as x --> inf you want u --> 0 which also means that u' --> 0, so the constant = -2

u'^2 = 2*(cosh(u)-1) = 4*sinh(u/2)^2
u' = -2*sinh(u/2)      % pick the negatatve root
du/(2*sinh(u/2)) = -dx
% integrate
log(tanh(u/4)) = -x + constant
tanh(u/4) = A*exp(-x)

you want u = 1 at x = 0 so

A = tanh(1/4)

For the derivative at the origin, take the derivative of both sides of [1] at the top

( (1/4) / cosh^2(u/4) ) du/dx = -A*exp(-x)
% evaluate this at x = 0 and u(0) = 1
( (1/4) / cosh^2(1/4) ) du/dx = -A
du/dx = -tanh(1/4)*4*cosh^2(1/4) = -4*sinh(1/4)*cosh(1/4) = -2*sinh(1/2)

2 个评论
显示无隐藏无

Alan Stevens 2020-10-24

Nicely done David!

David Goodmanson 2020-10-24

Hi Alan,

Thanks, it did work in this special case but I think your comments about a finite boundary and/or an asymptotic approach are going to be more useful in general. In this case out at infinity the answer has to be u = A*exp(-x), so both u and u' are known there. Then starting at large x, integrating to the origin and using the shooting method for various A should also work, I think. I have not tried it.

请先登录，再进行评论。

Answer 3

Roi Bar-On 2020-10-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/623238-poisson-boltzmann-eq-with-an-infinite-boundary-condition#answer_523273

Because my problem is for a half space by definition..

I don't mind using a finite value for ''infinity'', it's okay. In the case you've presented x=5 is infinity and that is okay. the thing is that in my code as I attached, even if I'll take x=5 as ''infinity'', the solution won't converge. Have you tried to run it and got something maybe?

2 个评论
显示无隐藏无

Alan Stevens 2020-10-24

No, but I'll have a look later today.

Roi Bar-On 2020-10-24

Thank you so much!

请先登录，再进行评论。

Answer 4

Roi Bar-On 2020-10-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/623238-poisson-boltzmann-eq-with-an-infinite-boundary-condition#answer_524038

Thank you very much guys!

Analytical solution is great!

As for the numerical one with ode45 that I have attached. How can I guess the two x's that secant method requires correctly in order to converge to the same analytical solution?

Roi

8 个评论
显示 6更早的评论隐藏 6更早的评论

Alan Stevens 2020-10-25

在 MATLAB Online 中打开

The single point secant method only requires one guess together with an initial increment. I've modified your code to use this (its "cleaner" than the two point secant method), and compared the result with the analytical solution:

tolerance = 10^-8;  % Set as desired
xspan = [0 5];  % inf = 5 here!
% Initial guess for gradients  dpsidx(x=0)
dp1 = -1;
delta = -10^-10;  
flag = 1; % Repeat indicator
itmax = 100; % Maximum allowable number of iterations
its = 0;
while (flag==1) && its<itmax 
    % Call ode solver twice with different guesses for initial gradient
    % and get values of psi at boundary in return psi(x=inf)
    [~, p1] = BP([1 dp1],xspan);
    P1 = p1(end);
    [x, p2] = BP([1 (dp1+delta)],xspan);
    P2 = p2(end);
       
    if abs(P2)<tolerance  % then we are finished
          flag = 0;
    else  % we need another iteration
          % Secant improvement
          dP = P2 - P1;  % Difference between end values of psi
          dp1 = dp1 - delta*P1/dP; 
    end
    
    its = its + 1;
end
% Display data
format long
disp('psi(x=inf) and dpsidx(x=0)')
disp([P2 dp1])
disp(['number of iterations = ' num2str(its)])
% Analytical solution
 A = tanh(1/4);
 psi = 4*atanh(A*exp(-x));
 % Plot numerical and analytical comparison
 plot(x,p2(:,1),x,psi),grid
 xlabel('x'),ylabel('\psi')
 legend('numerical','analytical')
function [x,p] = BP(y0,xspan)
%Shooting method for the BP eq
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[x,p] = ode45(@eqs, xspan, y0 ,options);
%n=size(p);p(n(1),1)
%plot(x,p(:,1),x,p(:,2)),grid
%legend('\psi','d\psi dx')
end
function dy = eqs(~,y)  % x not used within the function so replaced by ~
%PB equation
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = sinh(y(1));
end

This gives the following graph:

(Interestingly, if you just plot exp(-x), you get a curve extremely close to the analytical solution also!)

Alan Stevens 2020-10-30

Yes, one of your boundary conditions is that psi(x=infinity) = 0. So we try a value of dpsi/dx(x=0) and see what value of psi(x = infinity) we get. If the latter isn't zero, we try another value of dpsi/dx(x=0) and keep repeating this process until we hit on a value that makes psi(x=infintiy) zero. Instead of trying random guesses for dpsi/dx(x=0), the secant method provides a way of improving subsequent guesses (though the method isn't infallible!).

Roi Bar-On 2020-11-3

Thank you!

请先登录，再进行评论。

Poisson Boltzmann eq. with an infinite boundary condition

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（3 个）

2 个评论
显示无隐藏无

2 个评论
显示无隐藏无

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Poisson Boltzmann eq. with an infinite boundary condition

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（3 个）

2 个评论 显示 无隐藏 无

2 个评论 显示 无隐藏 无

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

2 个评论
显示无隐藏无

2 个评论
显示无隐藏无

8 个评论
显示 6更早的评论隐藏 6更早的评论