Second order ODE - Error in using ODE45?

2 次查看（过去 30 天）
Jennifer Yang
Jennifer Yang 2018-6-17
评论： Walter Roberson 2018-6-17
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!

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

回答（1 个）

Walter Roberson
Walter Roberson 2018-6-17
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 个评论
显示 4更早的评论
Jennifer Yang
Jennifer Yang 2018-6-17
Sorry. I'm still a bit confused. Is this for the original code where I had
dC_Ldx(2)=-R_L/D_ij
or for the situation where I had
dC_Ldx(2) = ((51.*C_L.*(((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1)./((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901) + (1146198855342281490208.*C_L.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(39069471042761828759.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) + (576460752303423488.*C_L.^2.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(6188225981639695.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) + (100.*C_L.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1))./(363.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901)) - 17./10))./20000 - (10487338329099335857321.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1))./(4611686018427387904000.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901)) + (409.*C_L.*(((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1)./((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901) + (100.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(901.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) + (100000.*C_L.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(27931.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) + (576460752303423488.*C_L.^2.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(6188225981639695.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) + (100.*C_L.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1))./(363.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901)) - 17./10))./50 - (10487338329099335857321.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(41551291026030765015040.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) - (333636569133360832851632851.*C_L.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(5000892293473514081152000.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) - (64927927453439194281.*C_L.^2.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1).^2)./(49505807853117560000.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901).^2) - (6554620466871470812811417.*C_L.*((2000.*C_L)./31 + (((2000.*C_L)./31 + 1).^2 - (160000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./1643 - 80./53).^(1./2) - 1))./(10462762654307136307200000.*((800000.*C_L.*((1000000.*C_L)./76867 + 1139./139))./27931 + 400./901)) + 178284751594688709574457./46116860184273879040000)./D_ij;
So what is the C = C_L(1); and L=C_L(2); supposed to do? Is this supposed to replace or fix the mismatch vector length for dC_Ldx(1) and dC_Ldx(2)?
You also mentioned to write the code in terms of C and L. I only really have everything in terms of C_L. So I'm not too sure what you mean.
Walter Roberson
Walter Roberson 2018-6-17
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;

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by