如何解决这个 ODE？

Question

无论是用odeint还是solve_ivp求解，我都遇到了这个方程式的问题。

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def ODE(E, p):
    u, v = p
    n = 3.25
    dudE = v
    dvdE = -(u**n)-(2*v/E)
    return [dudE, dvdE]

P0 = [1,0]
solve = solve_ivp(ODE, (0.001,10), P0, t_eval=np.linspace(0.001,10,500))

我不能让 n=3.25 或 n=0.25 等。它有一个像

这样的错误

Project_q2d.py:19: RuntimeWarning: invalid value encountered in double_scalars
  dvdE = -(u**n)-(2*v/E)

但是如果让 n 是一个整数，它将运行完美地没有任何问题。谁能帮帮我？

Answer 1

问题是u在求解过程中变成了负数。对于 n 的整数值，这很好，但对于非整数 n，取幂必须由

之一执行

将指数表示为约分式，取底数的分母阶根并将结果提高到分子。
使用类似于 this answer on Mathematics Stack Exchange

对于负底，如果分母是偶数，则第一种情况的结果将是虚数，只有当分母为奇数时才为实数，而第二种情况的结果将始终是虚数。

这个简单的解决方案是用虚部初始化 P0，

P0 = [1 + 0j, 0 + 0j]

注意 - 在Python 2.7中，pow函数必须用于求幂，在Python 3.x中要么可以使用 ** 运算符或 pow 函数。

Answer 2

考虑 ersatz 问题

u'' + 2*u'/E + u*abs(u)^(n-1)

只要 u 保持正数，它就会给出相同的解决方案，恰好达到大约 E=8。

如何解决这个 ODE？

How to solve this ODE?

python

numpy

scipy

ode