使用 Sympy 的 ODE 积分方法的解析表达式

Question

我正在尝试使用 Sympy 计算常微分方程 (ODE) 数值解误差的解析表达式，例如 Euler 和 Runge-Kutta 方法。

我们有 ODE：$\dot x = f(x)$，初始值为 $x(t_0) = x_0$，我们想找出近似解 $x(t+\Delta t) pprox x_0 的误差+ \Delta t f(x_0) + \Delta t^2 f'(x_0)f(x_0)/2 +...$.

程序是将泰勒级数中的解 $x(t)$ 扩展到比所分析的数值方法高两个数量级，然后将此扩展与数值解的扩展评估进行比较。对于 Euler 方法，这非常简单（容易），因为两个展开都不需要太多代数：

欧拉法：x1 = x0 + \Delta t f(x_0)

泰勒级数，直到$O(\Delta t^2)$：$x(t+\Delta t) = x_0 + \Delta t f(x_0) + \Delta t^2 f'(x_0)f(x_0)/ 2 + O(\Delta t^3) $

错误：$|E| = | \Delta t^2 f'(x_0)f(x_0)/2 | $.

上面的技巧是在$x(吨）$。然而，为了 2 或更高，事情变得更加复杂，所以我想使用 Sympy 来帮助开发表达式。

这是我的尝试：

from sympy import *
init_printing()

t, t0, dt, x, x0, x1, x_dot = symbols("t t_0 {\\Delta}t x x_0 x_1 \\dot{x}") 
xt, f  = symbols("x f", cls=Function) #x(t) e f(x(t))
Erro =symbols("E")
fp = [f] # A Python list to contain the symbols of df/dx, up to order 9
string_fp = "f"
for i in range(1,10):
    string_fp +="'"
    fp.append(symbols(string_fp, cls=Function))

def my_taylor(f, x, x0, order):
    Taylor_expansion = f(x0)
    for i in range(1, order+1):
        Taylor_expansion += (1/factorial(i))*((x-x0)**i)* diff(f(x), x, i)
        #Taylor_expansion += (1/factorial(i))*((x-x0)**i)* diff(f(x), x, i).subs(x,x0)
        #Taylor_expansion += (1/factorial(i))*((x-x0)**i)*fp[i](x0)
    Taylor_expansion += Order((x-x0)**(order+1),(x,x0))
    return Taylor_expansion

#teste = my_taylor(f, x, x0, 2)

def my_other_taylor(xt, t, t0, f, order):
    Taylor_expansion = xt(t0)
    Taylor_expansion += (t-t0)*f(x0)
    for i in range(2, order+1):
        #Taylor_expansion += (1/factorial(i))*((t-t0)**i)* diff(f(x), x, i-1).subs(x,x0)
        #Taylor_expansion += (1/factorial(i))*((t-t0)**i)*fp[i](x0)
        Taylor_expansion += (1/factorial(i))*((t-t0)**i)* diff(f(xt(t)), t, i-1).subs(xt(t),x0)
    #Taylor_expansion += Order((factor(t-t0)**(order+1)))
    Taylor_expansion += Order((t-t0)**(order+1),(t,t0))
    return Taylor_expansion
#teste2 = my_other_taylor(xt, t, t0, f, 5)

#Taylor_f  = f(x).series(x, x0=x0, n=3)
#Taylor_x  = xt(t).series(t, x0=t0, n=3)
#Taylor_x_2 = Taylor_x.subs(diff(xt(t),t), f(x))  ## Esta substituição não está funcionando
Taylor_f = my_taylor(f, x, x0, 2)
Taylor_x = my_other_taylor(xt, t, t0, f, 4)
Taylor_x_dt = Taylor_x.subs(t, dt+t0)

#Taylor_x_2 = Taylor_x.subs(t-t0, dt)
pprint(Taylor_x_dt, use_unicode=True)

### Calcular x1 usando método de Euler
x1 = x0 + dt*f(x0)
Erro = abs(x1-Taylor_x_dt.removeO().subs(xt(t0),x0))

我没有得到正确的结果，我知道为什么，但我不知道如何修复泰勒展开。 $f(x)$ 的导数应该在 $x_0$ 处“求值”，但我无法在 sympy 的文档和教程中找到如何对导数进行符号求值。如果我做

.subs(x, x0)

，它会替换微分中使用的变量，比如$df(x)/dx$变成$df(x_0)/dx_0$，而不是$df(x_0)/dx$。此外，我找不到一种方法来应用链式规则并将 $x(t)$ 的时间导数替换为 $f$ w.r.t $x$ 的导数。