Solve differential equations with bvp4c

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!

[Torsten]

If all boundary conditions are given at x=0, you should use ODE45 instead of BVP4C.

Best wishes

Torsten.

[Luka]

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]

How can you get an error message concerning "bvp4c" if you use "ode45" ?

Best wishes

Torsten.

[Luka]

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]

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]

Please insert the complete code you are using.

Best wishes

Torsten.

[Luka]

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

Answer 1

[Torsten]

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.

[Luka]

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