Solving a system of ODEs using ode45 in Matlab

Hello guys,

I have a set of four ODEs and need to solve them using ode45 in Matlab. In this ODEs, one of the parameters is a data matrix and other two parameters are vectors. The ODEs model needs to be repeated several time based on the second dimenstion of the data matrix. I tried to code this but getting an error saying that one of the ODEs is undefined once running the code.

The error I got is:

Undefined function or variable 'dmdt'.

Error in testmodel>model (line 120)

dxdt=[dmdt,dfdt,dpdt,dgdt];

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in testmodel (line 74)

[T,X] = ode45(@model,tspan,InitialCond);

close all; clear all; clc;

Kg = 0.1;
Kd = 1.1;
kg = 0.05;
kd = 0.05;
ka = 0.0455;
Kr = 4.5;
Ks = 1;
alpha = 300;
beta = 0.8;
w = 4;
l = 4;

m16= 7.26; 
m15 = m16; m14 = m15; m13 = m14; m12 = m13; m11 = m12; 
m10 = m11; m9 = m10; m8 = m9; m7 = m8; m6 = m7; m5 = m6; 
m4 = m5; m3 = m4; m2 = m3; m1 = m2;

f16= 0.465; 
f15 = f16; f14 = f15; f13 = f14; f12 = f13; f11 = f12; 
f10 = f11; f9 = f10; f8 = f9; f7 = f8; f6 = f7; f5 = f6; 
f4 = f5; f3 = f4; f2 = f3; f1 = f2;

p16= 0.016; 
p15 = p16; p14 = p15; p13 = p14; p12 = p13; p11 = p12; 
p10 = p11; p9 = p10; p8 = p9; p7 = p8; p6 = p7; p5 = p6; 
p4 = p5; p3 = p4; p2 = p3; p1 = p2;

g16= 0.127; 
g15 = g16; g14 = g15; g13 = g14; g12 = g13; g11 = g12; 
g10 = g11; g9 = g10; g8 = g9; g7 = g8; g6 = g7; g5 = g6; 
g4 = g5; g3 = g4; g2 = g3; g1 = g2;

ds= rand(1473,16);
sp= rand(1,16);
sg= rand(1,16);

d1 = ds(:,1) 
d2 = ds(:,2) 
d3 = ds(:,3) 
d4 = ds(:,4) 
d5 = ds(:,5) 
d6 = ds(:,6) 
d7 = ds(:,7) 
d8 = ds(:,8) 
d9 = ds(:,9) 
d10 = ds(:,10) 
d11 = ds(:,11) 
d12 = ds(:,12) 
d13 = ds(:,13) 
d14 = ds(:,14) 
d15 = ds(:,15) 
d16 = ds(:,16) 

rd = [d1,d2 ,d3 ,d4 ,d5 ,d6 ,d7 ,d8 ,d9 ,d10,d11, d12,d13,d14,d15,d16];

m= [m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16];
f= [f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16];
p= [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16];
g= [g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13,g14,g15,g16];

InitialCond= [m,f,p,g];
Tstart= 0;
StepSize= 0.5;
Tend= 470000;
tspan=(Tstart:StepSize:Tend);

[rds cds]=size(ds);

for ii=1:rds
    for jj=1:cds   
    [T,X] = ode45(@model,tspan,InitialCond);

figure(1)
plot(T, X(:,1));
hold on

figure(2)
plot(X(:,2),X(:,1))

hold on
   end
end

function dxdt = model(t,Y)
format long
global m
global f
global p
global g
global rd
global sp
global sg

n= 16

mn = Y(1:n);
fn = Y(n+1:2*n);
pn = Y(2*n+1:3*n);
gn = Y(3*n+1:4*n);

[rim_r,rim_c]= size(rd)

                for i=1:length(rim_r)
                    for j=1:length(rim_c)
                        for ii=1:length(sp)
                            for jj=1:length(sg)
                                dmdt = kg*(fn*gn)/(Kg+fn)-(mn*pn)/(1+mn)+rd(i,j)
                                dfdt = -kg*(fn*gn)/(Kg+fn)+(mn*pn)/(1+mn)-(fn*pn)/(1+fn)
                                dpdt = sp(ii)*alpha*(fn^w)/(Kr^w+fn^w)-ka*pn^2
                                dgdt = sg(jj)*beta*(fn^l)/(Ks^l+fn^l)-kd*(gn*pn)/(Kd+gn)
                            end
                        end
                    end
                    
                end
                
