How to solve this ODE?

user_13107678 0
Viewed 92

I have encountered a problem in this equation whether I use odeint or solve_ivp to solve. 

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))

I cannot let the n=3.25 or n=0.25 etc. It have an error like

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

But if let n is an integer it will run perfectly without any problem. Can anyone help me?

ode scipy numpy python
2 Answers

The problem is that u becomes negative during the solving. For integer values of n this is fine, but for a non-integer n the exponentiation will have to be performed by either

  1. Representing the exponent as a reduced fraction, taking a denominator-order root of the base and raising the result to the numerator.
  2. Using an approach like that in this answer on Mathematics Stack Exchange

For a negative base, the result in the first case will be imaginary if the denominator is even and real only if it is odd and the result in the second case will always be imaginary.

This simple solution is to initialize P0 with imaginary components, 

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

Note  -  In Python 2.7 the pow function must be used for the exponentiation, in Python 3.x either the ** operator or the pow function can be used.

user_106599100

Consider the ersatz problem

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

which gives the same solution as long as u stays positive, which happens to be up to about E=8.

enter image description here

user_30881380
Related