1

我正在尝试y'' + (epsilon-x^2)y = 0使用 odeint 数值求解方程。我知道解决方案（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

4

1 回答 1

0

看起来您的方程式没有正确解释。你有一个微分方程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

于 2013-12-07T20:43:47.193 回答