Using ode45 for solution curve

Equations:

df/dt= 4f(t) - 3f(t)p(t)

dp/dt= -2p(t) + f(t)p(t)

Question: For the solution curve that starts at (1,1), how many units of time does it take before the curve returns to (1,1)? Try to get it within two decimal places. Repeat for the curve that starts at (0.5,0.5).

Here is the code that creates the solution curve:

gh = @(t,x) [4*x(1) - 3*x(1).*x(2); -2*x(2) + x(1).*x(2)];
tspan = linspace(0, 5, 250);
figure; hold on
for c = 0:0.1:5
    [t,x] = ode45(gh, tspan, [c; c]);
    plot(x(:,1), x(:,2))
end
axis tight

I was told there is a way to use ode45 to find how many units of time, but I am not sure how?

[Bob]

Below is the code that I created and it gives an answer but I am not sure if it is the correct answer?

Code:

gh = @(t,x) [4*x(1) - 3*x(1).*x(2); -2*x(2) + x(1).*x(2)];
[t,x] = ode45(gh, tspan, [1; 1]);
t(end)
x(end)

Answer 1

[Star Strider]

These are getting more challenging (and interesting).

One way is to use an 'Event' function. (You can either create a separate function file (without input arguments) and put all of these in it and then run the function from the Command Window, or save ‘EventFcn’ to a separate .m file (as ‘EventFcn.m’ and run the ODE integration from its own script file.)

The code:

      function [value,isterminal,direction] = EventFcn(t,x)
      value = [(x(1) - 1); (x(2) - 1)];
      isterminal = [0; 0];
      direction = [0; 0];
      end

gh = @(t,x) [4*x(1) - 3*x(1).*x(2); -2*x(2) + x(1).*x(2)];
tspan = linspace(0, 5, 250);
c = 1;
opts = odeset('Events', @EventFcn);
[t,x,EventTime,EventXVal,ie] = ode45(gh, tspan, [c; c], opts);
% EventTime                                                           % Show Values (Delete)
% EventXVal                                                           % Show Values (Delete)

x_1_1 = abs(diff(EventXVal,1,2)) < 1.1;                             % Logical Vector Of Approximately Equal ‘x’ Values
EqualTimes = EventTime(x_1_1)                                       % Corresponding Event Times
UniqueTimes = diff([0; EqualTimes]) >= 0.2                          % Use ‘diff’ To Select Unique Times
ReturnTimes = EqualTimes(UniqueTimes)                               % Desired Result

When I run it, I get:

ReturnTimes =

       1.7980
       2.2888
       4.0924
       4.5863

Set ‘c = 0.5;’ for the second run to complete the exercise.

[Bob]

Thank you so much. I have one question though when you get the return times which one is the actual answer out of the 4 possible numbers that are printed out?

[Bob]

Also getting this error:

Error: File: Project_5_2_E.m Line: 9 Column: 1 This statement is not inside any function. (It follows the END that terminates the definition of the function "EventFcn".)

If I put end at the end of the code I get this error then: Not enough input arguments.

Error in Project_5_2_E (line 4) value = [(x(1) - 1); (x(2) - 1)];

[Star Strider]

My pleasure.

Those are all the times on the [0,5] interval that the function returns to approximately [1,1]. The first return, 1.7980, would be my choice if I had to choose one particular time.

[Star Strider]

About the error — you have to troubleshoot that. I have no idea how you set up your code, only that the code I posted ran without error for me.

I would not put the ‘value’ assignment line outside my ‘EventFcn’ function, since it has no context outside it and will throw errors. My code as I posted it gives you everything you need.

[Bob]

Ok that sounds good. I really only need the answer though so since it is working on your end, could you just tell me what you get for (0.5,0.5). I am not sure why I can not get it to work for me. Thanks again!

[Star Strider]

My pleasure.

Please keep troubleshooting it. You need to understand how the code works, and how to make it work. That’s the whole point of my helping you.

For [0.5, 0.5], I get:

ReturnTimes =

     253.6151e-003
       2.1916e+000
       2.8704e+000
       4.8043e+000

Those equalities aren’t quite as impressive as with the earlier initial conditions, but looking at the raw output as well, they’re close enough. I would choose either 0.2536 or 2.1916 here.

[Star Strider]

ERROR!

I just discovered that I forgot to update ‘EventFcn’ for different initial conditions. (The code for [1,1] is correct, as are the results.) This updated code allows for them to be passed to ‘EventFcn’, so it will now work for all initial conditions to detect the return. I defined the initial condition vector separately so it is now created and then passed as an additional parameter to ‘EventFcn’, and directly to your ode45 call:

      function [value,isterminal,direction] = EventFcn(t,x,ic)
      value = [(x(1) - ic(1)); (x(2) - ic(2))];
      isterminal = [0; 0];
      direction = [0; 0];
      end

gh = @(t,x) [4*x(1) - 3*x(1).*x(2); -2*x(2) + x(1).*x(2)];
tspan = linspace(0, 5, 1000);
c = 0.5;
ic = [c; c];
opts = odeset('Events', @(t,x)EventFcn(t,x,ic));
[t,x,EventTime,EventXVal,ie] = ode45(gh, tspan, ic, opts);
EventTime                                                           % Show Values (Delete)
EventXVal                                                           % Show Values (Delete)

x_1_1 = abs(diff(EventXVal,1,2)) < 1.1;                             % Logical Vector Of Approximately Equal ‘x’ Values
EqualTimes = EventTime(x_1_1)                                       % Corresponding Event Times
UniqueTimes = diff([0; EqualTimes]) >= 0.2                          % Use ‘diff’ To Select Unique Times
ReturnTimes = EqualTimes(UniqueTimes)                               % Desired Result

There is only one return time with [0.5, 0.5]:

ReturnTimes =

       2.6050

I apologise for the error. It’s fixed now, and my code is now robust to all initial conditions.