这是我解微分方程dy/dt=2/sqrt(pi)*exp(-x*x)绘制erf(x)的代码。在
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
import math
def euler(df, f0, x):
h = x[1] - x[0]
y = [f0]
for i in xrange(len(x) - 1):
y.append(y[i] + h * df(y[i], x[i]))
return y
def i(df, f0, x):
h = x[1] - x[0]
y = [f0]
y.append(y[0] + h * df(y[0], x[0]))
for i in xrange(1, len(x) - 1):
fn = df(y[i], x[i])
fn1 = df(y[i - 1], x[i - 1])
y.append(y[i] + (3 * fn - fn1) * h / 2)
return y
if __name__ == "__main__":
df = lambda y, x: 2.0 / math.sqrt(math.pi) * math.exp(-x * x)
f0 = 0.0
x = np.linspace(-10.0, 10.0, 10000)
y1 = euler(df, f0, x)
y2 = i(df, f0, x)
y3 = odeint(df, f0, x)
plt.plot(x, y1, x, y2, x, y3)
plt.legend(["euler", "modified", "odeint"], loc='best')
plt.grid(True)
plt.show()
下面是一个情节:
我是用错了odeint还是错误的错误?在
请注意,如果您将}非常小或非常大。[总体推测:也许odeint(lsoda)算法根据在
x
更改为x = np.linspace(-5.0, 5.0, 10000)
,那么您的代码就可以工作了。因此,我怀疑这个问题与exp(-x*x)
太小有关,而{x = -10
周围采样的值来调整其步长,并以这样的方式增加步长,从而忽略x = 0
周围的值?]在代码可以通过使用
tcrit
参数来修复,该参数告诉odeint
要特别注意某些临界点。在所以,通过设置
我们告诉
^{pr2}$odeint
在0附近更仔细地采样。在相关问题 更多 >
编程相关推荐