Численное решение оду в Python

Question

Численное решение оду в Python

Как численно решить оду в Python?

Рассмотрим

уравнение для решения

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

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

u(0) = 1.49907

И

dot{u}(0) = 0

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

0 <= phi <= 7pi.

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

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

Вот ошибка:

Odepack.ошибка: функция и ее Якобиан должны быть вызываемыми функциями

Который генерирует приведенный ниже код. Рассматриваемая линия-это 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()

Edit

Я пытаюсь построить сюжет на странице 9

Классическая Механика Тейлор

Вот сюжет с 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}}]

1793 4

python plot numerical-methods differential-equations

4 ответов:

Comments

Ничего не найдено.

unutbu · Accepted Answer · 2013-04-12 04:03:58

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

Введите описание изображения здесь

jorgeca · Accepted Answer · 2017-05-23 15:32:06

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

вы решали другую оду (потому что вы изменили два знака внутри функции 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()
Тем не менее, я предлагаю вам использовать код из ответа унутбу, потому что он сам документирует (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()

HenriV · Accepted Answer · 2013-05-03 16:33:32

Вы можете использовать scipy.интегрировать.ода. Для решения задачи dy/dt = f(t, y), с начальным условием y (t0)=y0, в момент времени=t1 с помощью Рунге-Кутты 4-го порядка можно сделать примерно следующее:
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
Edit: вы должны получить производную первого порядка, чтобы использовать численное интегрирование. Этого можно достичь, установив, например, 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()

Bitwise · Accepted Answer · 2013-04-10 17:41:54

Сципион.integrate () делает интеграцию ODE. Это то, что вы ищете?