Как решить три квадратичных дифференциальных уравнения в Python? - PullRequest
Новая
       71

Как решить три квадратичных дифференциальных уравнения в Python?

0 голосов
/ 27 мая 2018

Я только начал использовать Python для научного рисования для построения численных решений дифференциальных уравнений.Я знаю, как использовать модули для решения и построения дифференциальных уравнений, но понятия не имею о системах дифференциальных уравнений.Как можно построить следующую связанную систему?

Моя система дифференциальных уравнений:

dw/dx=y и
dy/dx=-a-3*H*y и
dz/dx=-H*(1+z)

что a = 0.1 и H=sqrt((1+z)**3+w+u**2/(2*a))

И мой код: 

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

def model(w,y,z,x,H): 
    dwdx=y
    dydx=-a-3*H*y
    dzdx=-H*(1+z)
    a=0.1
    H=sqrt((1+z)**3+w+u**2/(2*a))
    return [w,y,z,H]   

z0=1100      #initial condition
w0=-2.26e-8
y0=-.38e-4
H0=36532.63
b=0
c=10000
x=N.arange(b,c,0.01)

y=odeint(model,y0,x)    #f=Function name that returns derivative values at requested y and t values as dydt = f(y,t) 
w=odeint(model,w0,x)
z=odeint(model,z0,x)

plt.plot(w,x)
plt.plot(y,x)
plt.plot(z,x)
plt.legend(loc='best')
plt.show()

1 Ответ

0 голосов
/ 27 мая 2018

Интеграторы ODE общего назначения ожидают, что динамическая система сведена к абстрактной системе первого порядка.Такая система имеет векторное пространство состояний, а дифференциальное уравнение обеспечивает векторы скорости для этого пространства.Здесь состояние имеет 3 скалярных компонента, которые дают трехмерный вектор в качестве состояния.Если вы хотите использовать компоненты отдельно, первым шагом в функции ODE является извлечение этих компонентов из вектора состояния, а последним шагом является составление вектора возврата из производных компонентов в правильном порядке.

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

def model(u,t): 
    w, y, z = u
    a=0.1
    H=sqrt((1+z)**3+w+u**2/(2*a))
    dwdx=y
    dydx=-a-3*H*y
    dzdx=-H*(1+z)
    return [dwdx, dydx, dzdx]

, а затем вызвать интегратор один раз с комбинированным начальным состоянием 

u0 = [ w0, y0, z0]
u = odeint(model, u0, x)

w,y,z = u.T

Пожалуйста, проверьте также аргументыфункции графика, общая схема plot(x,y).

...