dxdt=[dmdt,dfdt,dpdt,dgdt];
end

Can you please help me with fixing my code.

[Cris LaPierre]

Without your code, data and equations, there's very little we can do to help you. Best I can think to do is point you to the ode45 examples in the documentation.

[James Tursa]

What have you done so far? Can you post your code?

[Shy]

Hi Cris and James,

Thanks for your reply. Here below is my code:

The error I got is:

Undefined function or variable 'dmdt'.

Error in testmodel>model (line 120)

dxdt=[dmdt,dfdt,dpdt,dgdt];

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in testmodel (line 74)

[T,X] = ode45(@model,tspan,InitialCond);

close all; clear all; clc;

Kg = 0.1;
Kd = 1.1;
kg = 0.05;
kd = 0.05;
ka = 0.0455;
Kr = 4.5;
Ks = 1;
alpha = 300;
beta = 0.8;
w = 4;
l = 4;

m16= 7.26; 
m15 = m16; m14 = m15; m13 = m14; m12 = m13; m11 = m12; 
m10 = m11; m9 = m10; m8 = m9; m7 = m8; m6 = m7; m5 = m6; 
m4 = m5; m3 = m4; m2 = m3; m1 = m2;

f16= 0.465; 
f15 = f16; f14 = f15; f13 = f14; f12 = f13; f11 = f12; 
f10 = f11; f9 = f10; f8 = f9; f7 = f8; f6 = f7; f5 = f6; 
f4 = f5; f3 = f4; f2 = f3; f1 = f2;

p16= 0.016; 
p15 = p16; p14 = p15; p13 = p14; p12 = p13; p11 = p12; 
p10 = p11; p9 = p10; p8 = p9; p7 = p8; p6 = p7; p5 = p6; 
p4 = p5; p3 = p4; p2 = p3; p1 = p2;

g16= 0.127; 
g15 = g16; g14 = g15; g13 = g14; g12 = g13; g11 = g12; 
g10 = g11; g9 = g10; g8 = g9; g7 = g8; g6 = g7; g5 = g6; 
g4 = g5; g3 = g4; g2 = g3; g1 = g2;

ds= rand(1473,16);
sp= rand(1,16);
sg= rand(1,16);

d1 = ds(:,1) 
d2 = ds(:,2) 
d3 = ds(:,3) 
d4 = ds(:,4) 
d5 = ds(:,5) 
d6 = ds(:,6) 
d7 = ds(:,7) 
d8 = ds(:,8) 
d9 = ds(:,9) 
d10 = ds(:,10) 
d11 = ds(:,11) 
d12 = ds(:,12) 
d13 = ds(:,13) 
d14 = ds(:,14) 
d15 = ds(:,15) 
d16 = ds(:,16) 

rd = [d1,d2 ,d3 ,d4 ,d5 ,d6 ,d7 ,d8 ,d9 ,d10,d11, d12,d13,d14,d15,d16];

m= [m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16];
f= [f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16];
p= [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16];
g= [g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13,g14,g15,g16];

InitialCond= [m,f,p,g];
Tstart= 0;
StepSize= 0.5;
Tend= 470000;
tspan=(Tstart:StepSize:Tend);

[rds cds]=size(ds);

for ii=1:rds
    for jj=1:cds   
    [T,X] = ode45(@model,tspan,InitialCond);

figure(1)
plot(T, X(:,1));
hold on

figure(2)
plot(X(:,2),X(:,1))

hold on
   end
end

function dxdt = model(t,Y)
format long
global m
global f
global p
global g
global rd
global sp
global sg

n= 16

mn = Y(1:n);
fn = Y(n+1:2*n);
pn = Y(2*n+1:3*n);
gn = Y(3*n+1:4*n);

[rim_r,rim_c]= size(rd)

                for i=1:length(rim_r)
                    for j=1:length(rim_c)
                        for ii=1:length(sp)
                            for jj=1:length(sg)
                                dmdt = kg*(fn*gn)/(Kg+fn)-(mn*pn)/(1+mn)+rd(i,j)
                                dfdt = -kg*(fn*gn)/(Kg+fn)+(mn*pn)/(1+mn)-(fn*pn)/(1+fn)
                                dpdt = sp(ii)*alpha*(fn^w)/(Kr^w+fn^w)-ka*pn^2
                                dgdt = sg(jj)*beta*(fn^l)/(Ks^l+fn^l)-kd*(gn*pn)/(Kd+gn)
                            end
                        end
                    end
                    
                end
                
