Python：odeint 求解具有可变系数的 ODE (QHO)

Question

我正在尝试使用 odeint 以数值方式求解方程 y'' + (epsilon-x^2)y = 0。我知道解决方案（QHO 的波函数），但 odeint 的输出与它没有明显关系。我可以很好地解决具有常数系数的 ODE，但是一旦我转向可变系数，我就无法解决我尝试过的任何一个。这是我的代码：

#!/usr/bin/python2
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi

x = np.linspace(-5,5,1e4)

n = 0
epsilon = 2*n+1 

def D(Y,x):
    return np.array([Y[1], (epsilon-x**2)*Y[0]])

Y0 = [0,1]

Y = spi.odeint(D,Y0,x)
# Y is an array with the first column being y(x) and the second y'(x) for all x
plt.plot(x,Y[:,0],label='num')
#plt.plot(x,Y[:,1],label='numderiv')

plt.legend()
plt.show()

还有剧情： [没有足够的代表：]https://drive.google.com/file/d/0B6840LH2NhNpdUVucUxzUGFpZUk/edit?usp=sharing

在这里查看解图：http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hosc5.html

Answer 1

您的方程式似乎没有正确解释。您有一个微分方程y'' + (epsilon-x^2)y = 0，但您忘记了向量形式中的减号。特别是应该是

y[0]' = y[1]
y[1]' = -(epsilon-x^2)y[0]

所以（在 ​​epsilon 项前添加减号

def D(Y,x):
    return np.array([Y[1], -(epsilon-x**2)*Y[0]])

事实上，你所拥有的情节与 DE y'' + (epsilon-x^2)y = 0 是一致的。看看：Wolphram Alpha