Мірошниченко_Артем_ПМ-11_Лр№8_Мет_комп_експ

import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Latex
import random

def my_time_(s):
    h = s // 3600
    time = s % 3600
    m = time // 60
    s = time % 60
    if h < 10:
        h = '0' + str(int(h))
    if m < 10:
        m = '0' + str(int(m))
    if s < 10:
        s = '0' + str(int(s))

    return fr'{h}:{m}:{s}'

def to_latex(text=''): return fr'${text}$'

def d2f_dc2(c, a):
    return 2 * (-a * c + 1)**2 - 8 * c * (-a * c + 1) * a + 2 * c**2 * a**2

def f(c, a):
    return c**2 * (1 - a*c)**2

def stability_index(k, c, a, z):
    return -M * k**2 * d2f_dc2(c, a) * (c - M*z*k**4)



M = 1
a = 1
z = 0.5
c0 = 0.5

k_values = np.linspace(0, 1.2, 100)
lambda_values = []

for k in k_values:
    lambda_val = stability_index(k, c0, a, z)
    lambda_values.append(lambda_val)

plt.figure(figsize=(9, 5))
plt.plot(k_values, lambda_values)
plt.xlabel(to_latex('k'))
plt.ylabel(to_latex('λ'))
plt.title(to_latex('Залежність\ λ(k)'))
plt.grid(True)
plt.show()









M = 1
z = 0.5
k_values = np.linspace(0, 1.2, 200)

a_values = np.linspace(0.1, 2, 200)
c0_values = np.linspace(0.1, 1, 200)

stability_results = np.zeros((len(a_values), len(c0_values)))

def d2f_dc2(c, a):
    return 2 * (-a * c + 1)**2 - 8 * c * (-a * c + 1) * a + 2 * c**2 * a**2

def stability_index(k, c, a, z):
    return -M * k**2 * d2f_dc2(c, a) - M * z * k**4

for i, a in enumerate(a_values):
    for j, c0 in enumerate(c0_values):
        lambda_values = []
        for k in k_values:
            lambda_val = stability_index(k, c0, a, z)
            if lambda_val > 0:
                lambda_values.append(lambda_val)
            else:
                lambda_values.append(np.nan)
        stability_results[i, j] = np.nanmean(lambda_values)

plt.figure(figsize=(9, 5))
plt.imshow(stability_results, extent=[min(c0_values), max(c0_values), min(a_values), max(a_values)], origin='lower', aspect='auto', cmap='magma')
plt.colorbar(label=to_latex('λ'))
plt.xlabel(to_latex('c_0'))
plt.ylabel(to_latex('a'))
plt.title(to_latex('Діаграма\ стійкості\ a(c_0)'))
plt.show()









def micro_ch_pre(Nx, Ny, co, iflag):
    np.random.seed(0)
    noise = 0.02
    con = np.zeros((Nx, Ny))
    if iflag == 1:
        for i in range(Nx):
            for j in range(Ny):
                con[i, j] = co + noise * (0.5 - random.random())
    else:
        for i in range(Nx):
            for j in range(Ny):
                ii = (i - 1) * Nx + j
                con[ii] = co + noise * (0.5 - random.random())
    return con

# Parameters
Nx = 64
Ny = 64
dx = 1.0
dy = 1.0
nstep = 10000
nprint = 50
dtime = 1.0e-2
ttime = 0.0
co = 0.40
mobility = 1.0
grad_coef = 0.5

# Prepare microstructure
iflag = 1
con = micro_ch_pre(Nx, Ny, co, iflag)

lap_con = np.zeros((Nx, Ny))
lap_dummy = np.zeros((Nx, Ny))
dummy = np.zeros((Nx, Ny))

# Time and concentration arrays for plotting
time_vals = []
concentration_vals = []

