Numeryczne rozwiązywanie ODE w Pythonie

Question

W jaki sposób mogę rozwiązać problem ODE w języku Python?

Rozważmy

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

z następujących warunków

u(0) = 1.49907 

i

\dot{u}(0) = 0 

z ograniczeniem

0 <= \phi <= 7\pi. 

W końcu chcę utworzyć wykres parametryczny, w którym współrzędne x i y są generowane w funkcji u.

Problem polega na tym, że muszę uruchomić dwa razy odiuint, ponieważ jest to równanie różniczkowe drugiego rzędu. Próbowałem go uruchomić ponownie po raz pierwszy, ale wraca z błędem jakobian. Musi istnieć sposób na uruchomienie go dwa razy na raz.

Tutaj jest błąd:

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

który poniżej generuje kod. Linia, o której mowa, to sol = odint.

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

Edit

próbuję skonstruować fabułę na stronie 9

Classical Mechanics Taylor

Oto działka z 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() integruje ODE. Czy tego właśnie szukasz?

Answer 3

Można użyć scipy.integrate.ode. Aby rozwiązać dy/dt = f (t, y), ze warunek początkowy y (t0) = y0, w czasie = T1 z 4 rzędu Runge-Kutta mógłby zrobić coś takiego:

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 

Edycja: You Muszę uzyskać twoją pochodną do pierwszego rzędu, aby użyć integracji numerycznej. Można to osiągnąć, ustawiając np. z1 = u i z2 = du/dt, po których masz dz1/dt = z2 i dz2/dt = d^2u/dt^2. Zastąp je oryginalnym równaniem i po prostu wykonaj iteracje po wektorach dZ/dt, czyli w pierwszym rzędzie.

Edit 2: Oto przykładowy kod do całej sprawy:

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

Kod z innych question jest bardzo blisko do tego, co chcesz. Potrzebne są dwie zmiany:

• Byłaś rozwiązywania inny ODE (bo zmieniłeś dwa znaki wewnątrz funkcji deriv)
• y składnikiem żądanej działki pochodzi z wartościami rozwiązania, a nie od wartości pierwszego pochodną rozwiązania, dlatego należy zastąpić u[:,0] (wartości funkcji) dla u[:, 1] (pochodne).

Jest to końcowy wynik:

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

Jednak uważam, że należy użyć kodu z odpowiedzią unutbu, ponieważ to samo dokumentowanie (u, udot = z) i wykorzystuje np.linspace zamiast np.arange. Następnie uruchom to, aby uzyskać żądaną liczbę:

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