1

$$x' = \lambda - \rho x - \beta xz \\ y' = \beta xz - \delta y \\ z' = py - cz$$

$x_0=43100, y_0 = 0, z_0 = 0.0033$

$\lambda = 388, \rho = 0.009, \beta =3.61e-8, \delta = 0.18, p = 50000, c = 23. $

my code is as follows:

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]);

Issue: The code keeps running without any result or termination. So I have tried using ode23s and ode15s but I get a number of warning messages like

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 incorrect? Please help!

Michael Medvinsky
  • 5,728
ChrisK
  • 13
  • I presume this is Matlab syntax? You didn't specify in your question... – Xoque55 Dec 05 '15 at 15:52
  • Yes, my apologies! – ChrisK Dec 05 '15 at 15:54
  • 1
    If it's code then don't use invalid special symbols. Show your actual code. If you want to show an equation/formula then use TeX to format it. Errors like you show can be due to many things, but one common one is typos in your ODE function. Check that you don't have any sign errors, etc. Also, the scaling of your problem may also be an issue, i.e., parameters like 0.0000000361 and 50000 used together. – horchler Dec 05 '15 at 19:39

1 Answers1

1

After trying to re-scale your system I've just found a little typo in your code in a last component of f(t,x). Now everything is working, but if you want you can use my version of your code, it looks different, but after all all the changes just minor.

lambda = 388;
rho = 0.009;
beta=3.61e-8;
delta=0.18;
p=50000;
c=23;

f = @(t,x) [lambda-rho*x(1)-beta*x(1)*x(3); beta*x(1)*x(3) - delta*x(2); p*x(2) - c*x(3)];
tspan = [0,1];
options = odeset('RelTol',1e-1,'AbsTol',[1e-4 1e-4 1e-4]);
[t,xa] = ode45(f,tspan,[43100 0 0.0033],options);
plot(t,xa(:,1),'-'); figure
plot(t,xa(:,2),'-.');figure
plot(t,xa(:,3),'.');
Michael Medvinsky
  • 5,728