scipy.odeint strange behavior

Question

Here is my code to solve differential equation dy / dt = 2 / sqrt(pi) * exp(-x * x) to plot 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()

And here is a plot:

plot

Am I using odeint in a wrong way or it's a bug?

unutbu · Accepted Answer · 2012-11-04 16:18:49Z

Notice that if you change x to x = np.linspace(-5.0, 5.0, 10000), then your code works. Therefore, I suspect the problem has something to do with exp(-x*x) being too small when x is very small or very large. [Total speculation: Perhaps the odeint (lsoda) algorithm adapts its stepsize based on values sampled around x = -10 and increases the stepsize in such a way that values around x = 0 are missed?]

The code can be fixed by using the tcrit parameter, which tells odeint to pay special attention around certain critical points.

So, by setting

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

we tell odeint to sample more carefully around 0.

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()

Collectives™ on Stack Overflow

scipy.odeint strange behavior

1 Answer 1

Comments

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Related