con_images = []
step_list = []
# Evolve
for istep in range(1, nstep + 1):
    ttime += dtime
    for i in range(Nx):
        for j in range(Ny):
            jp = j + 1 if j + 1 < Ny else 0
            jm = j - 1 if j - 1 >= 0 else Ny - 1
            ip = i + 1 if i + 1 < Nx else 0
            im = i - 1 if i - 1 >= 0 else Nx - 1

            hne = con[ip, j]
            hnw = con[im, j]
            hns = con[i, jm]
            hnn = con[i, jp]
            hnc = con[i, j]

            lap_con_ij = (hnw + hne + hns + hnn - 4 * hnc) / (dx * dy)

            dfdcon = 2 * (con[i, j] * (1 - con[i, j])**2 - 2 * con[i, j]**2 * (1 - con[i, j]))
            dummy[i, j] = dfdcon - grad_coef * lap_con_ij

    for i in range(Nx):
        for j in range(Ny):
            jp = j + 1 if j + 1 < Ny else 0
            jm = j - 1 if j - 1 >= 0 else Ny - 1
            ip = i + 1 if i + 1 < Nx else 0
            im = i - 1 if i - 1 >= 0 else Nx - 1

            hne = dummy[ip, j]
            hnw = dummy[im, j]
            hns = dummy[i, jm]
            hnn = dummy[i, jp]
            hnc = dummy[i, j]

            lap_dummy[i, j] = (hnw + hne + hns + hnn - 4 * hnc) / (dx * dy)

            con[i, j] = con[i, j] + dtime * mobility * lap_dummy[i, j]

            if con[i, j] >= 0.9999:
                con[i, j] = 0.9999
            if con[i, j] < 0.00001:
                con[i, j] = 0.00001

    if istep % nprint == 0 or istep == 1:
        print(f"done step: {istep}")
        time_vals.append(ttime)
        concentration_vals.append(con.copy())

        if istep in [1, 100, 1000, 2500, 5000, 10000]:
            con_images.append(np.copy(con))
            step_list.append(round(ttime, 4))


fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.ravel()

for i in range(len(con_images)):
    im = axes[i].imshow(con_images[i], cmap='coolwarm')
    axes[i].set_title(to_latex(f'Time: \ {step_list[i]}'))
    plt.colorbar(im, ax=axes[i])
    axes[i].axis('off')

plt.tight_layout()
plt.show()









def micro_ch_pre(Nx, Ny, co, iflag):
    np.random.seed(0)
    noise = 0.02
    con = np.zeros((Nx, Ny))
    if iflag == 1:
        for i in range(Nx):
            for j in range(Ny):
                con[i, j] = co + noise * (0.5 - random.random())
    else:
        for i in range(Nx):
            for j in range(Ny):
                ii = (i - 1) * Nx + j
                con[ii] = co + noise * (0.5 - random.random())
    return con

# Parameters
Nx = 64
Ny = 64
dx = 1.0
dy = 1.0
nstep = 1000
nprint = 50
dtime = 1.0e-2
ttime = 0.0
co = 0.40
mobility = 1.0
grad_coef = 0.5

# Prepare microstructure
iflag = 1
con = micro_ch_pre(Nx, Ny, co, iflag)

lap_con = np.zeros((Nx, Ny))
lap_dummy = np.zeros((Nx, Ny))
dummy = np.zeros((Nx, Ny))

# Time and concentration arrays for plotting
time_vals = []
concentration_vals = []
diff_history = []

