# Software for numerical solution of a non-linear ODE system?

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.

• Did you try maxima, sage, or maybe some on here en.wikipedia.org/wiki/List_of_computer_algebra_systems? – Amzoti Nov 22 '13 at 2:09
• 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. – 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.

Code:

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)]
endfunction
x0=[1;1;4]
t0=0
t=0:0.05:5
x = ode(x0,t0,t,f);
plot(t,x(1,:),'r')
plot(t,x(2,:),'g')
plot(t,x(3,:),'b')

Graphical output:

The numerical output is in the matrix x.

• 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. – 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;
endfunction

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

figure(1);

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

We get the similar plot: