Finding the maxima of a function using ODE45

Question

I'm trying to locate the locations of one of the equations in a system of differential equations in MATLAB.I'm trying to use the Events propety of odeset.How do I pick out the particular equation in my function?

options = odeset('Events',@event);
[t x tm xm ie] = ode45(@Lorenz,[0 200],I,options);


function X = Lorenz(t,x)
r = 15;
sigma = 10;
b = 8/3;
X(1,1) = sigma*(x(2,1)-x(1,1));
X(2,1) = r*(x(1,1)) - x(2,1) -x(1,1)*x(3,1);
X(3,1) = x(1,1)*x(2,1) - b*x(3,1);
end

function [value,isterminal,direction] = event(t,x)
value  = Lorenz(t,x); %I can't specify X(3,1) here
isterminal = 0;
direction = -1;
end

In particular I'm trying to record whenever X(3,1) = 0.

Thank you

Answer 1

If you're looking for maxima of an ODE, as the title of your question indicates, then you are very close. You're using the the roots of the differential equation itself to find these points, i.e., when the derivatives are zero. This is slightly different from the solution having zero (or some other) value, but related. The problem is that you're specifying value = Lorenz(t,x) and the ODE function returns a vector when you're only interested in x(3). But you have access to the state vector and have have three alternatives.

The simplest:

function [value,isterminal,direction] = event(t,x)
b = 8/3;
value = x(1)*x(2)-b*x(3); % The third equation from Lorenz(t,x) 
isterminal = 0;
direction = -1;

Or, less efficiently:

function [value,isterminal,direction] = event(t,x)
y = Lorenz(t,x); % Evaluate all three equations to get third one
value = y(3);
isterminal = 0;
direction = -1;

Or, if you want maxima for all three dimensions:

function [value,isterminal,direction] = event(t,x)
value = Lorenz(t,x);
isterminal = [0;0;0];
direction = [-1;-1;-1];

If you're interested in the global maximum then you'll need to process the outputs, xm. Or if you're in a regime where the system has certain oscillatory behaviors, then you may be able to just switch isterminal to 1 or just look at the first value in xm.

Lastly, you might consider passing your parameters via anonymous functions:

r = 15;
sigma = 10;
b = 8/3;
f = @(t,x)Lorenz(t,x,r,sigma,b);

I = [1;5;10];
options = odeset('Events',@(t,x)event(t,x,b));

[t,x,tm,xm,ie] = ode45(f,[0;10],I,options);

with:

function X = Lorenz(t,x,r,sigma,b)
...

function [value,isterminal,direction] = event(t,x,b)
...

Answer 2

Basically, looking at the documentation, if you're interested in looking at when x(3) = 0, then you need to rewrite your event function:

function [value,isterminal,direction] = event(t,x)
value  = x(3); %I can't specify X(3,1) here --> why not?? Lorenz(t,x) is going to return the differential. That's not what you want
isterminal = 0;
direction = 0; %The desired directionality should be "Detect all zero crossings"
end

Now I don't know how you defined I in

[t x tm xm ie] = ode45(@Lorenz,[0 200],I,options);

But your solution to the equation is very stable around a few points and you may only see one zero crossing if x(3) at time zero is negative.

[t x tm xm ie] = ode45(@khal,[0 5],[10 -1 -20],options);
tm =

    0.1085