Second order ODE - Error in using ODE45?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Walter Roberson
am 17 Jun. 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 Kommentare
Walter Roberson
am 17 Jun. 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;
Siehe auch
Kategorien
Mehr zu Symbolic Math Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!