Solve differential equations with bvp4c
Show older comments
I would like to solve 3 differential equations with bvp4c, but I kept getting error messages as:
Error in bvparguments (line 105)
testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 130)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in BPM_ode_Dec1417 (line 39)
sol=bvp4c(@odefun,@odebc,solinit,options);
Here's what I wrote:
solinit=bvpinit(linspace(0,LA,1000),[E_A U_A C3A]);
options=bvpset('RelTol',1e-3,'AbsTol',[1e-6 1e-6 1e-6]);
sol=bvp4c(@odefun,@odebc,solinit,options);
The function file is:
function [dydx]=odefun(y)
global F R T Z1 Z2 Z3 Z4 C1A C2A i U_A FixedCharge D3 e D4
C1=C1A.*exp(-Z1.*F/R/T.*(y(2)-U_A));
C2=C2A.*exp(-Z2.*F/R/T.*(y(2)-U_A));
dydx(1)=F/e*(C1-C2+10^(-8)/y(3)-y(3)+FixedCharge);
dydx(2)=-y(1);
J4=-D4*(-10^(-8)/y(3)^2*dydx(3)+Z4*F/R/T*10^(-8)/y(3)*dydx(2));
dydx(3)=-((i/F-Z4*J4)/Z3/D3)-Z3*F/R/T*y(3)*dydx(2);
end
The boundary condition is:
global C3A U_A E_A
res(1)=ya(1) - E_A;
res(2)=ya(2) - U_A;
res(3)=ya(3) - C3A;
end
I could not find the problem causing the error message here. Please help me out. Thank you very very very much!
7 Comments
Torsten
on 18 Dec 2017
If all boundary conditions are given at x=0, you should use ODE45 instead of BVP4C.
Best wishes
Torsten.
Luka
on 18 Dec 2017
Torsten
on 18 Dec 2017
How can you get an error message concerning "bvp4c" if you use "ode45" ?
Best wishes
Torsten.
Luka
on 18 Dec 2017
Luka
on 18 Dec 2017
Torsten
on 18 Dec 2017
Please insert the complete code you are using.
Best wishes
Torsten.
Luka
on 18 Dec 2017
Accepted Answer
Torsten
on 18 Dec 2017
global F R T e0 Z_Na Z_Cl Z_OH Z_H C_Na_A C_Cl_A C_OH_A C_H_A Cond_A i U_A E_A FixedCharge LA D_H D_OH e_r e
% Constants
F=96485; % Faradays number [C/mol]
R=8.3144621; % Ideal gas constant [J/(mol K)]
T=25+273.15; % Temperature [K]
e0=8.85419E-12; % Vacuum permittivity [F/m] or [C V/m]
Z_Na=1;
Z_Cl=-1;
Z_OH=-1;
Z_H=1;
C_Na_A=1.1*1000; % Conc of Na+ [mol/m3]
C_Cl_A=(1+10E-13))*1000; % Conc of Cl- [mol/m3]
C_OH_A=0.1*1000; % Conc of OH- [mol/m3]
C_H_A=10^(-13)*1000;
Cond_A=10.1; % Conductivity of the analyte [S/m] or [1/ohm m]
i=500; % current density [A/m2]
U_A=2; % Potential
E_A=i/Cond_A;
FixedCharge=1.71*1000;
LA=1.5E-4;
D_H=5.94/10^10;
D_OH=3.47/10^10;
e_r=20;
e=e_r*e0;
y0=[E_A;U_A;C_OH_A];
[T Y]=ode45(@odefun,linspace(0,LA,10),y0);
figure(1)
plot(T,Y(:,1),'g') %plot electric field
figure(2)
plot(T,Y(:,2),'g') %plot electric potential
figure(3)
plot(T,Y(:,3),'g') %plot conc OH-
function [dydx]= odefun(x,y)
global F R T Z_Na Z_Cl Z_OH Z_H C_Na_A C_Cl_A i U_A FixedCharge D_OH e D_H C_Na C_Cl
dydx=zeros(3,1);
C_Na=C_Na_A.*exp(-Z_Na.*F/(R*T).*(y(2)-U_A));
C_Cl=C_Cl_A.*exp(-Z_Cl.*F/(R*T).*(y(2)-U_A));
dydx(1)=F/e*(C_Na-C_Cl+10^(-8)/y(3)-y(3)+FixedCharge);
dydx(2)=-y(1);
J_H=-D_H*(-10^(-8)/y(3)^2*dydx(3)+Z_H*F/R/T*10^(-8)/y(3)*dydx(2));
dydx(3)=-((i/F-Z_H*J_H)/Z_OH/D_OH)-Z_OH*F/R/T*y(3)*dydx(2);
end
Note that you use dydx(3) in the calculation of J_H before it is defined. You will have to correct this error first.
Best wishes
Torsten.
More Answers (0)
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
