What is wrong with my use of ode45?

Good morning,

I am trying to numerically solve a system of ODEs from a research paper. I am new to MATLAB, so I am folloing along with this documentation (section titled "Pass Extra Parameters to ODE function"). I am failing miserably. I think one of my problem is not understanding the convoluted system of naming and calling and passing functions in MATLAB. The examples in the documentation don't seem to consistently use one method. If there is a reason for switching methods, it is not explained. I see functions defined with arguments called without them. I see functions with @ signs, and I see functions without @ signs. I see functions with two sets of parameters in two sets of parenthesis like @(t,y) func(t,y, otherstuff). I am having a hard time understanding why such an essential part of any programming language has been made this complex.

Anyway, end rant. Sorry.

I have created a .m file containing the following function definition:

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)
dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(3))) - bet*(y(1)^2/((alph*y(3))^2 + y(1)^2));
         r_e*y(2)*(1 - (y(2)/k_e)) - p*y(1)/y(3);
         r_s*y(3)*(1 - y(3)*k_e/(k_s*y(2)))
        ];
end

Then in the script I'm trying to write, I have the following code:

addpath('C:\Users\wesle\Documents\Math\564 Homework');

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.1;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

tspan = [0 500];
init = [1 1 1];

[t, y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);
plot(t,y(:,1),'-');

%define system
%dB = diff(b) == r_b*b*(1 - b/(k_b*s)) - bet*(b^2/(alph^2 + b^2));
%dE = diff(e) == r_e*e*(1 - e/k_e) - p*b/s;
%dS = diff(s) == r_s*s*(1 - (s*k_e)/(k_s*e));
%odes = [dB ; dE ; dS];

The current error I am recieving is: "Failure at t=1.209318e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.273737e-13) at time t."

[madhan ravi]

I am really impressed by your work so good work!

So coming to your errors:

1) In your function definition change your output as dYdt (MATLAB is case sensitive)

2) Just copy the input arguments in your function definition and paste it exactly in the input arguments when you call the function.

[Wesley Neill]

Sir,

Thank you. I saw those errors right before you posted. I have fixed them and updated my question with the current error that I am recieving:

"Failure at t=1.209318e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.273737e-13) at time t."

Regards,

Wes

[Steven Lord]

Others have addressed the specific ODE function concerns in more or less detail. I want to address some of the other concerns / questions / puzzles you called out.

I see functions defined with arguments called without them.

That can be correct. Some functions that are defined to take input arguments require those input arguments to function correctly. For example if you try to call the sin function with no input arguments, you will receive an error. MATLAB needs to know the angle whose sine you want to compute in order to return an answer.

On the other hand, some function define some default behavior to use if you don't specify the optional input arguments. An example of this is the eps function that can accept input arguments but iuses a default (x = 1) if you don't specify an input.

I see functions with @ signs, and I see functions without @ signs.

"functions with @ signs" makes something called a function handle. Function handles let you refer to and call a function without hard-coding a call to that specific function into your code. See this documentation page and/or the example below.

function y = evaluateFunctionHandleAtZero(fh)
y = fh(0);
end

If I wanted to call the sin function at x = 0, I can do that using the evaluateFunctionHandleAtZero function above.

y = evaluateFunctionHandleAtZero(@sin);

Note that nowhere in the evaluateFunctionHandleAtZero function do I make any reference to the sin function, so I can reuse that same function to evaluate the cos function simply by changing what I pass into it.

z = evaluateFunctionHandleAtZero(@cos);

Using a function name without the @ sign means you're calling it.

q = plus(sin(1), cos(2));

This calls the sin and cos functions to obtain numbers then passes those numbers into the plus function.

I see functions with two sets of parameters in two sets of parenthesis like @(t,y) func(t,y, otherstuff).

That defines a special type of function handle called an anonymous function. These are useful when you want to write a very simple function without going to the trouble of creating a separate function file to contain it. I could have written the evaluateFunctionHandleAtZero function above as an anonymous function rather than needing to create evaluateFunctionHandleAtZero.m.

evaluateFunctionHandleAtZero2 = @(fh) fh(0);
y2 = evaluateFunctionHandleAtZero2(@sin);
z2 = evaluateFunctionHandleAtZero2(@cos);

Anonymous functions are also used, as in the specific example you called out, as an adapter to bridge the gap between two functions with different expectations for signatures. ode45 will pass exactly two inputs into your ODE function. If your function needs a third input too, you need an adapter that will forward the two inputs received from ode45 and the third (the "otherstuff") from somewhere else into your ODE function.

When you specify @(t, y) func(t, y, otherstuff) ode45 sees a function that accepts two inputs (t and y) and it's happy. Your ODE function func receives three inputs (t, y, and otherstuff) and it's happy.

I am having a hard time understanding why such an essential part of any programming language has been made this complex.

Flexibility. In order to allow functions like the ODE solvers (that accept an unevaluated "reference" to a function rather than requiring the code author to specify the name of the function at write-time) to exist we need to be able to specify that unevaluated reference. Function handles are that unevaluated reference.

Answer 1

[Star Strider]

There were a few problems in your code that I fixed.

First, be very careful creating matrices such as you did in ‘odefcn’ because MATLAB interprets spaces as delimiting different columns. I suspect you would have gotten a concatenation error with that. I changed it slightly to eliminate that problem, however your matrix likely would have worked if you had put parentheses around each row, so MATLAB would consider them each one element.

Second, one of the additional arguments you are passing ‘r_k’ does not exist. I substitute ‘k_b’ for it here.

With those changes, your code runs:

function dYdt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)
dYdt = zeros(3,1);
dYdt(1) =  r_b*y(1)*(1 - y(1)/(k_b*y(3))) - bet*(y(1)^2/((alph*y(3))^2 + y(1)^2));
dYdt(2) =  r_e*y(2)*(1 - (y(2)/k_e)) - p*y(1)/y(3);
dYdt(3) =  r_s*y(3)*(1 - y(3)*k_e/(k_s*y(2)));
end

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.1;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;
tspan = [0 500];
init = [1 1 1];
[t, y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);
plot(t,y(:,1),'-');

[Star Strider]

@Wesley Neill — The error you are seeing is because your function generates a singularity (Inf or NaN) at that point.

[Wesley Neill]

Thank you for pointing out the typographical errors. As is normal for me, I did a forehead slap as soon as I posted because I noticed them myself.

As for the comment about matrices, I used semi-colons, exactly as was done in the documentation. I assumed that semi-colons delimit rows. Is that incorrect?

The code runs after my fixing of typographical errors, but it encounters an error after a certain number of time iterations. I have updated my question to reflect that error.

I will also try to run your code and see what happens!

Thank you very much for your answer.

[Star Strider]

My pleasure.

My code throws the same Warning (it is not an Error). The problem is that your function encounters a singularity at t=1.209318e+02. If it is supposed to be defined for the full 500 time units, check your code to be sure you coded your differential equations correctly.