Diffusion_2D_FTCS_v0.6.jl

# x
Δx  = 0.1
X   = 10.0
xlm = [0, X]
x   = collect(xlm[1]:Δx:xlm[2])
I   = length(x)

# y
Δy  = 0.1
Y   = 10.0
ylm = [0, Y]
y   = collect(ylm[1]:Δy:ylm[2])
J   = length(y)

# t
Δt  = 0.0025
T   = 2.0
tlm = [0,T]
t   = collect(tlm[1]:Δt:tlm[2])
N   = length(t)

# 初期化
u = zeros(I, J, N)

# 境界条件
u[  1,   :, :] = 0.0
u[end,   :, :] = 0.0
u[  :,   1, :] = 0.0
u[  :, end, :] = 0.0

# 初期条件
xx = repmat(x',J,1)
yy = repmat(y,1,I)
φ(x,y) = exp( - ( ( x - X / 2.0 )^2 + ( y - Y / 2.0 )^2) )
u[:,:,1] = φ.(xx, yy)

# 数値計算
for n = 1:N-1
    for i = 2:I-1
        for j = 2:J-1
            u[i,j,n+1] = ( 1.0 - 2.0 * Δt * ( 1.0 / Δx^2 + 1.0 / Δy^2 ) ) * u[i,j,n] + 
                         Δt * ( u[i-1,j,n] + u[i+1,j,n] ) / Δx^2 + 
                         Δt * ( u[i,j-1,n] + u[i,j+1,n] ) / Δy^2
        end
    end
end

# CFL = Δt*(1/Δx**2 + 1/Δy**2)
# print(CFL)
# plot
#----------------------------------------------
using PyPlot
using PyCall
@pyimport matplotlib.animation as anim

# figure
fig = figure(figsize=(8, 4))
subplot(1,2,1)
# image
image = imshow(u[:,:,1],cmap="hot",interpolation="none",origin="lower",extent=[xlm[1],xlm[2],ylm[1],ylm[2]])
# image = pcolormesh(xx,yy,u[:,:,0])
axhline(y[Int(round(J/2))],ls="--",color="C0",lw=1)
xlabel("x")
ylabel("y")
xticks(xlm[1]:2.5:xlm[2])
yticks(ylm[1]:2.5:ylm[2])
tight_layout()

subplot(1,2,2)
# line
line = plot([],[])[1]

# setting
xlim(xlm)
ylim([0,1])
grid()
xlabel("x")
ylabel("u")
xticks(xlm[1]:2.5:xlm[2])
grid("on",linestyle="--")
title(@sprintf("%4.2f",t[1]))
tight_layout()

interval = 5
# update
function update(i)
    i = i * interval
    line[:set_data](x, u[:,Int(round(J/2)),i+1])
    image[:set_data](u[:,:,i+1])
    # image.set_array(u[:-1,:-1,i].flatten())
    # title
    title(@sprintf("%4.2f",t[i+1]))
end
# frame number
frames = Int(round(N/interval))
# animation
ani = anim.FuncAnimation(fig, update, frames=frames, interval=10)
# save
ani[:save]("Diffusion_2D_FTCS.gif", writer="imagemagick")

Diffusion_2D_FTCS_v1.jl

# x
Δx  = 0.1
X   = 10.0
xlm = [0, X]
x   = collect(xlm[1]:Δx:xlm[2])
I   = length(x)

# y
Δy  = 0.1
Y   = 10.0
ylm = [0, Y]
y   = collect(ylm[1]:Δy:ylm[2])
J   = length(y)

# t
Δt  = 0.0025
T   = 2.0
tlm = [0,T]
t   = collect(tlm[1]:Δt:tlm[2])
N   = length(t)

# 初期化
u = zeros(I, J, N)

# 境界条件
u[  1,   :, :] .= 0.0
u[end,   :, :] .= 0.0
u[  :,   1, :] .= 0.0
u[  :, end, :] .= 0.0

# 初期条件
xx = repeat(x',J,1)
yy = repeat(y,1,I)
φ(x,y) = exp( - ( ( x - X / 2.0 )^2 + ( y - Y / 2.0 )^2) )
u[:,:,1] = φ.(xx, yy)

# 数値計算
for n = 1:N-1
    for i = 2:I-1
        for j = 2:J-1
            u[i,j,n+1] = ( 1.0 - 2.0 * Δt * ( 1.0 / Δx^2 + 1.0 / Δy^2 ) ) * u[i,j,n] + 
                         Δt * ( u[i-1,j,n] + u[i+1,j,n] ) / Δx^2 + 
                         Δt * ( u[i,j-1,n] + u[i,j+1,n] ) / Δy^2
        end
    end
end

# CFL = Δt*(1/Δx**2 + 1/Δy**2)
# print(CFL)
# plot
#----------------------------------------------
using PyPlot
using PyCall
@pyimport matplotlib.animation as anim

# @spritnf
using Printf

# figure
fig = figure(figsize=(8, 4))
subplot(1,2,1)
# image
image = imshow(u[:,:,1],cmap="hot",interpolation="none",origin="lower",extent=[xlm[1],xlm[2],ylm[1],ylm[2]])
axhline(y[Int(round(J/2))],ls="--",color="C0",lw=1)
xlabel("x")
ylabel("y")
xticks(xlm[1]:2.5:xlm[2])
yticks(ylm[1]:2.5:ylm[2])
tight_layout()

subplot(1,2,2)
# line
line = plot([],[])[1]

# setting
xlim(xlm)
ylim([0,1])
grid()
xlabel("x")
ylabel("u")
xticks(xlm[1]:2.5:xlm[2])
grid("on",linestyle="--")
title(@sprintf("%4.2f",t[1]))
tight_layout()

interval = 5
# update
function update(i)
    i = i * interval
    line[:set_data](x, u[:,Int(round(J/2)),i+1])
    image[:set_data](u[:,:,i+1])
    # title
    title(@sprintf("%4.2f",t[i+1]))
end
# frame number
frames = Int(round(N/interval))
# animation
myanim = anim.FuncAnimation(fig, update, frames=frames, interval=10)
# save
myanim[:save]("Diffusion_2D_FTCS_v1.gif", writer="imagemagick")

gistfile1.txt

# x
Δx  = 0.1
X   = 10.0
xlm = [0, X]
x   = collect(xlm[1]:Δx:xlm[2])
I   = length(x)

# y
Δy  = 0.1
Y   = 10.0
ylm = [0, Y]
y   = collect(ylm[1]:Δy:ylm[2])
J   = length(y)

# t
Δt  = 0.0025
T   = 2.0
tlm = [0,T]
t   = collect(tlm[1]:Δt:tlm[2])
N   = length(t)

# 初期化
u = zeros(I, J, N)

# 境界条件
u[  1,   :, :] = 0.0
u[end,   :, :] = 0.0
u[  :,   1, :] = 0.0
u[  :, end, :] = 0.0
# 初期条件
xx = repmat(x',J,1)
yy = repmat(y,1,I)
φ(x,y) = exp( - ( ( x - X / 2.0 )^2 + ( y - Y / 2.0 )^2) )
u[:,:,1] = φ.(xx, yy)

# 数値計算
for n = 1:N-1
    for i = 2:I-1
        for j = 2:J-1
            u[i,j,n+1] = ( 1.0 - 2.0 * Δt * ( 1.0 / Δx^2 + 1.0 / Δy^2 ) ) * u[i,j,n] + 
                         Δt * ( u[i-1,j,n] + u[i+1,j,n] ) / Δx^2 + 
                         Δt * ( u[i,j-1,n] + u[i,j+1,n] ) / Δy^2
        end
    end
end

# CFL = Δt*(1/Δx**2 + 1/Δy**2)
# print(CFL)
# plot
#----------------------------------------------
using PyPlot
using PyCall
@pyimport matplotlib.animation as anim

# figure
fig = figure(figsize=(8, 4))
subplot(1,2,1)
# image
image = imshow(u[:,:,1],cmap="hot",interpolation="none",origin="lower",extent=[xlm[1],xlm[2],ylm[1],ylm[2]])
# image = pcolormesh(xx,yy,u[:,:,0])
axhline(y[Int(round(J/2))],ls="--",color="C0",lw=1)
xlabel("x")
ylabel("y")
xticks(xlm[1]:2.5:xlm[2])
yticks(ylm[1]:2.5:ylm[2])
tight_layout()

subplot(1,2,2)
# line
line = plot([],[])[1]

# setting
xlim(xlm)
ylim([0,1])
grid()
xlabel("x")
ylabel("u")
xticks(xlm[1]:2.5:xlm[2])
grid("on",linestyle="--")
title(@sprintf("%4.2f",t[1]))
tight_layout()

interval = 5
# update
function update(i)
    i = i * interval
    line[:set_data](x, u[:,Int(round(J/2)),i+1])
    image[:set_data](u[:,:,i+1])
    # image.set_array(u[:-1,:-1,i].flatten())
    # title
    title(@sprintf("%4.2f",t[i+1]))
end
# frame number
frames = Int(round(N/interval))
# animation
ani = anim.FuncAnimation(fig, update, frames=frames, interval=10)
# save
ani[:save]("Diffusion_2D_FTCS.gif", writer="imagemagick")