Лабы / lab4

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

var = 13

# Функция, отвечающая за расчет ОДУ:
# Y - вектор переменных состояния,
# t - время,
# var - вариант
def ode(Y, t, b, c):
#1
    x, y = Y

# 1.1
    #dydt = [-4*x-4*y, 1.5*x + y]
    #dydt = [-4 * x + 0.1 * var * x**2 - 4 * y, 1.5 * x + y - 0.2 * var * y**3] #1 (не линеариз)
# 1.2
    # dydt = [1 * x + 0.5 * y, 0.5 * x + y]
    #dydt = [x + 0.5 * y - 0.1 * var * y**2, 0.5 * x - 0.2 * var * x**2 + y] #2 (не линеариз)
# 1.3
    # dydt = [2 * x + y, x - 3 * y]
    #dydt = [2 * x + 0.2 * var * x**2 + y - 0.1 * var * y**2, x - 3 * y]
# 1.4
    # dydt = [2 * y, -3 * x - y]
    #dydt = [-0.1 * var * x**2 + 2 * y, -3 * x - y]
# 1.5
    # dydt = [0.1 * x - 4 * y, 4 * x + 0.1 * y]
    #dydt = [0.1 * x - 4 * y, 4 * x - 0.2 * var * x**2 + 0.1 * y]
# 1.6
    # dydt = [x - 4 * y, 2 * x - y]
    #dydt = [x - 0.1 * var * x**2 - 4 * y + 0.3 * var * y**2, 2 * x + 0.2 * var * x**2 - y - 0.3 * var * y**3]



# 2
x, y = Y
def dead_zone_scalar(x, width=0.5):
    if np.abs(x) < width:
        return 0
    elif x > 0:
        return x - width
    else:
        return x + width

dead_zone = np.vectorize(dead_zone_scalar, otypes = [np.float], excluded = ['width'])
# dydt = [y, -y - np.sign(dead_zone(x, 0.2 * var + 0.2))]




# 3
m = 0.2 * var # масса
l = 5 / var # длинна маятника
b = 0.1 + 0.015 * var # коэффициент трения
theta, omega = Y
# 3.1
# dydt = [omega, - m * 9.81 * l * np.sin(theta)]
# 3.2
# dydt = [omega, -b * omega - m * 9.81 * l * np.sin(theta)]

#4
x, y = Y
# U = var # 4.1
U = var / 2 # 4.2
# U = 2 * var # 4.3
dydt = [y, U * (1 - x * x) * y - x]
return dydt

# Функция для получения решения ОДУ с заданными начальными условиями:
# args - параметры ОДУ (см. определение функции с ОДУ),
# y0 - начальные условия для первой переменной состояния,
# dy0 - начальные условия для второй переменной состояния
# (или в нашем случае ее производной),
# ts - длительность решения,
# nt - число шагов в решении (=время интегрирования*шаг времени).
def calcODE(args, y0, dy0, ts = 10, nt = 101):
    y0 = [y0, dy0]
    t = np.linspace(0, ts, nt)
    sol = odeint(ode, y0, t, args)
    return sol

# args - парметры ОДУ (см. шаг 1),
# deltaX - шаг начальных условий по горизонтальной оси
# (переменной состояния)),
# deltaDX - шаг начальных условий по вертикальной оси
# (производной переменной состояния),
# startX - начальное значение интервала начальных условий,
# stopX - конечное значение интервала начальных условий,
# startDX - начальное значение интервала начальных условий,
# stopDX - конечное значение интервала начальных условий,
# ts - длительность решения,
# nt - число шагов в решении (= время интегрирования * шаг времени).
def drawPhasePortrait(args, deltaX = 1, deltaDX = 1, startX = 0, stopX = 5, startDX = 0, stopDX = 5, ts = 10, nt = 101):
    for y0 in range(startX, stopX, deltaX):
        for dy0 in range(startDX, stopDX, deltaDX):
            sol = calcODE(args, y0, dy0, ts, nt)
            plt.plot(sol[:, 0], sol[:, 1], 'b')
            plt.xlabel('x')
            plt.ylabel('dx/dt')
            plt.grid()
            plt.show()

#Задание 1. Изучение фазовых портретов типовых особых точек
b = 0
c = 0
args=(b, c)
# drawPhasePortrait(args, 1, 1, -5, 5, -4, 4, ts = 0.5, nt = 301)
#Задание 2. Построение фазовых портретов для нелинейных систем
#Задание 3. Посторение фазового портерта для маятника
#Задание 4. Построение фазового портетра осциллятора Ван дер Поля
#Задание 5. Построение фазового портрета аттрактора Лоренца (3 порядок)
drawPhasePortrait(args, 1, 1, -5, 5, -4, 4, ts = 20, nt = 301)