Как решить дифференциальное уравнение, используя встроенную функцию Python odeint?

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

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

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


person Physicist    schedule 07.01.2015    source источник


Ответы (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), y[1] (y') и вернуть их производные, y' и y''.

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

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

t : массив Последовательность моментов времени, для которых необходимо найти y. Точка начального значения должна быть первым элементом этой последовательности.

Таким образом, вы должны интегрировать, начиная с 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()

введите здесь описание изображения

person xnx    schedule 07.01.2015
comment
Для дифференциального уравнения второго порядка init должна иметь длину 2, а не 3 (и g должна возвращать массив длины 2). - person Warren Weckesser; 07.01.2015
comment
Вы правы: я запутался. Я отредактировал, чтобы исправить это. - person xnx; 07.01.2015