Second order ODE - Error in using ODE45?

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!

Answer 1

[Walter Roberson]

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

[Jennifer Yang]

So would me having a separate script and manually plugging in the answer to vectorize(R_L) help prevent this error?

Where the script is now

function dC_Ldx=odefcn(x,C_L)

dC_Ldx = zeros(2,1);

%Diffusion coefficient
D_ij= 1*10^-6;

dC_Ldx(1) = C_L(2);

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;

Although in this case I now get the error:

In an assignment A(I) = B, the number of elements in B and I must be the same.

Error in odefcn (line 54) dC_Ldx(2) = ((51.*C_L.*(((2000.*C_L...

I'm not quite sure what's the right way to set up the plot as well as whether I'm setting up the second order differential equation correctly.

[Walter Roberson]

When your C_L is a vector of length 2, then C_L is a vector on your dC_Ldx(2) line, so your would be calculating a vector of results and trying to store the vector into the single location.

[Walter Roberson]

Note: I was wrong about it being a scalar being passed into the function; it is a vector of length 2.

If you delete the "syms" call, then everything works okay numerically until the last line, which has a problem because R_L is a vector of length 2. Probably in the places where you used C_L in the code, you should have indexed according to which of the two entries you wanted.

Consider putting something like

C = C_L(1);
L = C_L(2);

at the top of your odefcn and then writing the code in terms of C and L. (With no syms)

[Walter Roberson]

Okay, I see that separating C_L into C and L is not appropriate, but you can do

CL = C_L(1);
d_CL = C_L(2);

and then refer to the appropriate variable in the appropriate formula.

I also suggest that you have a look at https://www.mathworks.com/help/symbolic/odefunction.html which shows how to work through symbolic expressions of an ODE to get to a numeric solution.

[Jennifer Yang]

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]

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;