Численное решение ODE в Python

Question

Как я могу численно решить ODE в Python?

Рассмотрим

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

со следующими условиями

u(0) = 1.49907 

и

\dot{u}(0) = 0 

с ограничением

0 <= \phi <= 7\pi. 

Затем, наконец, я хочу создать параметрический график, в котором координаты x и y генерируются как функция u.

Проблема в том, что мне нужно запустить odeint дважды, так как это дифференциальное уравнение второго порядка. Я попытался запустить его после первого раза, но он вернулся с ошибкой якобиана. Должен быть способ запустить его дважды все сразу.

Здесь ошибка:

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

, который ниже генерирует код. Линия, о которой идет речь, - золь = 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() 

Редактировать

Я пытаюсь построить график на странице 9

Classical Mechanics Taylor

Вот сюжет с 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() Интеграция ODE. Это то, что вы ищете?

Answer 3

Вы можете использовать scipy.integrate.ode. Для решения ду/дТ = F (T, Y), с начальным условием у (t0) = y0, в момент времени = t1 с четвёртого порядка Рунге-Кутта вы могли бы сделать что-то вроде этого:

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 

Редактировать: должны получить ваш производный для первого порядка использования численного интегрирования. Это можно достичь, установив, например, z1 = u и z2 = du/dt, после чего у вас есть dz1/dt = z2 и dz2/dt = d^2u/dt^2. Подставим их в исходное уравнение и просто перейдем к вектору dZ/dt, который является первым порядком.

Edit 2: Вот пример кода для всего этого:

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

код с другого question очень близко к тому, что вы хотите. Необходимы два изменения:

• Вы решали различную ОДУ (потому что вы изменили два знак внутри функций deriv)
• y компонента нужного участка происходит от значений решения, а не от значений первого производная от решения, поэтому вам нужно заменить u[:,0] (значения функции) на u[:, 1] (производные).

Это конечный результат:

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

Однако, я полагаю, что вы используете код из ответа unutbu, потому что это документирован (u, udot = z) и использует np.linspace вместо np.arange. Затем запустите это, чтобы получить нужную цифру:

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