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?
采纳的回答
Star Strider
2016-8-16
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.
7 个评论
Bob
2016-8-16
Star Strider
2016-8-16
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
2016-8-16
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.
Star Strider
2016-8-16
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
2016-8-16
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
