0

これは、微分方程式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を間違った方法で使用していますか、それともバグですか？

4

1 に答える 1

2

xに変更するx = np.linspace(-5.0, 5.0, 10000)と、コードが機能することに注意してください。したがって、問題が非常に小さいか非常に大きいexp(-x*x)場合、問題は小さすぎることに関係しているのではないかと思います。x[完全な推測：おそらく、odeint（lsoda）アルゴリズムは、サンプリングされた値に基づいてステップサイズを適応させ、周りx = -10の値が失われるようにステップサイズを増やしx = 0ますか？]

このコードは、特定の重要なポイントに特別な注意を払うようtcritに指示するパラメーターを使用して修正できます。odeint

だから、設定することによって

y3 = integrate.odeint(df, f0, x, tcrit = [0])

odeint0付近でより注意深くサンプリングするように指示します。

import matplotlib.pyplot as plt
import scipy.integrate as integrate
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

def df(y, x):
   return 2.0 / np.sqrt(np.pi) * np.exp(-x * x)

if __name__ == "__main__":
    f0 = 0.0
    x = np.linspace(-10.0, 10.0, 10000)

    y1 = euler(df, f0, x)
    y2 = i(df, f0, x)
    y3 = integrate.odeint(df, f0, x, tcrit = [0])

    plt.plot(x, y1)
    plt.plot(x, y2)
    plt.plot(x, y3)
    plt.legend(["euler", "modified", "odeint"], loc='best')
    plt.grid(True)
    plt.show()
于 2012-11-04T16:18:49.050 に答える