我在double_scalars中遇到溢出，在向后欧拉中遇到double_scalars错误时遇到无效值

Question

我正在尝试为不同的dts绘制向后的欧拉方法，代码如下所示：

import numpy as np
import matplotlib.pyplot as plt

def newton(x_old, t_new, dt, precision=1e-12):
""" Takes values x_old and t_new, and finds the root of the
function f(x_new), returning x_new. """
# initial guess:
x_new = x_old
f = func(x_new, x_old, t_new, dt)
dfdx = dfuncdx(t_new, dt)
while abs(f / dfdx) > precision:
x_new = x_new - f / dfdx
f = func(x_new, x_old, t_new, dt)
dfdx = dfuncdx(t_new, dt)
return x_new

def func(x_new, x_old, t_new, dt):
""" The function f(x) we want the root of."""
return x_new - x_old - dt * (-2*x_new*t_new)

def dfuncdx(t_new, dt):
""" The derivative of f(x) with respect to x_new."""
return 1 + dt * (-2*t_new)

dts = np.array([0.2, 0.15, 0.1, 0.01])
plt.figure()
for i, dt in enumerate(dts):
t_max = 10  
N = int(t_max / dt)
t = np.linspace(0, 10, N + 1)
x = np.zeros(N + 1)
t[0] = 0
x[0] = 1

for n in range(N):
x[n + 1] = newton(x[n], t[n + 1], dt)

plt.plot(t, x)
t = np.linspace(0, 10, 10000)
plt.plot(t, np.exp(-t**2))
plt.xlabel(r'$t$')
plt.ylabel(r'$x(t)$')
plt.legend([r'$Delta t=$ %0.2f' % dts[0], r'$Delta t=$ %0.2f' % dts[1], r'$Delta t=$ %0.2f' % dts[2],
r'$Delta t=$ %0.2f' % dts[3], r'Exact solution'])
plt.grid();
plt.show()

代码完成，我也得到了看起来正确的解决方案，但是图表在它们应该被切断之前就被切断了，我不明白为什么。我还遇到double_scalars溢出double_scalars错误中遇到的无效值。

Answer 1

我很愚蠢，在导数中犯了一个错误：/这现在有效

import numpy as np
import matplotlib.pyplot as plt

def newton(x_old, t_new, dt, precision=1e-6):
""" Takes values x_old and t_new, and finds the root of the
function f(x_new), returning x_new. """
# initial guess:
x_new = x_old
f = func(x_new, x_old, t_new, dt)
dfdx = dfuncdx(t_new, dt)
while abs(f / dfdx) > precision:
x_new = x_new - f / dfdx
f = func(x_new, x_old, t_new, dt)
dfdx = dfuncdx(t_new, dt)
return x_new

def func(x_new, x_old, t_new, dt):
""" The function f(x) we want the root of."""
return x_new - x_old - dt * (-2*x_new*t_new)

def dfuncdx(t_new, dt):
""" The derivative of f(x) with respect to x_new."""
return 1 + dt * (2*t_new)
dts = np.array([0.3, 0.2, 0.1, 0.0001])

plt.figure()
for i, dt in enumerate(dts):
t_max = 10  # maximum time
N = int(t_max / dt)
t = np.linspace(0, 10, N + 1)
x = np.zeros(N + 1)
t[0] = 0
x[0] = 1
t_old = t[0]
x_old = x[0]
for n in range(N):
x[n + 1] = newton(x[n], t[n + 1], dt)
t[n + 1] = t[n] + dt
# Plot the solution
plt.plot(t, x)
t = np.linspace(0, 10, 10000)
plt.plot(t, np.exp(-t**2))
plt.xlabel(r'$t$')
plt.ylabel(r'$x(t)$')
plt.legend([r'$Delta t=$ %0.2f' % dts[0], r'$Delta t=$ %0.2f' % dts[1], r'$Delta t=$ %0.2f' % dts[2],
r'$Delta t=$ %0.2f' % dts[3], r'Exact solution'])
plt.grid();
plt.show()