python - library - cómo resolver diff. eq. usando scipy.integrate.odeint?

solve differential equation odeint

Question

Quiero resolver estas ecuaciones diferenciales con las condiciones iniciales dadas:

(3x-1)y''''-(3x+2)y''+(6x-8)y=0, y(0)=2, y''(0)=3

los ans deben ser

y=2*exp(2*x)-x*exp(-x)

aquí está mi código:

def g(y,x):
    y0 = y[0]
    y1 = y[1]
    y2 = (6*x-8)*y0/(3*x-1)+(3*x+2)*y1/(3*x-1)
    return [y1,y2]

init = [2.0, 3.0]
x=np.linspace(-2,2,100)
sol=spi.odeint(g,init,x)
plt.plot(x,sol[:,0])
plt.show()

pero lo que obtengo es diferente de la respuesta. ¿qué he hecho mal?

Answer 1

Hay varias cosas mal aquí. En primer lugar, su ecuación es aparentemente

(3x-1) y '''' - (3x + 2) y ''- (6x-8) y = 0; y (0) = 2, y ''(0) = 3

(observe el signo del término en y). Para esta ecuación, su solución analítica y la definición de y2 son correctas.

En segundo lugar, como dice @Warren Weckesser, debe pasar 2 parámetros como y a g : y[0] (y), y[1] (y '') y devolver sus derivadas, y'' e y ''''.

En tercer lugar, sus condiciones iniciales se dan para x = 0, pero su cuadrícula x para integrar comienza en -2. De los documentos para odeint , este parámetro, t en su descripción de la firma de llamada:

odeint(func, y0, t, args=(),...) :

t: array Una secuencia de puntos de tiempo para los cuales resolver para y. El punto de valor inicial debe ser el primer elemento de esta secuencia.

Por lo tanto, debe integrar desde 0 o proporcionar condiciones iniciales comenzando en -2.

Finalmente, su rango de integración cubre una singularidad en x = 1/3. odeint puede tener un mal momento aquí (pero aparentemente no).

Aquí hay un enfoque que parece funcionar:

import numpy as np
import scipy as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def g(y, x):
    y0 = y[0]
    y1 = y[1]
    y2 = ((3*x+2)*y1 + (6*x-8)*y0)/(3*x-1)
    return y1, y2

# Initial conditions on y, y'' at x=0
init = 2.0, 3.0
# First integrate from 0 to 2
x = np.linspace(0,2,100)
sol=odeint(g, init, x)
# Then integrate from 0 to -2
plt.plot(x, sol[:,0], color=''b'')
x = np.linspace(0,-2,100)
sol=odeint(g, init, x)
plt.plot(x, sol[:,0], color=''b'')

# The analytical answer in red dots
exact_x = np.linspace(-2,2,10)
exact_y = 2*np.exp(2*exact_x)-exact_x*np.exp(-exact_x)
plt.plot(exact_x,exact_y, ''o'', color=''r'', label=''exact'')
plt.legend()

plt.show()