# Evolve
for istep in range(1, nstep + 1):
    ttime += dtime
    for i in range(Nx):
        for j in range(Ny):
            jp = j + 1 if j + 1 < Ny else 0
            jm = j - 1 if j - 1 >= 0 else Ny - 1
            ip = i + 1 if i + 1 < Nx else 0
            im = i - 1 if i - 1 >= 0 else Nx - 1

            hne = con[ip, j]
            hnw = con[im, j]
            hns = con[i, jm]
            hnn = con[i, jp]
            hnc = con[i, j]

            lap_con_ij = (hnw + hne + hns + hnn - 4 * hnc) / (dx * dy)

            dfdcon = 2 * (con[i, j] * (1 - con[i, j])**2 - 2 * con[i, j]**2 * (1 - con[i, j]))
            dummy[i, j] = dfdcon - grad_coef * lap_con_ij

    for i in range(Nx):
        for j in range(Ny):
            jp = j + 1 if j + 1 < Ny else 0
            jm = j - 1 if j - 1 >= 0 else Ny - 1
            ip = i + 1 if i + 1 < Nx else 0
            im = i - 1 if i - 1 >= 0 else Nx - 1

            hne = dummy[ip, j]
            hnw = dummy[im, j]
            hns = dummy[i, jm]
            hnn = dummy[i, jp]
            hnc = dummy[i, j]

            lap_dummy[i, j] = (hnw + hne + hns + hnn - 4 * hnc) / (dx * dy)

            con[i, j] = con[i, j] + dtime * mobility * lap_dummy[i, j]

            if con[i, j] >= 0.9999:
                con[i, j] = 0.9999
            if con[i, j] < 0.00001:
                con[i, j] = 0.00001

    if istep % nprint == 0 or istep == 1:
        print(f"done step: {istep}")
        average_temp = np.mean(con)
        concentration_vals.append(average_temp)
        time_vals.append(ttime)
        diff = np.mean(con**2) - np.mean(con)**2
        diff_history.append(diff)

# Plot Average Concentration vs. Time
plt.figure(figsize=(9, 5))
plt.plot(time_vals, concentration_vals)
plt.xlim(0, 5)
plt.xlabel(to_latex("Time"))
plt.ylabel(to_latex("Average\ Concentration"))
plt.title(to_latex("Average\ Concentration vs \ Time"))
plt.grid(True)
plt.show()

# Plot <c**2> - <c>**2 vs. Time
plt.figure(figsize=(9, 5))
plt.plot(time_vals, diff_history)
plt.xlabel(to_latex("Time"))
plt.ylabel(to_latex("<c^2> - <c>^2"))
plt.title(to_latex("<c^2> - <c>^2 \ vs\ Time"))
plt.grid(True)
plt.show()










# Задані параметри
M = 1
z = 0.5
k_values = np.linspace(0, 1.2, 200)

# Генерація значень a та c0
a_values1 = np.array([0.4, 0.6, 0.8, 1.2, 1.4, 1.6, 0.5, 1.0, 1.5, 1.0, 1.0])
c0_values1 = np.array([0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.5, 0.5, 0.5, 0.3, 0.7])


# Генерація значень a та c0
a_values = np.linspace(0.1, 2, 200)
c0_values = np.linspace(0.1, 1, 200)

# Створення масиву для результатів
stability_results = np.zeros((len(a_values), len(c0_values)))

# Визначення функції d2f_dc2
def d2f_dc2(c, a):
    return 2 * (-a * c + 1)**2 - 8 * c * (-a * c + 1) * a + 2 * c**2 * a**2

# Визначення функції f(c, a)
def f(c, a):
    return c**2 * (1 - a * c)**2

# Визначення показника стійкості lambda
def stability_index(k, c, a, z):
    return -M * k**2 * d2f_dc2(c, a) - M * z * k**4

# Побудова діаграми стійкості
for i, a in enumerate(a_values):
    for j, c0 in enumerate(c0_values):
        lambda_values = []
        for k in k_values:
            lambda_val = stability_index(k, c0, a, z)
            if lambda_val > 0:
                lambda_values.append(lambda_val)
            else:
                lambda_values.append(np.nan)
        stability_results[i, j] = np.nanmean(lambda_values)

