I have been given a nonlinear system of ODEs which has arisen out of a colleague's engineering research:

$$\begin{array}{rcl} \dot{x}_0&=&x_1\\ \dot{x}_1&=&-\frac{\lambda}{(x_2)^n-k^2\lambda}x_0\\ \dot{x}_2&=&\frac{2nx_0x_1\lambda(x_2)^n}{2n(x_2)^{2n-1}-2n\lambda k^2(x_2)^{n-1}-n(n-1)(x_0)^2\lambda(x_2)^{n-1}} \end{array}$$

(I have values for $\lambda$, $k$ and $n$, and also some initial conditions).

Basically all I want is an (open source) computer system to fling them at, and see what happens.

I have so far tried Octave and its "lsode" command (no good; gave errors), Python with "odeint" from sympy.integrate (gave a solution); I need to test out a few others.

I am far from being an expert (or even vaguely competent) at non-linear ODE systems, and I'm hoping for advice as to which system I can use with confidence to generate trustworthy numerical solutions.

  • $\begingroup$ Did you try maxima, sage, or maybe some on here en.wikipedia.org/wiki/List_of_computer_algebra_systems? $\endgroup$ – Amzoti Nov 22 '13 at 2:09
  • $\begingroup$ I've given maxima (its "rk" command) a go; Sage I'm not sure about. Anyway, for numerical work Sage would just call numpy/scipy. $\endgroup$ – Alasdair Nov 22 '13 at 2:28

I would use Scilab for this job. There are several ODE solvers implemented in Scilab; in the code given below I left the choice of solver to the program, but you can compare the results of different solvers by specifying them in the command.

For numerical values I used $\lambda=k=1$ and $n=3$. Starting point is the column vector $(1,1,4)^T$. Note that Scilab enumerates vector entries starting with $1$, not with $0$ as in your post.


function xdot=f(t, x)
    xdot=[x(2); -x(1)/(x(3)^3-1); 6*x(1)*x(2)*x(3)^3/(6*x(3)^5-6*x(3)^2-6*x(1)^2*x(3)^2)]
x = ode(x0,t0,t,f);

Graphical output:


The numerical output is in the matrix x.

  • $\begingroup$ Thanks for that. Scilab just wraps the ODEPACK routines, so I should be able to do exactly the same thing in Octave, which also does. $\endgroup$ – Alasdair Nov 25 '13 at 7:43

Mimicking user110985 s code in octave, first define a function in it's own file:

function lRet = xdot_se1 (x, t)
  lambda = 1;
  k = 1;
  n = 3;
  xdot    = zeros(3,1);
  xdot(1) = x(2);
  xdot(2) = -lambda/(x(3).^(n) - k^2*lambda + eps);
  xdot(3) = 2*n*x(1)*x(2)*lambda*x(3)^n ./ (eps + 2*n*x(3)^(2*n-1) - 2*n*lambda*k^2*x(3)^(n-1) -n*(n-1)*x(1)^2*lambda*x(3)^(n-1));
  lRet = xdot;

Then a script file doing the following:

close all; clear all; clc;

t = (0:0.05:5)';
a = lsode(@xdot_se1,[1,1,4]',t');


plot(t',a(:,1),'r','linewidth',2); hold on;

We get the similar plot: enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.