我正试图使用Sympy来计算常微分方程(ODEs)数值解误差的解析表达式，例如Euler和Runge-Kutta方法。
我们有初始值为$x(t_0)=x_0$的ODE:$\dot x=f(x)$，并且我们想要找到近似解$x(t+\Delta t)\approxy x_0+\Delta tf(x_0)+\Delt t^2 f'(x_0。
该程序是将泰勒级数中的解$x(t)$扩展到比分析中的数值方法高两个阶，然后将此扩展与数值解的扩展评估进行比较。对于欧拉方法来说，这是非常简单的，因为这两种展开都不需要太多代数：
欧拉法：x1=x0+\Δt f(x_0)
泰勒级数，至多$O(\Δt^2)$：$x(t+\Δt)=x_0+\Δt f(x_0)+\Δt^2 f'(x_0
错误：$|E|=|\Δt^2 f'(x_0)f(x_0。
上面的技巧是在$x(t)$的泰勒展开中使用$dx/dt=f(x)$和$d^2x/dt^2=(df/dx)(dx/dt)$(链式规则)。然而，按照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$。此外，我找不到一种方法来应用链式规则，并用$f$w.r.t$x$的导数代替$x(t)$的时间导数。
您想比较x(t+h)(x'(t)=f(x(t))和stepper(f,t,x(t),h)的解决方案)在h中的Taylor展开式，以获得在stepper中实现的方法的截断误差。对于高阶方法，可以使用基于根树的B序列来表示f的衍生物的组成。
实现这一点的最简单方法是在过程中将f(x+dx)实现为多项式或截断泰勒级数
import sympy as sy
sy.init_printing()
x,h = sy.symbols("x h")
f = sy.symbols("f_{:20}", commutative=False)
def eval_f(y,n):
res = f[n]
dx = y-x
for k in range(n):
res = res/(n-k)*dx+f[n-1-k]
return res
然后通过Picard迭代可以找到精确的解决方案，这是一个缓慢的过程
xe = x+sy.Order(h,h);
for k in range(5):
Dxe = eval_f(xe,k+1)
xe = x+sy.integrate(Dxe,h)
print("nexact ",k,": ",xe.expand().simplify())
其具有输出
exact 0 : x + h*f_{0} + O(h**2)
exact 1 : x + h*f_{0} + h**2*f_{1}*f_{0}/2 + O(h**3)
exact 2 : x + h*f_{0} + h**2*f_{1}*f_{0}/2 + h**3*f_{2}*f_{0}**2/6 + h**3*f_{1}**2*f_{0}/6 + O(h**4)
exact 3 : x + h*f_{0} + h**2*f_{1}*f_{0}/2 + h**3*f_{2}*f_{0}**2/6 + h**3*f_{1}**2*f_{0}/6 + h**4*f_{3}*f_{0}**3/24 + h**4*f_{2}*f_{1}*f_{0}**2/16 + h**4*f_{2}*f_{0}*f_{1}*f_{0}/16 + h**4*f_{1}**3*f_{0}/24 + h**4*f_{1}*f_{2}*f_{0}**2/24 + O(h**5)
exact 4 : x + h*f_{0} + h**2*f_{1}*f_{0}/2 + h**3*f_{2}*f_{0}**2/6 + h**3*f_{1}**2*f_{0}/6 + h**4*f_{3}*f_{0}**3/24 + h**4*f_{2}*f_{1}*f_{0}**2/16 + h**4*f_{2}*f_{0}*f_{1}*f_{0}/16 + h**4*f_{1}**3*f_{0}/24 + h**4*f_{1}*f_{2}*f_{0}**2/24 + h**5*f_{4}*f_{0}**4/120 + h**5*f_{3}*f_{1}*f_{0}**3/60 + h**5*f_{3}*f_{0}**2*f_{1}*f_{0}/60 + h**5*f_{3}*f_{0}*f_{1}*f_{0}**2/60 + h**5*f_{2}**2*f_{0}**3/60 + h**5*f_{2}*(f_{1}*f_{0})**2/40 + h**5*f_{2}*f_{1}**2*f_{0}**2/60 + h**5*f_{2}*f_{0}*f_{2}*f_{0}**2/60 + h**5*f_{2}*f_{0}*f_{1}**2*f_{0}/60 + h**5*f_{1}**4*f_{0}/120 + h**5*f_{1}**2*f_{2}*f_{0}**2/120 + h**5*f_{1}*f_{3}*f_{0}**3/120 + h**5*f_{1}*f_{2}*f_{1}*f_{0}**2/80 + h**5*f_{1}*f_{2}*f_{0}*f_{1}*f_{0}/80 + O(h**6)
现在可以实现常用的数值方法，如二阶Heun方法
k1 = eval_f(x,0)
k2 = eval_f(x+h*k1,3) +sy.Order(h**3,h);
xheun = x+((k1+k2)*h/2).expand()
print("Heun: ", xheun)
print("error: ", xheun-xe)
结果
Heun: x + h*f_{0} + h**2*f_{1}*f_{0}/2 + h**3*f_{2}*f_{0}**2/4 + O(h**4)
error: h**3*f_{2}*f_{0}**2/12 - h**3*f_{1}**2*f_{0}/6 + O(h**4)
或二阶显式中点法
k1 = eval_f(x,0)
k2 = eval_f(x+h*k1,3) +sy.Order(h**3,h);
xheun = x+((k1+k2)*h/2).expand()
print("Heun: ", xheun)
print("error: ", xheun-xe)
print()
k2 = eval_f(x+k1*h/2,3)+sy.Order(h**3,h);
xmid = x+(k2*h).expand()
xmid = xmid.expand()
print("midpoint: ", xmid)
print("error: ", xmid-xe)
结果
midpoint: x + h*f_{0} + h**2*f_{1}*f_{0}/2 + h**3*f_{2}*f_{0}**2/8 + O(h**4)
error: -h**3*f_{2}*f_{0}**2/24 - h**3*f_{1}**2*f_{0}/6 + O(h**4)
最后，标准RK4方法
k2 = eval_f(x+k1*h/2,4).expand() +sy.Order(h**5,h);
k3 = eval_f(x+k2*h/2,4).expand() +sy.Order(h**5,h);
k4 = eval_f(x+k3*h ,4).expand() +sy.Order(h**5,h);
xrk4 = x+((k1+2*k2+2*k3+k4)*h/6).expand()
print("nRK4:")
print("k1 = ",k1)
print("k2 = ",k2)
print("k3 = ",k3)
print("k4 = ",k4)
print("nxRK4 = ",xrk4)
print("error = ",(xrk4-xe))
带输出
k1 = f_{0}
k2 = f_{0} + h*f_{1}*f_{0}/2 + h**2*f_{2}*f_{0}**2/8 + h**3*f_{3}*f_{0}**3/48 + h**4*f_{4}*f_{0}**4/384 + O(h**5)
k3 = f_{0} + h*f_{1}*f_{0}/2 + h**2*f_{2}*f_{0}**2/8 + h**2*f_{1}**2*f_{0}/4 + h**3*f_{3}*f_{0}**3/48 + h**3*f_{2}*f_{1}*f_{0}**2/16 + h**3*f_{2}*f_{0}*f_{1}*f_{0}/16 + h**3*f_{1}*f_{2}*f_{0}**2/16 + h**4*f_{4}*f_{0}**4/384 + h**4*f_{3}*f_{1}*f_{0}**3/96 + h**4*f_{3}*f_{0}**2*f_{1}*f_{0}/96 + h**4*f_{3}*f_{0}*f_{1}*f_{0}**2/96 + h**4*f_{2}**2*f_{0}**3/64 + h**4*f_{2}*f_{1}*f_{0}*f_{1}*f_{0}/32 + h**4*f_{2}*f_{0}*f_{2}*f_{0}**2/64 + h**4*f_{1}*f_{3}*f_{0}**3/96 + O(h**5)
k4 = f_{0} + h*f_{1}*f_{0} + h**2*f_{2}*f_{0}**2/2 + h**2*f_{1}**2*f_{0}/2 + h**3*f_{3}*f_{0}**3/6 + h**3*f_{2}*f_{1}*f_{0}**2/4 + h**3*f_{2}*f_{0}*f_{1}*f_{0}/4 + h**3*f_{1}**3*f_{0}/4 + h**3*f_{1}*f_{2}*f_{0}**2/8 + h**4*f_{4}*f_{0}**4/24 + h**4*f_{3}*f_{1}*f_{0}**3/12 + h**4*f_{3}*f_{0}**2*f_{1}*f_{0}/12 + h**4*f_{3}*f_{0}*f_{1}*f_{0}**2/12 + h**4*f_{2}**2*f_{0}**3/16 + h**4*f_{2}*f_{1}**2*f_{0}**2/8 + h**4*f_{2}*f_{1}*f_{0}*f_{1}*f_{0}/8 + h**4*f_{2}*f_{0}*f_{2}*f_{0}**2/16 + h**4*f_{2}*f_{0}*f_{1}**2*f_{0}/8 + h**4*f_{1}**2*f_{2}*f_{0}**2/16 + h**4*f_{1}*f_{3}*f_{0}**3/48 + h**4*f_{1}*f_{2}*f_{1}*f_{0}**2/16 + h**4*f_{1}*f_{2}*f_{0}*f_{1}*f_{0}/16 + O(h**5)
xRK4 = x + h*f_{0} + h**2*f_{1}*f_{0}/2 + h**3*f_{2}*f_{0}**2/6 + h**3*f_{1}**2*f_{0}/6 + h**4*f_{3}*f_{0}**3/24 + h**4*f_{2}*f_{1}*f_{0}**2/16 + h**4*f_{2}*f_{0}*f_{1}*f_{0}/16 + h**4*f_{1}**3*f_{0}/24 + h**4*f_{1}*f_{2}*f_{0}**2/24 + 5*h**5*f_{4}*f_{0}**4/576 + 5*h**5*f_{3}*f_{1}*f_{0}**3/288 + 5*h**5*f_{3}*f_{0}**2*f_{1}*f_{0}/288 + 5*h**5*f_{3}*f_{0}*f_{1}*f_{0}**2/288 + h**5*f_{2}**2*f_{0}**3/64 + h**5*f_{2}*f_{1}**2*f_{0}**2/48 + h**5*f_{2}*f_{1}*f_{0}*f_{1}*f_{0}/32 + h**5*f_{2}*f_{0}*f_{2}*f_{0}**2/64 + h**5*f_{2}*f_{0}*f_{1}**2*f_{0}/48 + h**5*f_{1}**2*f_{2}*f_{0}**2/96 + h**5*f_{1}*f_{3}*f_{0}**3/144 + h**5*f_{1}*f_{2}*f_{1}*f_{0}**2/96 + h**5*f_{1}*f_{2}*f_{0}*f_{1}*f_{0}/96 + O(h**6)
error = h**5*f_{4}*f_{0}**4/2880 + h**5*f_{3}*f_{1}*f_{0}**3/1440 + h**5*f_{3}*f_{0}**2*f_{1}*f_{0}/1440 + h**5*f_{3}*f_{0}*f_{1}*f_{0}**2/1440 - h**5*f_{2}**2*f_{0}**3/960 - h**5*f_{2}*(f_{1}*f_{0})**2/40 + h**5*f_{2}*f_{1}**2*f_{0}**2/240 + h**5*f_{2}*f_{1}*f_{0}*f_{1}*f_{0}/32 - h**5*f_{2}*f_{0}*f_{2}*f_{0}**2/960 + h**5*f_{2}*f_{0}*f_{1}**2*f_{0}/240 - h**5*f_{1}**4*f_{0}/120 + h**5*f_{1}**2*f_{2}*f_{0}**2/480 - h**5*f_{1}*f_{3}*f_{0}**3/720 - h**5*f_{1}*f_{2}*f_{1}*f_{0}**2/480 - h**5*f_{1}*f_{2}*f_{0}*f_{1}*f_{0}/480 + O(h**6)
有了最后一个主要错误项，人们可以开始看到这种方法的局限性。在将导数CCD_ 9设置为"0"的同时；可交换的"；会结合一些术语，会结合太多术语。将它们设置为非交换项会使一些项分离，这些项实际上与向量值系统函数的导数在其向量自变量中是对称的相同。因此，需要引入严格遵循这些规则的特殊代数类型。