Как решить это 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 et c. Он имеет ошибку типа

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

Но если пусть n является целым числом, он будет работать без проблем. Кто-нибудь может мне помочь?

Answer 1

Рассмотрим проблему эрзаца

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

, которая дает одно и то же решение, пока u остается положительным, что составляет примерно E=8.

Answer 2

Проблема в том, что u становится отрицательным во время решения. Для целочисленных значений n это хорошо, но для нецелого n возведение в степень должно выполняться либо:

1. Представление показателя в виде сокращенной дроби с использованием знаменателя- порядок root основания и возведение результата в числитель.
2. Используя такой подход, как в этот ответ на Математический стек обмена

Для Отрицательное основание, результат в первом случае будет мнимым, если знаменатель четен и действителен, только если он нечетен, а результат во втором случае всегда будет мнимым.

Это простое решение заключается в инициализации P0 с мнимыми компонентами,

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

---

Примечание - В Python 2.7 функция pow должна для возведения в степень, в Python 3.x можно использовать либо оператор **, либо функцию pow.