Solve differential equations with bvp4c

2017 12 17
1 个回答
回答已采纳
更新时间：2017 12 18
11 次查看（30 天）

0 个投票

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 个评论
显示 5更早的评论 隐藏 5更早的评论

Torsten
Torsten 2017-12-18
If all boundary conditions are given at x=0, you should use ODE45 instead of BVP4C.
Best wishes
Torsten.
Luka
Luka 2017-12-18
Thank you very much Torsten. Follow your advice, I've tried ODE45, and simplified one of the equations as below:
y0=[E_A;U_A;C3A];
xspan=[0 10^(-4)];
sol=ode45(@odefun,xspan,y0);
function [dydx]= odefun(x,y)
global F R T Z1 Z2 Z3 C1A C2A i U_A FixedCharge D3 e
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);
dydx(3)=-i/Z3/D3-Z3*F/R/T*y(3)*dydx(2);
end
I got the error message of
Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in ODEtest (line 37)
sol=bvp4c(@odefun,@odebc,solinit);
Could you please give me more advice about what's happening here? Really appreciate your help.
Torsten
Torsten 2017-12-18
How can you get an error message concerning "bvp4c" if you use "ode45" ?
Best wishes
Torsten.
Luka
Luka 2017-12-18
I'm so sorry I messed up here. The error message for ode45 is
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in ode_Dec1717 (line 38)
sol=ode45(@odefun,xspan,y0);
Luka
Luka 2017-12-18
When I simplified one of the equation and used bvp4c, the error message I got is
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in ode_Dec1717 (line 38)
sol=ode45(@odefun,xspan,y0);
Torsten
Torsten 2017-12-18
Please insert the complete code you are using.
Best wishes
Torsten.
Luka
Luka 2017-12-18
clc
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;
solinit =bvpinit(linspace(0,LA,10),[E_A U_A C_OH_A]);
sol=bvp4c(@odefun,@odebc,solinit);
xint=linspace(0,LA);
Sxint=deval(sol,xint);
figure(1)
plot(xint,Sxint(1,:),'g') %plot electric field
figure(2)
plot(xint,Sxint(2,:),'g') %plot electric potential
figure(3)
plot(xint,Sxint(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
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
function res= odebc(ya,yb)
global C_OH_A U_A E_A
res=[ya(1)-E_A;ya(2)-U_A;ya(3)-C_OH_A];
end

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Torsten
Torsten 2017-12-18

0 个投票

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.

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Luka
Luka 2017-12-18
Thank you so much Torsten!
Now I can get results for x below 1e-10. For x higher than 1e-10, it says ’Unable to solve the collocation equations -- a singular Jacobian encountered.‘ I might need to fix my equations also.
I really appreciate your help!
Best,
Luka

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by