dxdt=[dmdt,dfdt,dpdt,dgdt];
end

Answer 1

[Walter Roberson]

You do not initialize any of your global variables, so when MATLAB sees the global statements it creates the global variables and initializes them to empty. But length(empty) is 0, so your for loops do not do any work, and so do not assign to any of the variables you need to create dxdt .

Best thing to do is not use global. See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Side note: you need to return a column vector of values, not a row vector.

Also, you are overwriting all of your variables in every iteration through the loops.

Note that you need to produce one output element for each element of the boundary conditions, which is 64 elements. But you have a nested loop that is 4 deep, so you are at risk of producing too many outputs.

[rim_r,rim_c]= size(rd)

rim_r and rim_d are going to be scalars

                for i=1:length(rim_r)
                    for j=1:length(rim_c)

length() of a scalar is 1.

[Shy]

Thanks Walter. I have initialised the global variables but still getting the same error. Also, I adjusted the length(rim_r) and length(rim_c) into size(rd,1) and size(rd,2) respectively. But getting the same error. Can you please help me more in this? I apprecited your help.

[Walter Roberson]

How exactly do you initialize your global variables?

Note that in your script you initialize variables that have the same name as the global variables but which are not global. A variable is only global inside a workspace that uses "global"

[Shy]

@Walter: can you please help me? I am stuck and lost with the code. Appreciate your help

[Walter Roberson]

How exactly do you initialize your global variables?

[Walter Roberson]

I recommend that you rewrite the code without global variables.

[Walter Roberson]

P.Kg = 0.1;
P.Kd = 1.1;
P.kg = 0.05;
P.kd = 0.05;
P.ka = 0.0455;
P.Kr = 4.5;
P.Ks = 1;
P.alpha = 300;
P.beta = 0.8;
P.w = 4;
P.l = 4;

m16= 7.26; 
m15 = m16; m14 = m15; m13 = m14; m12 = m13; m11 = m12; 
m10 = m11; m9 = m10; m8 = m9; m7 = m8; m6 = m7; m5 = m6; 
m4 = m5; m3 = m4; m2 = m3; m1 = m2;

f16= 0.465; 
f15 = f16; f14 = f15; f13 = f14; f12 = f13; f11 = f12; 
f10 = f11; f9 = f10; f8 = f9; f7 = f8; f6 = f7; f5 = f6; 
f4 = f5; f3 = f4; f2 = f3; f1 = f2;

p16= 0.016; 
p15 = p16; p14 = p15; p13 = p14; p12 = p13; p11 = p12; 
p10 = p11; p9 = p10; p8 = p9; p7 = p8; p6 = p7; p5 = p6; 
p4 = p5; p3 = p4; p2 = p3; p1 = p2;

g16= 0.127; 
g15 = g16; g14 = g15; g13 = g14; g12 = g13; g11 = g12; 
g10 = g11; g9 = g10; g8 = g9; g7 = g8; g6 = g7; g5 = g6; 
g4 = g5; g3 = g4; g2 = g3; g1 = g2;

ds = rand(1473,16);

%why do you split ds into columns and then re-assemble exactly the
%same columns into rd ? Why not just rd = ds ?

d1 = ds(:,1) 
d2 = ds(:,2) 
d3 = ds(:,3) 
d4 = ds(:,4) 
d5 = ds(:,5) 
d6 = ds(:,6) 
d7 = ds(:,7) 
d8 = ds(:,8) 
d9 = ds(:,9) 
d10 = ds(:,10) 
d11 = ds(:,11) 
d12 = ds(:,12) 
d13 = ds(:,13) 
d14 = ds(:,14) 
d15 = ds(:,15) 
d16 = ds(:,16) 

P.rd = [d1, d2 ,d3 ,d4 ,d5 ,d6 ,d7 ,d8 ,d9 ,d10, d11, d12, d13, d14, d15, d16];

P.sp = rand(1,16);
P.sg = rand(1,16);

P.m = [m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16];
P.f = [f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16];
P.p = [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16];
P.g = [g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13,g14,g15,g16];

InitialCond= [P.m, P.f, P.p, P.g];
Tstart= 0;
StepSize= 0.5;
Tend= 470000;
tspan=(Tstart:StepSize:Tend);

[rds cds]=size(ds);

