Second order ODE - Error in using ODE45?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/191312/image.png)
Hello, I'm trying to solve the second order differential equation (Fick's second law) with a reaction term (which is a function of C_L) in MATLAB.
I created one script setting up my equations (odefcn.m) and another script to call to call odefcn (solodefcn.m).
function dC_Ldx=odefcn(x,C_L)
dC_Ldx = zeros(2,1);
syms N_r N_R N_r2 N_rR N_pR N_R2 R_L C_L
%Diffusion coefficient
D_ij= 1*10^-6;
%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;
%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
%K_5 = 0.077;
%K_6 = 0.000818;
K_i = 0.139;
%
N_T = 1.7;
N_r = (-(1+(C_L./K_2))+sqrt((1+(C_L./K_2)).^2-((4.*((2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3)))))).*N_T)))./...
((4./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3)))));
N_R = (C_L./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((C_L.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;
R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*C_L.*N_r))+...
((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
(d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));
R_L = vectorize(R_L);
dC_Ldx(1) = C_L(2);
dC_Ldx(2) = -R_L/D_ij;
and (solodefcn.m)
clear all; close all; clc;
x0=0;
xf=100;
C_L0 = [1 0];
[x,C_L] = ode45(@odefcn, [x0,xf], C_L0);
plot(x,C_L(:,1),x,C_L(:,2))
I keep on getting the errors:
Error using symengine (line 58) Index exceeds matrix dimensions.
Error in sym/subsref (line 696) B = mupadmex('symobj::subsref',A.s,inds{:});
Error in odefcn (line 45) dC_Ldx(1) = C_L(2);
Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in solodefcn (line 6) [x,C_L] = ode45(@odefcn, [x0,xf], C_L0);
and am unsure what I am doing wrong.
Am I setting up my second order differential equation incorrectly? I would appreciate any help I can get. Thank you!
0 commentaires
Réponses (1)
Walter Roberson
le 17 Juin 2018
Your declaration
syms N_r N_R N_r2 N_rR N_pR N_R2 R_L C_L
is overriding the C_L parameter name you are using . When you
syms C_L
that is the same as
C_L = sym('C_L')
so C_L ends up being a symbolic scalar, not something of length 2.
But even if you were preserving it correctly, what you are passing in is C_L(:,2) which is a scalar itself, so it would not have two elements inside ode function
6 commentaires
Walter Roberson
le 17 Juin 2018
Hmmm, I think,
function dC_Ldx=odefcn(x, C_L_derivs)
C_L = C_L_derivs(1);
....
dC_Ldx(1) = C_L_derivs(1);
dC_Ldx(2) = -R_L/D_ij;
Voir également
Catégories
En savoir plus sur Symbolic Math Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!