Cahn1ddiflog

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import rfft, irfft, rfftfreq
import os
def cahnexp1d(cc, c1, c2):
    cc1 = np.copy(cc)*c1
    cc2 = np.copy(cc)*c2
    cc3 = np.copy(cc**3)*c1
    for i in range(l):
            n = i - 2
            cc[n] =   (   2*(cc1[n] - cc3[n]) + (cc3[n-1] + cc3[n+1] -1*(cc1[n-1] + cc1[n+1]))
                            -1*(cc2[n-2] + cc2[n+2] + 6*cc2[n])
                            +(4*cc2[n-1] + 4*cc2[n+1])
                        ) + cc[n]
    return cc

def cahnfourier1d(cc, k2, k4):
    cct = rfft(cc)
    cct3 = rfft(cc**3 - cc)
    cct = cct + difd*dt*(-k2*(cct3) - k4*cct)
    cc = irfft(cct)
    return cc

np.random.seed(4)

xmin = 0

xmax = 1

tmin = 0

tmax = 0.5

dt = 1/22000000

dx = 1/128

difd = 1

gamma = (3.4/128)**2

l = int((xmax/dx))

k  = rfftfreq(l, dx/(2*np.pi))

c1 = difd*dt/(dx**2)

c2 = gamma*c1/(dx**2)

k2 = k**2

k4 = (k2**2)*gamma

ccexp = np.random.rand(l)*2 - 1

ccfourier = np.copy(ccexp)

t = 0

tamanho = l/2

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

plt.yscale("log")

os.makedirs(f".\\1d", exist_ok = True)

dif1, dif2, dif3, dif4 = [], [], [], []

tt = np.arange(dt, tmax + dt, dt)
u = 0
while (t < tmax):
    ccfourier = cahnfourier1d(ccfourier, k2, k4)
    ccexp = cahnexp1d(ccexp, c1, c2)
    if u%1000 == 0:
        diferenca = abs(ccexp-ccfourier)
        diferencas = [diferenca[0], diferenca[42], diferenca[75], diferenca[103]]
        dif1.append(diferencas[0])
        dif2.append(diferencas[1])
        dif3.append(diferencas[2])
        dif4.append(diferencas[3])
    t = round(t + dt, int(-np.log10(dt) + 5))
    u+=1
tt = np.linspace(dt, tmax, len(dif1))
plt.title(f"Tempo final: {t:.7f}")
plt.plot(tt, dif1)
plt.plot(tt, dif2)
plt.plot(tt, dif3)
plt.plot(tt, dif4)
plt.yscale("log")
plt.savefig(f".\\1d\\aaaaa{t:.6f}.png", dpi = tamanho)