Мірошниченко_КМЗПМ_ЛР№6

import math
import numpy as np
import matplotlib.pyplot as plt
from numba import njit

@njit(cache=True, fastmath=True)
def step_f(
    phi,
    x,
    phi0,
    new_phi,
    new_x,
    dt,
    dx,
    eps,
    omega,
    lam,
    F,
):
    n = phi.shape[0]
    inv_dx2 = 1.0 / (dx * dx)

    for i in range(n):
        im = (i - 1) % n
        ip = (i + 1) % n
        for j in range(n):
            jm = (j - 1) % n
            jp = (j + 1) % n

            lap_phi = (
                phi[ip, j] + phi[im, j] + phi[i, jp] + phi[i, jm] - 4.0 * phi[i, j]
            ) * inv_dx2

            phase = math.pi * (phi[i, j] - phi0[i, j])
            new_phi[i, j] = phi[i, j] + dt * (
                (omega * omega) * lap_phi
                + math.sin(phase)
                + lam * x[i, j] * (1.0 + math.cos(phase))
            )

    for i in range(n):
        im = (i - 1) % n
        ip = (i + 1) % n
        for j in range(n):
            jm = (j - 1) % n
            jp = (j + 1) % n

            lap_x = (
                x[ip, j] + x[im, j] + x[i, jp] + x[i, jm] - 4.0 * x[i, j]
            ) * inv_dx2

            valx = x[i, j]

            reaction = valx * math.exp(-eps * valx)

            diff_factor = 1.0 - eps * valx * (1.0 - valx)
            diffusion = diff_factor * lap_x

            new_x[i, j] = (
                x[i, j]
                + dt * (F - reaction + diffusion)
                - 0.5 * (new_phi[i, j] - phi[i, j])
            )

def simulate_numba(
    N=256,
    dx=1.0,
    dt=0.01,
    t_end=500.0,
    eps=4.0,
    D0=1.0,
    omega=2.0,
    lam=10.0,
    F=2.0,
    x0=0.40,
    phi0_level=0.0,
    noise_amp_x=0.01,
    noise_amp_phi=0.01,
    seed=0,
    snapshot_times=(0.0, 20.0, 40.0, 60.0, 100.0, 500.0),
):
    rng = np.random.default_rng(seed)
    x = (x0 + noise_amp_x * rng.standard_normal((N, N))).astype(np.float64)
    phi = (phi0_level + noise_amp_phi * rng.standard_normal((N, N))).astype(np.float64)

    Phi0 = phi.copy()

    new_phi = np.empty_like(phi)
    new_x = np.empty_like(x)

    n_steps = int(round(t_end / dt))
    snap_steps = {int(round(t / dt)): float(t) for t in snapshot_times}
    snaps = {0.0: phi.copy()}

    progress_every = max(1, n_steps // 10)

    for step in range(1, n_steps + 1):
        t = step * dt

        step_f(
            phi, x, Phi0, new_phi, new_x,
            dt, dx, eps, omega, lam, F,
        )

        phi, new_phi = new_phi, phi
        x, new_x = new_x, x

        if step in snap_steps:
            snaps[snap_steps[step]] = phi.copy()

        if step % progress_every == 0:
            print(
                f"t={t:7.2f} | min/max phi={phi.min(): .4f}/{phi.max(): .4f} | "
                f"min/max x={x.min(): .4f}/{x.max(): .4f}"
            )

        if not (np.isfinite(phi).all() and np.isfinite(x).all()):
            raise FloatingPointError("NaN/Inf detected")

    for tt in snapshot_times:
        if float(tt) not in snaps:
            raise RuntimeError(f"Snapshot at t={tt} not recorded (check dt).")

    return snaps

def plot_snapshots(phi_snaps, dx):
    times = sorted(phi_snaps.keys())
    cols = 2
    rows = int(math.ceil(len(times) / cols))

    N = next(iter(phi_snaps.values())).shape[0]
    xs = np.arange(N) * dx
    ys = np.arange(N) * dx
    X, Y = np.meshgrid(xs, ys)

    fig = plt.figure(figsize=(10, 12))
    for i, tt in enumerate(times, start=1):
        ax = fig.add_subplot(rows, cols, i, projection="3d")
        Z = phi_snaps[tt]
        surf = ax.plot_surface(
            X, Y, Z, rstride=2, cstride=2,
            cmap="coolwarm", linewidth=0, antialiased=True
        )
        ax.set_title(f"t = {int(tt)}")
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_zlabel("ϕ")
        ax.view_init(elev=25, azim=-60)
        fig.colorbar(surf, ax=ax, shrink=0.55, pad=0.05)

    fig.suptitle("ϕ(x,y,t)", fontsize=14)
    plt.tight_layout()
    plt.show()

def main():
    snaps = simulate_numba(
        N=256, dx=1.0, dt=0.001, t_end=500.0,
        eps=4.0, omega=2.0, lam=10.0, F=2.0
    )
    plot_snapshots(snaps, dx=1.0)

if __name__ == "__main__":
    main()