MATLAB: Poisson Boltzmann eq. with an infinite boundary condition

ode45optimization

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

Best Answer

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.
Related Question