Как решить дифференциальное уравнение с помощью встроенной функции Python odeint?

Question

Я хочу, чтобы решить эти дифференциальные уравнения с заданными начальными условиями:

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

АНС должен быть

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

вот мой код:

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

но то, что я получаю отличается от ответа. что я сделал не так?

Answer 1

Здесь несколько вещей не так. Во-первых, ваше уравнение, по-видимому,

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

(обратите внимание на знак термина в y). Для этого уравнения ваше аналитическое решение и определение y2 верны.

Во-вторых, как говорит @Warren Weckesser, вы должны передать 2 параметра, как y к g: y[0] (у), y[1] (у ') и возвращать их производные, у' и у '.

В-третьих, ваши начальные условия задаются для x = 0, но ваша x-сетка для интеграции начинается с -2. Из документов для odeint, этот параметр, t в их описании подписи вызова:

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

т: массив последовательность моментов времени, для которых решить для у. Начальная точка должна быть первым элементом этой последовательности.

Таким образом, вы должны интегрировать начало в 0 или обеспечить начальные условия, начиная с -2.

Наконец, ваш диапазон интеграции охватывает особенность при x = 1/3. odeint может быть плохой день (но, видимо, нет).

Вот один подход, который, кажется, работает:

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