Using ODE45 with a matrix as input to my function

I want to repeatedly solve a set of ODEs for different values of some parameters, which will be encoded in a matrix. The equation I want to solve is

where Γcontains the parameters I wish to vary.

I made a function file:

function dxdt = UseGamma(t,x,Gamma)
Gamma
y = Gamma*x;
dxdt = zeros(size(x));
dxdt(1) = y(1);
dxdt(2) = y(2);
dxdt(3) = y(3);
dxdt(4) = y(4); 
dxdt(5) = y(5) ;

If I create the variables tspan, xinit, Gamma of appropriate size and call the function, I get:

UseGamma(tspan,xinit,Gamma)

Gamma =

   -0.7500         0   -0.2500    1.0000    0.2500
         0   -0.2500   -0.2500    1.0000    0.7500
   -0.1250   -0.1250   -1.0000         0         0
         0         0         0   -2.0000         0
    0.7500    0.2500    0.5000         0   -1.0000


ans =

   -0.7400
    0.0075
   -0.1237
         0
    0.7325

However when I try to use this function as the RHS of my differential equation by calling it with ode45, I get:

    [t,x] = ode45('UseGamma', tspan, xinit, Gamma);

Gamma =

  0×0 empty char array

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise
multiplication, use '.*'.

Error in UseGamma (line 3)
y = Gamma*x;

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);
  

Is it possible to use the matrix as the input to ode45? If so, what am I doing wrong? If not, what is the simplest way to input these 25 values in Gamma to ode45? Thanks to everybody who is always so helpful!

Answer 1

[Jan]

You can abbreviate:

function dxdt = UseGamma(t,x,Gamma)
Gamma
y = Gamma*x;
dxdt = zeros(size(x));
dxdt(1) = y(1);
dxdt(2) = y(2);
dxdt(3) = y(3);
dxdt(4) = y(4); 
dxdt(5) = y(5) ;

to

function dxdt = UseGamma(t,x,Gamma)
dxdt = Gamma * x;

This style of calling an ODE integrator is outdated for over 20 years now (since Matlab R6.5):

[t,x] = ode45('UseGamma', tspan, xinit, Gamma);

In Matlab 2022b this is treated as error and I get a message about an unknown function "UseGamma" when I run this in the forum. Prefer:

[t,x] = ode45(@(t,x) UseGamma(t, x, Gamma), tspan, xinit);

But the main problem is, that somewhere Gamma is defined as empty CHAR ''. Please post a copy of the code.

[Jan]

Gamma = [   -0.7500         0   -0.2500    1.0000    0.2500
         0   -0.2500   -0.2500    1.0000    0.7500
   -0.1250   -0.1250   -1.0000         0         0
         0         0         0   -2.0000         0
    0.7500    0.2500    0.5000         0   -1.0000];
tspan = [0, 1];
xinit = rand(5, 1);

% [t,x] = ode45('UseGamma', tspan, xinit, Gamma); % Failing with message:

Unrecognized function or variable 'UseGamma'.

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Working correctly:
[t,x] = ode45(@(t,x) UseGamma(t, x, Gamma), tspan, xinit);

function dxdt = UseGamma(t,x,Gamma)
dxdt = Gamma * x;
end

[Emily]

Thanks Jan,

I had originally used your shorter version of the function, but had been trying the longer more explicit version to troubleshoot.

Yes, the problem seems to be the way I was calling the integrator. I was following an example on the internet. Here is a code that works:

% A MWE calling UseGamma in ODE45

ParameterA = 1
ParameterB = 10
sum = ParameterA + ParameterB

Gamma=zeros(5,5);
Gamma(1,1) = -(1/2)*(1+ParameterA^2/sum^2);
Gamma(1,3) = -(ParameterB*ParameterA/(2*sum^2));
Gamma(1,4) = 1;
Gamma(1,5) = ParameterB^2/(2*sum^2);
Gamma(2,2) = -1*ParameterB^2/(2*sum^2);
Gamma(2,3) = -ParameterA*ParameterB/(2*sum^2);
Gamma(2,4) = 1;
Gamma(2,5) = (1/2)*(1+ParameterA^2/sum^2);
Gamma(3,1) = -ParameterA*ParameterB/(4*sum^2);
Gamma(3,2) =  -ParameterA*ParameterB/(4*sum^2); 
Gamma(3,3) = -1;
Gamma(4,4) = -2;
Gamma(5,5) = -1;
Gamma(5,1) = (1/2)*(1+ParameterA^2/sum^2);
Gamma(5,2) = 1*ParameterB^2/(2*sum^2);
Gamma(5,3) = ParameterB*ParameterA/sum^2;

xinit = [0.99; 0; 0;0;0.01];
tspan = [0 50];

UseGamma(tspan,xinit,Gamma)
[t,x] = ode45(@(t,x) UseGamma(t, x, Gamma), tspan, xinit);

function dxdt = UseGamma(t,x,Gamma)
dxdt = Gamma * x;
end

However when I call it using the old format, it doesn't work:

[t,x] = ode45('UseGamma', tspan, xinit, Gamma);

Gamma =

  0×0 empty char array

[Jan]

@Emily: Fine.

Just a note: Using the names of built-in functions as variables is a frequent cause of bugs. See:

sum(1:10)  % Working

ans = 55

sum = 17;
sum(1:10)  % Error:

Index exceeds the number of array elements. Index must not exceed 1.

'sum' appears to be both a function and a variable. If this is unintentional, use 'clear sum' to remove the variable 'sum' from the workspace.