Risoluzione ODE numerica in Python

Question

Come risolvo numericamente un ODE in Python?

consideri

\ddot{u}(\phi) = -u + \sqrt{u} 

con le seguenti condizioni

u(0) = 1.49907 

e

\dot{u}(0) = 0 

con il vincolo

0 <= \phi <= 7\pi. 

Infine, voglio produrre un grafico parametrico in cui le coordinate xey sono generate come una funzione di u.

Il problema è che ho bisogno di eseguire due volte odeint poiché questa è un'equazione differenziale del secondo ordine. Ho provato a farlo funzionare di nuovo dopo la prima volta ma torna con un errore Jacobiano. Ci deve essere un modo per eseguirlo due volte in una volta.

Ecco l'errore:

odepack.error: The function and its Jacobian must be callable functions

che il codice qui sotto genera. La linea in questione è il sol = odeint.

import numpy as np 
from scipy.integrate import odeint 
import matplotlib.pyplot as plt 
from numpy import linspace 


def f(u, t): 
    return -u + np.sqrt(u) 


times = linspace(0.0001, 7 * np.pi, 1000) 
y0 = 1.49907 
yprime0 = 0 
yvals = odeint(f, yprime0, times) 

sol = odeint(yvals, y0, times) 

x = 1/sol * np.cos(times) 
y = 1/sol * np.sin(times) 

plot(x,y) 

plt.show() 

Modifica

Sto cercando di costruire la trama a pagina 9

Classical Mechanics Taylor

Ecco la trama con Mathematica

In[27]:= sol = 
NDSolve[{y''[t] == -y[t] + Sqrt[y[t]], y[0] == 1/.66707928, 
    y'[0] == 0}, y, {t, 0, 10*\[Pi]}]; 

In[28]:= ysol = y[t] /. sol[[1]]; 

In[30]:= ParametricPlot[{1/ysol*Cos[t], 1/ysol*Sin[t]}, {t, 0, 
    7 \[Pi]}, PlotRange -> {{-2, 2}, {-2.5, 2.5}}] 

Answer 1

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

pi = np.pi 
sqrt = np.sqrt 
cos = np.cos 
sin = np.sin 

def deriv_z(z, phi): 
    u, udot = z 
    return [udot, -u + sqrt(u)] 

phi = np.linspace(0, 7.0*pi, 2000) 
zinit = [1.49907, 0] 
z = integrate.odeint(deriv_z, zinit, phi) 
u, udot = z.T 
# plt.plot(phi, u) 
fig, ax = plt.subplots() 
ax.plot(1/u*cos(phi), 1/u*sin(phi)) 
ax.set_aspect('equal') 
plt.grid(True) 
plt.show() 

Answer 2

scipy.integrate() esegue l'integrazione ODE. È quello che stai cercando?

Answer 3

È possibile utilizzare scipy.integrate.ode. Per risolvere dy/dt = f (t, y), con condizione iniziale y (t0) = y0, al tempo = T1 con 4 ° ordine Runge-Kutta si potrebbe fare qualcosa di simile:

from scipy.integrate import ode 
solver = ode(f).set_integrator('dopri5') 
solver.set_initial_value(y0, t0) 
dt = 0.1 
while t < t1: 
    y = solver.integrate(t+dt) 
    t += dt 

Modifica: devi ottenere il tuo derivato al primo ordine per utilizzare l'integrazione numerica. Questo è possibile ottenere impostando per es. z1 = u e z2 = du/dt, dopo di che hai dz1/dt = z2 e dz2/dt = d^2u/dt^2. Sostituiscili nella tua equazione originale e semplicemente esegui un'iterazione sul vettore dZ/dt, che è il primo ordine.

Edit 2: Ecco un codice di esempio per l'intera faccenda:

import numpy as np 
import matplotlib.pyplot as plt 

from numpy import sqrt, pi, sin, cos 
from scipy.integrate import ode 

# use z = [z1, z2] = [u, u'] 
# and then f = z' = [u', u''] = [z2, -z1+sqrt(z1)] 
def f(phi, z): 
    return [z[1], -z[0]+sqrt(z[0])] 


# initialize the 4th order Runge-Kutta solver 
solver = ode(f).set_integrator('dopri5') 

# initial value 
z0 = [1.49907, 0.] 
solver.set_initial_value(z0) 

values = 1000 
phi = np.linspace(0.0001, 7.*pi, values) 
u = np.zeros(values) 

for ii in range(values): 
    u[ii] = solver.integrate(phi[ii])[0] #z[0]=u 

x = 1./u * cos(phi) 
y = 1./u * sin(phi) 

plt.figure() 
plt.plot(x,y) 
plt.grid() 
plt.show() 

Answer 4

Il codice dal vostro altro question è molto vicino a ciò che si desidera. Sono necessari due cambiamenti:

• Stavi risolvere un ODE diversa (perché hai cambiato due segni funzione deriv all'interno)
• La componente y del vostro diagramma desiderato viene dai valori di soluzione, non dai valori della prima derivato della soluzione, quindi è necessario sostituire u[:,0] (valori di funzione) per u[:, 1] (derivati).

Questo è il risultato finale:

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

def deriv(u, t): 
    return np.array([u[1], -u[0] + np.sqrt(u[0])]) 

time = np.arange(0.01, 7 * np.pi, 0.0001) 
uinit = np.array([1.49907, 0]) 
u = odeint(deriv, uinit, time) 

x = 1/u[:, 0] * np.cos(time) 
y = 1/u[:, 0] * np.sin(time) 

plt.plot(x, y) 
plt.show() 

Tuttavia, suggerisco di utilizzare il codice di risposta di unutbu perché è auto documentazione (u, udot = z) e usa np.linspace invece di np.arange. Quindi, esegui questo per ottenere la cifra desiderata:

x = 1/u * np.cos(phi) 
y = 1/u * np.sin(phi) 
plt.plot(x, y) 
plt.show()