- I put all code into one MATLAB file (named my_pde.m)
- You need to pass all required parameters to relevant functions
- The output vector in pbebc needs to be a 2-by-1, I simply "fixed" it as shown above. Please verify.
Help regarding use of BVP4C in solving an ODE.
3 views (last 30 days)
Show older comments
This is a script I have written :-
clear all
e = 1.602*10^-19; %electron charge
eps = 8.85*10^-12; %permittivity of free space
k = 1.38*10^-23; %boltzmann constant
T = 293 ;%absolute temperature
a3 = 4.6 ; %steric packing
czero = 0.1 ; %bulk counterion concentration
nu = 2*czero*a3;%steric size parameter
l = e^2/(4*pi*eps*k*T);%Bjerrum length
H = 3.8*10^-9;
psival = -10;
newpbe1(l,H,psival)
*Firstly a few constants are defined. Then in the end a function "newpbe1" is called.
The code for the function "newpbe1" is given below. It is saved in a different file. Also the two different functions "pbeode " and "pbebc" are written in two different files and are called in the last line of "newpbe1" where the "bvp4c" function is used.*
function [x,psi] = newpbe1(l,H,psival)
x = [0,logspace(-15,log10(H/2),2e3)]*1/l; % Exponentially distributed grid
solinit = bvpinit(x,[psival 0]); % Initial mesh
options = bvpset('RelTol',1e-7,'AbsTol',1e-7); % Tolerance values
psi = bvp4c(@pbeode,@pbebc,solinit,options);
end
function dydx = pbeode(x,psi,l,czero,nu)
dydx = [ psi(2);(2*l^2*((czero*sinh(psi(1)))/(1+2*nu*(sinh(psi(1)/2))^2)))];
end
function bc = pbebc(psia,psib,psival)
bc = [ psia(1)-psival ; psib(1)];
end
- When I am running this code I am getting the following error message. *
Error using pbeode (line 2)
Not enough input arguments.
Error in bvparguments (line 106)
testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 130)
[n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] = ...
Error in newpbe1 (line 7)
psi = bvp4c(@pbeode,@pbebc,solinit,options);
Error in newpbe (line 14)
newpbe1(l,H,psival)
Kindly help me in resolving the problem. I am unable to fix it. Thanks.
0 Comments
Accepted Answer
Mischa Kim
on 11 Jun 2014
AGNIVO, check out:
function my_pde()
e = 1.602*10^-19; %electron charge
[...] %%remove code for better readability
psival = -10;
newpbe1(l,H,psival,czero,nu) %%need to pass all parameter
end
function [x,psi] = newpbe1(l,H,psival,czero,nu)
x = [0,logspace(-15,log10(H/2),2e3)]*1/l; % Exponentially distributed grid
solinit = bvpinit(x,[psival 0]); % Initial mesh
options = bvpset('RelTol',1e-7,'AbsTol',1e-7); % Tolerance values
par_vec = [l; czero; nu]; %%package parameter in vector
psi = bvp4c(@pbeode,@pbebc,solinit,options,par_vec);
end
function dydx = pbeode(x,psi,par_vec)
l = par_vec(1);
czero = par_vec(2);
nu = par_vec(3);
dydx = [ psi(2);(2*l^2*((czero*sinh(psi(1)))/(1+2*nu*(sinh(psi(1)/2))^2)))];
end
function bc = pbebc(psia,psib,psival)
bc = [psia(1) - psival(1); psib(1)]; %%this needs to be a 2-by-1
end
More Answers (0)
See Also
Categories
Find more on Boundary Value Problems 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!