# Побудова графіка
plt.figure(figsize=(9, 5))
plt.imshow(stability_results, extent=[min(c0_values), max(c0_values), min(a_values), max(a_values)], origin='lower', aspect='auto', cmap='magma')
plt.colorbar(label=to_latex('λ'))
plt.xlabel(to_latex('c_0'))
plt.scatter(c0_values1, a_values1, color='red')
plt.ylabel(to_latex('a'))
plt.title(to_latex('Діаграма\ стійкості\ a(c_0)'))
plt.show()









def micro_ch_pre(Nx, Ny, co, iflag):
    np.random.seed(0)
    noise = 0.02
    con = np.zeros((Nx, Ny))
    if iflag == 1:
        for i in range(Nx):
            for j in range(Ny):
                con[i, j] = co + noise * (0.5 - random.random())
    else:
        for i in range(Nx):
            for j in range(Ny):
                ii = (i - 1) * Nx + j
                con[ii] = co + noise * (0.5 - random.random())
    return con

# Parameters
Nx = 64
Ny = 64
dx = 1.0
dy = 1.0
nstep = 25000
nprint = 50
dtime = 1.0e-2
ttime = 0.0
mobility = 1.0
grad_coef = 0.5

# Values to iterate over
a_values = np.array([0.4, 0.6, 0.8, 1.2, 1.4, 1.6, 0.5, 1.0, 1.5, 1.0, 1.0])
co_values = np.array([0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.5, 0.5, 0.5, 0.3, 0.7])

for a, co in zip(a_values, co_values):
    print(f"\n\nRunning for a = {a} and co = {co}\n\n")
    ttime = 0
    con = micro_ch_pre(Nx, Ny, co, iflag=1)

    lap_con = np.zeros((Nx, Ny))
    lap_dummy = np.zeros((Nx, Ny))
    dummy = np.zeros((Nx, Ny))

    # Time and concentration arrays for plotting
    time_vals = []
    concentration_vals = []

    # Evolve
    for istep in range(1, nstep + 1):
        ttime += dtime
        for i in range(Nx):
            for j in range(Ny):
                jp = j + 1 if j + 1 < Ny else 0
                jm = j - 1 if j - 1 >= 0 else Ny - 1
                ip = i + 1 if i + 1 < Nx else 0
                im = i - 1 if i - 1 >= 0 else Nx - 1

                hne = con[ip, j]
                hnw = con[im, j]
                hns = con[i, jm]
                hnn = con[i, jp]
                hnc = con[i, j]

                lap_con_ij = (hnw + hne + hns + hnn - 4 * hnc) / (dx * dy)

                dfdcon =  2*con[i, j]*(-a*con[i, j]+1)**2-2*((con[i, j])**2)*(-a*con[i, j]+1)*a
                dummy[i, j] = dfdcon - grad_coef * lap_con_ij

        for i in range(Nx):
            for j in range(Ny):
                jp = j + 1 if j + 1 < Ny else 0
                jm = j - 1 if j - 1 >= 0 else Ny - 1
                ip = i + 1 if i + 1 < Nx else 0
                im = i - 1 if i - 1 >= 0 else Nx - 1

                hne = dummy[ip, j]
                hnw = dummy[im, j]
                hns = dummy[i, jm]
                hnn = dummy[i, jp]
                hnc = dummy[i, j]

                lap_dummy[i, j] = (hnw + hne + hns + hnn - 4 * hnc) / (dx * dy)

                con[i, j] = con[i, j] + dtime * mobility * lap_dummy[i, j]

                if con[i, j] >= 0.9999:
                    con[i, j] = 0.9999
                if con[i, j] < 0.00001:
                    con[i, j] = 0.00001

        if istep % nprint == 0 or istep == 1:
            if istep % 5000 == 0:
                print(f"done step: {istep}")
            time_vals.append(ttime)
            concentration_vals.append(con.copy())

            if istep == nstep:
                plt.figure(figsize=(9, 5))
                plt.imshow(con, cmap='coolwarm')
                plt.colorbar()
                plt.title(to_latex(f"a={a},\ c_0={co},\ Time: {round(ttime, 4)}"))
                plt.show()