Why is bvp4c not able to faithfully solve the Poisson Boltzmann equation ?

3 views (last 30 days)
Hi ,
I am trying to numerically solve the Poisson-Boltzmann equation for electrolytes in a solution using MATLAB. The equation can be written as 2nd order non linear ordinary differential equation.
Some of the relevant terms and symbols (as used in my code are explained below.
I have got 2 boundary conditions : psi = 0 at x = 0 and psi = 0 at x = inf ( or a very large distance away from x = 0)
Below is my MATLAB code :-
function [sol] = pbeexp()
el = 4.8*10^-10; %electron charge in StatCoulomb in CGS
eps = 80; %permittivity of water at 20 degrees C
k = 1.3807*10^-16; %Boltzmann constant in CGS
T = 293 ; %absolute temperature at 20 degrees C
l = el^2/(eps*k*T);%Bjerrum length in cm
pn = k*T/el; %non dimensional potential
psiwall =300*10^-3*0.00333564; %Electric field at the electrode in StatVolts
%Formulating the problem
x=linspace(0,35*10^-7,10000)/l;%Linearly distributed grid
solinit = bvpinit(x,[psiwall/pn 0]); % Initial mesh
options = bvpset('RelTol',1e-5,'AbsTol',1e-5); % Tolerance values
sol = bvp4c(@pbeode,@pbebc,solinit,options);%Initial solution structure
xint=linspace(sol.x(1),sol.x(end),size(sol.x,2));%Linearly distributed grid
sxint=deval(sol,xint);
dist = xint;
pot = sxint(1,:);
plot(dist,pot)
function dydx = pbeode(x,psi)
el = 4.8*10^-10; %electron charge in StatCoulomb in CGS
eps = 80; %permittivity of water at 20 degrees C
k = 1.3807*10^-16; %Boltzmann constant in CGS
T = 293 ; %absolute temperature at 20 degrees C
l = el^2/(eps*k*T);%Bjerrum length in cm
H = 7.8*10^-7; %length of probe in cm
a3 = (1/4.6)*10^3 ; %steric packing in cc
czero = 6.023*10^20*0.1 ; %bulk counterion concentration
nu = 2*czero*a3; %steric size parameter for use in solution
%%Non - dimensional quantities
pn = k*T/el; %non dimensional potential
%%PBE terms for different components of the model
A1 = (4*pi*czero*l*sinh(psi(1)))/(1+2*nu*((sinh(psi(1)/2)^2)));% expression for counterions in the solution
dydx = [ psi(2)
(2*l^2*A1)];
end
%%Providing the B.C.s
function bc = pbebc(psia,psib)
psiwall = 300*10^-3*0.00333564; %Electric field at the electrode in StatVolts
pn = k*T/el; %non dimensional potential
bc = [ psia(1)-psiwall/pn
psib(1)];
end
end
When I run the code the resulting plot of "psi" vs "x" is coming as linear.
However according to published data it should decay down to zero in an exponential manner!! One link is provided below :- http://www.kirbyresearch.com/index.cfm/wrap/textbook/microfluidicsnanofluidicsse55.html.
Could anyone please run the code and instruct me in troubleshooting ? What is wrong with my modelling ?

Answers (1)

AG4102M
AG4102M on 31 Jul 2014
Can anyone please help ? If I want to supply analytical Jacobians then how to go about it ?? I have read the documentation but it is not clear to me.
Thanks.

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!