for ii=1:rds
    for jj=1:cds   
    [T,X] = ode45(@(t,y)model(t,y,P), tspan, InitialCond);

    figure(1)
    plot(T, X(:,1));
    hold on
    
    figure(2)
    plot(X(:,2),X(:,1))
    
    hold on
   end
end

function dxdt = model(t, Y, P)
    format long
    m = P.m;
    f = P.f;
    p = P.p;
    g = P.g;
    rd = P.rd;
    sp = P.sp;
    sg = P.sg;
    
    n = 16;
    
    mn = Y(1:n);
    fn = Y(n+1:2*n);
    pn = Y(2*n+1:3*n);
    gn = Y(3*n+1:4*n);
    
    [rim_r,rim_c] = size(rd);
    
    for i=1:rim_r
        for j=1:rim_c
            for ii=1:length(sp)
                for jj=1:length(sg)
                    dmdt = P.kg .* (fn.*gn) ./ (P.Kg+fn) - (mn.*pn) ./ (1+mn) + rd(i,j);
                    dfdt = -P.kg .* (fn.*gn) ./ (P.Kg+fn) + (mn.*pn) ./ (1+mn) - (fn.*pn) ./ (1+fn);
                    dpdt = sp(ii) .* P.alpha .* (fn.^P.w) ./ (P.Kr.^P.w + fn.^P.w) - P.ka .* pn.^2;
                    dgdt = sg(jj) .* P.beta .* (fn.^P.l) ./ (P.Ks.^P.l + fn.^P.l) - P.kd .* (gn.*pn) ./ (P.Kd + gn);
                end
            end
        end
    end
                    
    dxdt=[dmdt;dfdt;dpdt;dgdt];
end

[Walter Roberson]

Notice that your mn and fn and pn and gn are all vectors, and in your nested for loop that you are using the entire vectors in each case, never indexing them. So each of the four variables comes out the right size when you assign to them... but each time you are overwriting the entire variable, so you might as well only do the last iteration (if overwriting was your intent.)

You do index rd and sp and sg, but it is not clear why.

[Shy]

Thanks Walter for your help and answer. I really appreciate your nice doing.

In my code, I have four ODEs and I need to iterated the ODEs model 16 times. This is because the number of columns of rd=ds is equal to 16. In each time, I need to feed in the ODEs model with d1 at iteration 1, and then with d2 at iteration 2 and so on till d16. Finally the ODEs solution should be the updated outcome after all 16 iterations (after feeding in the ODEs with all 16 columns of ds or rd).

As you see but I am not sure if I have done it the correct way, I have mn, fn, pn and gn variables as vectors with length of 16 each. I made them as it is because I need to find the ODEs solution (m1,m2,...m16), (f1,f2,...f16), (p1,p2,...p16), and (g1,g2,...g16). So, m1 is the solution of ODEs at iteration 1 and m16 is the solution of ODEs at iteration 16. The same applies for f, p and g. So imagine, we have four ODEs without iteration, the solution should be (4 by 1) vector. But, if we have four ODEs with 16 iteration, the solution vector should be (4 by 16), if I am not mistaken!

Also, in my ODEs, I have sp and sg as vectors of length of 16.

Can you reflect on my code again? I really appreciate your help and time. The code you shared above is still running! It takes so long and still running.

[Walter Roberson]

You lost me. rd columns might explain your for j loop, but none of the other loops, and even then would require indexing of results locations. I do not understand at all about the iterations.

[Shy]

@Walter: the other for-loops are established for sp and sg vectors. Can you help me with the code? Thanks

[Walter Roberson]

I need to iterated the ODEs model 16 times

I need to feed in the ODEs model with d1 at iteration 1, and then with d2 at iteration 2 and so on till d16

Each iteration you are working with a vector of 16 elements, so with 16 iterations you would get 16*16 = 256 results at least -- unless, that is, you are initializing some variable before the first round, and the results of the first round get assigned to that variable, and the results are used to calculate values for the second round, which is used to update the variable, and so on, so at the end you get 16 results that incorporate information from all 16 rounds. For example in some cases it might make sense to total the results over 16 round.

Unless you do have some kind of totalling or feedback, I am unable to match up the number of iterations you are doing with the number of output values you are expecting.

Not unless where you have (for example)

dmdt = P.kg .* (fn.*gn) ./ (P.Kg+fn) - (mn.*pn) ./ (1+mn) + rd(i,j);

that the vector variables being used should be indexed so that you are producing a scalar output each time ??