SymPy的dsolve的返回值中的r()函数是啥意思？

Posted 2023-02-18

【中文标题】SymPy的dsolve的返回值中的r()函数是啥意思？

【英文标题】：What does r() function mean in the return value of SymPy's dsolve?SymPy的dsolve的返回值中的r()函数是什么意思？

【发布时间】：2020-11-08 20:43:01

【问题描述】：

我想评估phi(+oo)的价值 其中phi(xi)是ODE的解

Eq(Derivative(phi(xi), (xi, 2)), (-K + xi**2)*phi(xi))

而K 是一个已知的实变量。 通过dsolve，我得到了解决方案：

Eq(phi(xi), -K*xi**5*r(3)/20 + C2*(K**2*xi**4/24 - K*xi**2/2 + xi**4/12 + 1) + C1*xi*(xi**4/20 + 1) + O(xi**6))

在右侧的第一项中具有未知函数r()。 这是我的代码：

import numpy as np
import matplotlib.pyplot as plt
import sympy
from sympy import I, pi, oo

sympy.init_printing()

def apply_ics(sol, ics, x, known_params):
    """
    Apply the initial conditions (ics), given as a dictionary on
    the form ics = y(0): y0, y(x).diff(x).subs(x, 0): yp0, ...,
    to the solution of the ODE with independent variable x.
    The undetermined integration constants C1, C2, ... are extracted
    from the free symbols of the ODE solution, excluding symbols in
    the known_params list.
    """

    free_params = sol.free_symbols - set(known_params)
    eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics)
            for n in range(len(ics))]
    sol_params = sympy.solve(eqs, free_params)

    return sol.subs(sol_params)

K = sympy.Symbol('K', positive = True)
xi = sympy.Symbol('xi',real = True)
phi = sympy.Function('phi')
ode = sympy.Eq( phi(xi).diff(xi, 2), (xi**2-K)*phi(xi))

ode_sol = sympy.dsolve(ode)
ics =  phi(0):1, phi(xi).diff(xi).subs(xi,0): 0
phi_xi_sol = apply_ics(ode_sol, ics, xi, [K])

其中ode_sol 是解，phi_xi_sol 是应用初始条件后的解。 由于 r() 在 NumPy 中未定义，我无法评估结果

for g in [0.9, 0.95, 1, 1.05, 1.2]:
    phi_xi = sympy.lambdify(xi, phi_xi_sol.rhs.subs(K:g), 'numpy')

有谁知道r()这个函数是什么意思，我应该怎么处理？

【问题讨论】：

有了足够多的symbols 命令，我能够重新创建您的解决方案。要进一步进行，您/我们需要研究 dsolve 文档 - 可能会添加一个或多个参数。

谢谢你，@hpaulj。通常，在求解二阶 ODE 时，只需定义两个初始/边界即可定义未知参数。但是在我应用它们之后（参见添加的代码），我得到的“phi_xi_sol”中仍然存在函数r（）。

您使用的是最新的 sympy 版本吗？第二个公式的解在数学上看起来是错误的，第一项应该属于系数为C1的奇次幂基解，看来解的过程提前结束了。

用sympy 1.5.1检查，结果是一样的。

还是和sympy 1.7.1一样

【参考方案1】：

从结果的形式中可以看出，求解器回退到幂级数解（而不是像 WolframAlpha 那样根据抛物柱面函数搜索解）。

所以让我们设置phi(xi)=sum a[k]*xi^k 导致系数方程（使用a[k]=0 代表k<0）

(k+2)(k+1)a[k+2] = -K*a[k] + a[k-2]

a[0] = C2  
a[1] = C1  
a[2] = -K/2*C2
a[3] = -K/6*C1
a[4] = (K^2/2 + 1)/12*C2
a[5] = (K^2/6 + 1)/20*C1

插入幂级数解应该是

C2*(1-K/2*xi**2+(K**2/24+1/12)*xi**4) + C1*xi*(1-K/6*xi**2+(K/120+1/20)*xi**4) + O(xi**6)

与 sympy 解决方案相比，所有包含 C1 和 K 的术语都丢失了，尤其是缺少的 3 级术语无法解释。求解过程似乎过早结束，或者某些方程变换未正确反转。

请注意，sympy 中的 ODE 求解器例程是实验性的和初级的。此外，幂级数解决方案仅提供 xi 的小值的有效信息，无法推导出 +oo 的限制的任何确切值。

【参考方案2】：

sol_params 是一个包含单个字典的列表。传递该字典而不是列表给出了解决方案phi_xi_sol，而没有r(3)：

Eq(rho(s), (-K*s**2/2 + s**2*xi**2/2 + 1)*(-6*rho(s) + 6*C2*s - C2*K*s**3 + O(s**5))/
            (3*(K*s**2 - 2)) + C2*(-K*s**3/6 + s**3*xi**2/6 + s) + O(s**5))