Help regarding use of BVP4C in solving an ODE.

2 次查看(过去 30 天)
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.

采纳的回答

Mischa Kim
Mischa Kim 2014-6-11
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
  • 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.

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Boundary Value Problems 的更多信息

产品

Community Treasure Hunt

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

Start Hunting!

Translated by