- 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.
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.
0 个评论
采纳的回答
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!