The following code doesn't seem to work - it keeps running with no result (RK4 method)

4 views (last 30 days)
There are 3 differential equations as defined in 'f' below separated by semi-colon. The time span is from 0 to 84 by 7.
f = @(t,x) [388-0.009*x(1)-0.0000000361*x(1)*x(3); 0.0000000361*x(1)*x(3) - 0.18*x(2); 50000*x(3) - 23*x(3)];
tspan = 0:7:84;
[t,xa] = ode45(f,tspan,[43100 0 0.0033]);
  3 Comments
KC
KC on 5 Dec 2015
Edited: KC on 5 Dec 2015
x' = λ - ρx - βxz;
y' = βxz - δy;
z' = py - cz;
x0=43100; y0 = 0, z0 = 0.0033, λ = 388, ρ = 0.009 δ = 0.18, p = 50000, c = 23, β=3.61e-8

Sign in to comment.

Accepted Answer

Marc
Marc on 6 Dec 2015
You figure out your mistake yet??
Your code is just "fine"... But your function as you have it above is "wrong".
Take a look at your third equation. Change that first x(3) to a x(2) and you will be in business.
No issues here using ODE45 or ODE15s. Also, you don't need to set tspan as such... [0:7:84] You can simply use [0 84] and you get a smoother answer. This took well under a second for both 45 and 15s....
Interestingly, you should try both and compare. ODE45 gives you a solution in ~0.1-0.2 seconds (depending on your machine) and the solution is ~2400 points. 15s gives you an answer in under 0.1s and the solution is 163 points. Just goes to show you how much more efficient ODE15s can be with stiff problems.

More Answers (1)

John D'Errico
John D'Errico on 5 Dec 2015
Edited: John D'Errico on 5 Dec 2015
This is surely what would be called a stiff problem. It will cause ODE45 (or any non-stiff solver like it) to run very poorly. The first clue to the problem being a stiff one is the fact that ODE45 DID work poorly. A better clue is the wide disparity of magnitudes of the coefficients in that ODE. That tends to create stiff systems.
So instead, try one of the stiff solvers. All of the stiff solvers have names that end in s, thus ode23s, ode15s.
  1 Comment
KC
KC on 5 Dec 2015
Edited: KC on 5 Dec 2015
I just tried using both ode23s and ode15s, but I get many warning messages like these: > In ode23s (line 379) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.484584e-246. Warning: Failure at t=7.014056e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.491893e-14) at time t. I am wondering if my code is correct?

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!