python - 求解非线性微分方程 Sympy

Question

如果 v(t) 未平方，此代码仅适用于求解微分方程 v_equation。当我对它求平方时，它返回错误 PolynomialDivisionFailed。是否有另一种方法可以使用 Sympy 执行此操作，或者我是否应该找到一个不同的 python 包来执行这些类型的计算。

from sympy import *
from matplotlib import pyplot as plt
import numpy as np

m = float(raw_input('Mass:\n> '))
g = 9.8
k = float(raw_input('Drag Coefficient:\n> '))
f1 = g * m
t = Symbol('t')
v = Function('v')
v_equation = dsolve(f1 - k * (v(t) ** 2) - m * Derivative(v(t)), 0)
C1 = Symbol('C1')
C1_ic = solve(v_equation.rhs.subs({t:0}),C1)[0]
v_equation = v_equation.subs({C1:C1_ic})

func = lambdify(t, v_equation.rhs,'numpy')

Answer 1

根据我使用符号数学包的经验，我不建议使用浮点常量执行(符号)计算。最好用符号常数定义方程，尽可能进行计算，然后用数值代入。

通过这种方法，Sympy 可以为这个 D.E. 提供解决方案。

首先，定义符号常量。为了帮助计算，请注意我们可以提供有关这些常量的附加信息(例如，实数、正数等)

import sympy as sp

t = sp.symbols('t', real = True)
g, k, m = sp.symbols('g, k, m', real = True, positive = True)
v = sp.Function('v')

DE的符号解可以得到如下

f1 = g * m
eq = f1 - k * (v(t) ** 2) - m * sp.Derivative(v(t))
sol = sp.dsolve(eq,v(t)).simplify()

解 sol 将是 k、m、g 和常数 的函数>C1。通常，初始条件对应两个复数 C1 值。但是，C1 的两个值在 sol 中替换时会产生相同的(实值)解。

请注意，如果您不需要符号解，则可以使用数值 ODE 求解器，例如 Scipy 的 odeint。代码如下(对于初始条件 0):

from scipy.integrate import odeint

def fun(v, t, m, k, g):
    return (g*m - k*v**2)/m

tn = np.linspace(0, 10, 101)
soln = odeint(fun, 0, tn, args=(1000, 0.2, 9.8))

soln 是样本数组 v(t) 对应于 tn 元素