¿Software para la solución numérica de un sistema EDO no lineal?

Question

Me han dado un sistema no lineal de ecuaciones diferenciales ordinarias que ha surgido de un colega de investigación de ingeniería:

$$\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}$$

(He valores para $\lambda$, $k$ y $n$, y también algunas de las condiciones iniciales).

Básicamente todo lo que quiero es un (open source) que el sistema informático de aventura, y a ver qué pasa.

Hasta ahora he tratado de Octava y de su "lsode" comando (no es bueno; le dio errores), Python con "odeint" de sympy.integrar (le dio una solución); tengo que probar un par de los demás.

Yo estoy lejos de ser un experto (o incluso vagamente competente) en no-lineal de los sistemas de educación a distancia, y estoy esperando para el consejo como para que el sistema se puede utilizar con confianza para generar confianza soluciones numéricas.

Answer 1

Me gustaría utilizar el Scilab para este trabajo. Hay varias soluciones ODE implementado en Scilab; en el código que se indican a continuación me dejó la opción de solver para el programa, pero puede comparar los resultados de diferentes solucionadores de problemas mediante la especificación de ellas en el comando.

Para valores numéricos solía $\lambda=k=1$$n=3$. El punto de partida es el vector columna $(1,1,4)^T$. Tenga en cuenta que Scilab enumera vector de entradas de inicio con el $1$, no con $0$ como en tu post.

Código:

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

Salida gráfica:

Los resultados numéricos se encuentra en la matriz x.

Answer 2

Imitando user110985 s código en la octava, primero debe definir una función en el archivo:

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

A continuación, un archivo de script de la siguiente manera:

`cerca de todos; clear all; clc;

t = (0:0.05:5)';

a = lsode(@xdot_se1,[1,1,4]',t');

figura(1);

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

plot(t,a(:,2),'g','linewidth',2);

plot(t,a(:,3),'b','linewidth',2);`

Tenemos la similar de la trama: