Cahn1dcompare: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Criou página com '<source lang = "python"> 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(...')
 
Sem resumo de edição
 
Linha 22: Linha 22:
     cc = irfft(cct)
     cc = irfft(cct)
     return cc
     return cc
 
#SEMENTE UTILIZADA NA FIGURA DA PÁGINA: 4
np.random.seed(4)
np.random.seed(4)



Edição atual tal como às 18h51min de 29 de setembro de 2022

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
#SEMENTE UTILIZADA NA FIGURA DA PÁGINA: 4
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