Google ColabのJulia言語で搭載されているGeminiを使って一行もコードを書かずに2次元拡散方程式を差分法で計算してGIFアニメにする。次に同じように2次元波動方程式もやってもらう。

今回は生成AI、Geminiを使ってみよう。

最初の指示は「Julia言語で2次元拡散方程式を差分法で計算するコードを書いてください。」

とした。最初からまあまあいい感じだが境界条件がない。

ノイマン型で値は0として、とするとちゃんと付け加えた。

今度は解が不安定でNaNになる。ここも指摘すると、

dt = dx^2 * dy^2 / (2 * D * (dx^2 + dy^2)) * 0.5 # 安定性条件を満たすようにdtを計算、0.5倍することで余裕を持たせる

という安定性条件をちゃんと出してきた。

最終的なコードは

using Plots

using Printf

# パラメータ設定

nx = 100 # x方向の格子点数

ny = 100 # y方向の格子点数

dx = 1.0 / nx # x方向の格子間隔

dy = 1.0 / ny # y方向の格子間隔

# 安定性条件を満たすようにdtを計算

D = 1.0 # 拡散係数

dt = dx^2 * dy^2 / (2 * D * (dx^2 + dy^2)) * 0.5 # 0.5倍することで余裕を持たせる

nsteps = 5000 # 計算ステップ数

# 初期条件設定

u = zeros(nx, ny)

for i in 1:nx

    for j in 1:ny

        if (i - nx / 2)^2 + (j - ny / 2)^2 < (nx / 4)^2

            u[i, j] = 1.0

        end

    end

end

# 結果を格納する配列

results = []

# 時間発展計算

for t in 1:nsteps

    u_new = copy(u)

    for i in 2:nx-1

        for j in 2:ny-1

            u_new[i, j] = u[i, j] + D * dt * (

                (u[i+1, j] - 2u[i, j] + u[i-1, j]) / dx^2 +

                (u[i, j+1] - 2u[i, j] + u[i, j-1]) / dy^2

            )

        end

    end

    # ノイマン境界条件を適用

    u_new[1, :] = u_new[2, :] # x = 0

    u_new[nx, :] = u_new[nx-1, :] # x = 1

    u_new[:, 1] = u_new[:, 2] # y = 0

    u_new[:, ny] = u_new[:, ny-1] # y = 1

    u = u_new

    # 100ステップごとに結果を配列に格納

    if t % 100 == 0

        push!(results, copy(u))

    end

end

# 計算結果をアニメーションで表示、heatmapのスケールを0から1に固定

anim = @animate for i in 1:length(results)

    heatmap(results[i], title="t = $(@sprintf("%.3f", i * 100 * dt))", clim=(0, 1))  

end

gif(anim, "diffusion_neumann.gif", fps = 10)

実行結果はこちら。ちゃんと意図したようなグラフができた。

次は波動方程式。

using Plots

using Printf

# パラメータ設定

nx = 200 # x方向の格子点数

ny = 200 # y方向の格子点数

dx = 1.0 / nx # x方向の格子間隔

dy = 1.0 / ny # y方向の格子間隔

c = 1.0 # 波の速度

dt = minimum([dx, dy]) / c * 0.5 # 安定性条件を満たすようにdtを計算 (0.5倍することで余裕を持たせる)

nsteps = 500 # 計算ステップ数

# 初期条件設定

u = zeros(nx, ny)

u_old = zeros(nx, ny)

for i in 1:nx

    for j in 1:ny

        if (i - nx / 2)^2 + (j - ny / 2)^2 < (nx / 16)^2

            u[i, j] = 1.0

            u_old[i, j] = 1.0

        end

    end

end

# 結果を格納する配列

results = []

# 時間発展計算

for t in 1:nsteps

    u_new = zeros(nx, ny)

    for i in 2:nx-1

        for j in 2:ny-1

            u_new[i, j] = 2u[i, j] - u_old[i, j] + c^2 * dt^2 * (

                (u[i+1, j] - 2u[i, j] + u[i-1, j]) / dx^2 +

                (u[i, j+1] - 2u[i, j] + u[i, j-1]) / dy^2

            )

        end

    end

    # ノイマン境界条件を適用 (値0)

    u_new[1, :] = u_new[2, :]  # x = 0

    u_new[nx, :] = u_new[nx-1, :] # x = 1

    u_new[:, 1] = u_new[:, 2]  # y = 0

    u_new[:, ny] = u_new[:, ny-1] # y = 1

    u_old = u

    u = u_new

    # 10ステップごとに結果を配列に格納

    if t % 10 == 0

        push!(results, copy(u))

    end

end

# 計算結果をアニメーションで表示、heatmapのスケールを-1から1に固定

anim = @animate for i in 1:length(results)

    heatmap(results[i], title="t = $(@sprintf("%.3f", i  *10 * dt))", clim=(-1, 1))

end

gif(anim, "wave2d_neumann.gif", fps = 10)

これも予期した動きをしている。