在含平方函数的微分方程上使用solve_ivp

Question

我在用Python尝试数值解这个方程 (x’(t))^2 + ax^2(t) - b = 0。对于简单的方程，我知道可以用scipy的solve_ivp来解决，只需要像这样定义函数：

def f(t, y):
   return y ** 2

这表示方程 y’(t) = y**2。那我该怎么把第一个方程写成这样呢？我知道这个方程可能可以通过某种替换简化（比如 y=x^2 或其他），但这次我只是想直接（而且是数值上）解决它。

我试着直接求解第一个方程的 x’(t)，这意味着在Python中定义的 f 包含了 math.sqrt()。但是这让我遇到了一个域错误（在某个阶段输入值可能是负的）。这个错误在许多不同的 a 和 b 值下都出现，所以我不知道该怎么继续了。

Answer 1

我觉得你可以继续用开平方的方法，但要确保一开始的值是复数，这样就能避免出现域错误。

import numpy as np
from scipy.integrate import solve_ivp

def f(t, x, a, b):
    return np.sqrt(b - a*x**2)
    # return -np.sqrt(b - a*x**2)

a = 2 + 0j
b = 1 + 0j
t0, tf = 0, 10
y0 = [1 + 0j]

res = solve_ivp(f, (t0, tf), y0, args=(a, b))  # solves successfully

不过，根据在Mathematica里的实验，这样做似乎得不到正确的结果。这样的话，这就更像是一个数学问题，你可以去数学交流网站问问。

Answer 2

带有不光滑右侧的常微分方程在数值解法下可能会出现问题。

在这种情况下，确切的解可能会有一些地方是常数，也就是x'=0，还有一些地方是解x''+ax=0，而这些状态可以在某些点上切换，这些点是两个方程同时成立的地方。这就导致了无穷多的解。

所以，最好是从一开始就选择好解的行为，这样可以解一个奇点更少的常微分方程。