这里有几个问题 . 首先,你的等式显然是
(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 ') and return their derivatives, y'和y'' .
第三,你的初始条件是x = 0,但你要整合的x网格从-2开始 . 从 odeint 的文档中,此参数 t 在其调用签名描述中:
odeint(func, y0, t, args=(),...) :
t:array要求解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()