Cahn1dcompare

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 = 1

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)

xx = np.linspace(0, xmax, l)

t = 0

u = 0

tamanho = l/2
plt.figure(figsize=(8, 8))
os.makedirs(f".\\1d", exist_ok = True)
while (t < tmax):
    plt.cla()
    plt.title(f"Tempo: {t:.7f}, maior erro: {np.amax(abs(ccexp-ccfourier)):.7f}")
    plt.axhline(0, c = "red", alpha = 0.5)
    plt.ylim(-1.1, 1.1)
    plt.xlim(xmin, xmax)
    plt.scatter(xx, ccexp, c = "orange", alpha = 0.6, label = "FTCS")
    plt.plot(xx, ccexp, c = "orange", alpha = 0.4)
    plt.scatter(xx, ccfourier, c = "blue", alpha = 0.4, label = "FFT")
    plt.plot(xx, ccfourier, c = "blue", alpha = 0.6)
    plt.scatter(xx, abs(ccexp-ccfourier), c = "red", alpha = 0.4, label = "Erro")
    plt.legend(loc = 1)
    if u%100 == 0:
        plt.savefig(f".\\1d\\{t:.6f}.png", dpi = tamanho)
    if u ==10000:
        break
    ccfourier = cahnfourier1d(ccfourier, k2, k4)
    ccexp = cahnexp1d(ccexp, c1, c2)
    t = round(t + dt, int(-np.log10(dt) + 1))
    u+=1