python - 我在反向欧拉中遇到双标量溢出和双标量错误中遇到的无效值

Question

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

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 错误中遇到的无效值。

标签： pythonnumpy

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