MATLAB Answers

Poisson Boltzmann eq. with an infinite boundary condition

20 views (last 30 days)
Roi Bar-On
Roi Bar-On on 22 Oct 2020
Commented: Roi Bar-On on 3 Nov 2020
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:
  1. I don't know how to incorperate the inifinite boundary condition.
  2. For some reason my solution does not converge. I'm working dimensionless.
I would appreciate some help guys.
Thanks,
Roi

Accepted Answer

Alan Stevens
Alan Stevens on 23 Oct 2020
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 Comments
Alan Stevens
Alan Stevens on 24 Oct 2020
So why do you need to set an infinite boundary condition, rather than use the distance between the surfaces?

Sign in to comment.

More Answers (3)

David Goodmanson
David Goodmanson on 24 Oct 2020
Edited: David Goodmanson on 24 Oct 2020
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 Comments
David Goodmanson
David Goodmanson on 24 Oct 2020
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.

Sign in to comment.


Roi Bar-On
Roi Bar-On on 24 Oct 2020
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?

Roi Bar-On
Roi Bar-On on 25 Oct 2020
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 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!