Мірошниченко_ГЙМ_ПР№11-13

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import cumulative_trapezoid, trapezoid

TAU = 1.5
ZETA_E = 0.0
I_S_ADDITIVE = 0.0

CURVES_PARAMS = [
    {"label": "1: SOC (I_sigma=0, I_zeta=30)",    "I_sigma": 0.0, "I_zeta": 30.0, "style": "k-",  "marker": None},
    {"label": "2: A+N (I_sigma=0.5, I_zeta=30)",  "I_sigma": 0.5, "I_zeta": 30.0, "style": "r-",  "marker": "o"},
    {"label": "3: N (I_sigma=1, I_zeta=5)",       "I_sigma": 1.0, "I_zeta": 5.0,  "style": "g--", "marker": "^"},
    {"label": "4: A (I_sigma=2, I_zeta=0.5)",     "I_sigma": 2.0, "I_zeta": 0.5,  "style": "b-.", "marker": "s"}
]

S_MIN_LOG = -4
S_MAX_LOG = 0
POINTS = 2000

def d_tau(s):
    s = np.abs(s)
    return (1 + s**TAU)**(-1)

def func_f(s):
    term2 = ZETA_E * (np.abs(s)**(TAU/2)) * d_tau(s)
    return -s + term2

def func_I(s, I_sigma, I_zeta):
    s = np.abs(s)
    dt = d_tau(s)
    term_brackets = I_sigma + I_zeta * (s**TAU)
    return I_S_ADDITIVE + term_brackets * (dt**2)

plt.figure(figsize=(10, 8))

s_values = np.logspace(S_MIN_LOG, S_MAX_LOG, POINTS)
s_calc = np.insert(s_values, 0, 0.0)

for params in CURVES_PARAMS:
    I_sigma = params["I_sigma"]
    I_zeta = params["I_zeta"]

    I_vals = func_I(s_calc, I_sigma, I_zeta)

    with np.errstate(divide='ignore', invalid='ignore'):
        G_calc = func_f(s_calc) / I_vals

    if np.isinf(G_calc[0]) or np.isnan(G_calc[0]):
        G_calc[0] = 0

    integral_values = cumulative_trapezoid(G_calc, s_calc, initial=0)
    integral_res = integral_values[1:]

    I_vals_res = func_I(s_values, I_sigma, I_zeta)
    P_raw = (1.0 / I_vals_res) * np.exp(integral_res)

    norm_const_Z = trapezoid(P_raw, s_values)
    P_normalized = P_raw / norm_const_Z

    plt.plot(s_values, P_normalized, params["style"], label=params["label"],
             linewidth=2, marker=params["marker"], markevery=150, markersize=8)

plt.xscale('log')
plt.yscale('log')
plt.xlim(1e-4, 1e0)
plt.ylim(1e-3, 1e4)
plt.xlabel('Розмір лавини, s', fontsize=12)
plt.ylabel('Функція розподілу, P(s)', fontsize=12)
plt.title('Функції розподілу', fontsize=14)
plt.legend(loc="best")
plt.grid(True, which="both", ls="-", alpha=0.4)
plt.tight_layout()


DT = 0.01
STEPS = 5000
times = np.arange(STEPS) * DT

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

np.random.seed(42)

for i, params in enumerate(CURVES_PARAMS):
    ax = axes[i]
    I_sigma = params["I_sigma"]
    I_zeta = params["I_zeta"]
    label = params["label"]

    s_t = np.zeros(STEPS)
    s_t[0] = 0.1

    for t in range(STEPS - 1):
        s_curr = s_t[t]

        I_val = func_I(s_curr, I_sigma, I_zeta)
        f_val = func_f(s_curr)

        drift = f_val * DT
        diffusion = np.sqrt(2 * I_val * DT) * np.random.normal(0, 1)

        s_next = s_curr + drift + diffusion

        if s_next < 1e-6:
            s_next = 1e-6

        s_t[t+1] = s_next

    color = params["style"][0]
    if color == 'k': color = 'black'
    elif color == 'r': color = 'red'
    elif color == 'g': color = 'green'
    elif color == 'b': color = 'blue'

    ax.plot(times, s_t, color=color, linewidth=0.8)
    ax.set_title(f"Варіант {label}", fontsize=10)
    ax.set_xlabel("Час, t")
    ax.set_ylabel("s(t)")
    ax.grid(True, alpha=0.3)
    ax.set_ylim(bottom=0)

plt.suptitle("Часові ряди стохастичної змінної", fontsize=16)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()