ChatGPTに拡散方程式のプログラムを書いてと頼んだ。

ChatGPT使ってみた

プログラムを書いてくれるというので拡散方程式の差分のFortranプログラムを書いてもらいました。

結論からいうと、合ってそうで正しくないコードが出てきました。質問を日本語で適当にしたからなのか、陰解法でとかめんどくさいことをいったからなのか。

最初のお願い

拡散方程式の空間を2次精度の中心差分，時間を1次精度の陰的Euler法で離散化するFortranのプログラム書いて

ちなみに想定される離散化した式が以下。

\frac{u^{n+1}_i - u^n_i}{\Delta t} = \alpha \frac{u^{n+1}_{i+1} - 2u^{n+1}_i + u^{n+1}_{i-1}}{\Delta x^2}

上の添字が時刻，下の添字が空間座標のインデックス。既知の$u^n$から次の時刻の$u^{n+1}$を計算します。行列でかくと

\begin{pmatrix}
  1+2\nu  & -\nu      & 0 \dots &       &        0                         \\
   -\nu      & 1+2\nu          &  -\nu      & 0 \dots &0 \\
         &     \ddots            & \ddots            & &         \\
        &-\nu &1+2\nu & -\nu & \\
       &  &  \ddots      & \ddots      &      \\
   0 \dots      & 0 & -\nu& 1+2\nu & -\nu \\
   0  \dots    &                &    0    &     -\nu      & 1+2\nu
\end{pmatrix}
\begin{pmatrix}
u^{n+1}_1\\
u^{n+1}_2 \\
\vdots\\
u^{n+1}_i\\
\vdots\\
u^{n+1}_{imax-1} \\
u^{n+1}_{imax}
\end{pmatrix} = 
\begin{pmatrix}
u^{n}_1\\
u^n_2\\
\vdots\\
u^{n}_i\\
\vdots\\
u^{n}_{imax-1} \\
u^{n}_{imax}
\end{pmatrix}

ここで$\nu=\alpha \Delta t/\Delta x^2$。この連立1次方程式を解いてくれればOKです。

回答

以下に、拡散方程式の空間を2次精度の中心差分、時間を1次精度の陰的Euler法で離散化するFortranのプログラムを示します。このプログラムは、1次元の場合を想定しています。

program diffusion
  implicit none
  integer, parameter :: nx = 100 ! 空間の分割数
  integer, parameter :: nt = 100 ! 時間ステップ数
  real, parameter :: dx = 0.01 ! 空間刻み幅
  real, parameter :: dt = 0.01 ! 時間刻み幅
  real, parameter :: alpha = 0.1 ! 拡散係数
  real, dimension(0:nx) :: u ! 解の配列
  integer :: i, j ! インデックス用変数
  
  ! 初期条件
  do i = 0, nx
    u(i) = 0.0
    if (i >= 45 .and. i <= 55) u(i) = 1.0
  end do
  
  ! 時間発展
  do j = 1, nt
    do i = 1, nx-1
      u(i) = u(i) + alpha * dt / dx**2 * (u(i+1) - 2*u(i) + u(i-1))
    end do
  end do
  
  ! 結果の出力
  do i = 0, nx
    write(*,*) i*dx, u(i)
  end do
end program diffusion

陰解法を頼んだのに陽解法。。。。よく見ると

u(i) = u(i) + alpha * dt / dx**2 * (u(i+1) - 2*u(i) + u(i-1))

だと，u(i-1)が次のタイムステップの量になって，計算がおかしくなる。

次

これは陽的Euler法のプログラムです。陰解法のプログラムを書いて下さい。

回答

陰的Euler法の場合、更新式が差分方程式を解く必要があります。一般的には、連立方程式を解く必要がありますが、今回は簡単のために反復法を用いたガウス・ザイデル法を使って、差分方程式を解きます。以下に、陰的Euler法による拡散方程式のFortranプログラムを示します。

program diffusion
  implicit none
  integer, parameter :: nx = 100 ! 空間の分割数
  integer, parameter :: nt = 100 ! 時間ステップ数
  real, parameter :: dx = 0.01 ! 空間刻み幅
  real, parameter :: dt = 0.01 ! 時間刻み幅
  real, parameter :: alpha = 0.1 ! 拡散係数
  real, dimension(0:nx) :: u ! 解の配列
  real, dimension(0:nx) :: u_new ! 新しい解の配列
  integer :: i, j, k ! インデックス用変数
  real :: error, tol ! 反復の収束判定用変数
  
  ! 初期条件
  do i = 0, nx
    u(i) = 0.0
    u_new(i) = 0.0
    if (i >= 45 .and. i <= 55) u(i) = 1.0
  end do
  
  ! 時間発展
  do j = 1, nt
    do i = 1, nx-1
      ! 初期値の設定
      u_new(i) = u(i)
      error = 1.0
      tol = 1.0e-6
      
      ! 反復法による解法
      do k = 1, 1000
        u_new(i) = (u(i+1) + u_new(i-1) + alpha * dt / dx**2 * u(i-1) + alpha * dt / dx**2 * u(i+1)) / (2.0*(1.0 + alpha*dt/dx**2))
        error = max(error, abs(u_new(i) - u(i)))
        u(i) = u_new(i)
        if (error < tol) exit
      end do
    end do
  end do
  
  ! 結果の出力
  do i = 0, nx
    write(*,*) i*dx, u(i)
  end do
end program diffusion

ガウス・ザイデル法ってどんなんだっけと思いながら，教科書見返したら，パット見あってそうだけど，よく見るとこれもヘン。ガウス・ザイデル法で離散化すると

(1+\frac{2\alpha \Delta t}{\Delta x^2})u^{(k+1)}_i = \frac{\alpha \Delta t}{\Delta x^2} u^{(k+1)}_{i-1} + \frac{\alpha \Delta t}{\Delta x^2} u^{(k)}_{i+1} +u^n_i

kがガウス・ザイデルの反復の回数。

u_new(i) = (u(i+1) + u_new(i-1) + alpha * dt / dx**2 * u(i-1) + alpha * dt / dx**2 * u(i+1)) / (2.0*(1.0 + alpha*dt/dx**2))

余計な項がある。これをもとにして書き直したコードがこれ。

program diffusion
  implicit none
  integer, parameter :: nx = 100 ! 空間の分割数
  integer, parameter :: nt = 30 ! 時間ステップ数
  real, parameter :: dx = 0.01 ! 空間刻み幅
  real, parameter :: dt = 0.01 ! 時間刻み幅
  real, parameter :: alpha = 0.1 ! 拡散係数
  real, dimension(0:nx) :: u ! 解の配列
  real, dimension(0:nx) :: u_new ! 新しい解の配列
  real, dimension(0:nx) :: u_pre ! 
  integer :: i, j, k ! インデックス用変数
  real :: error, tol ! 反復の収束判定用変数
  real :: nu
  
  ! 初期条件
  do i = 0, nx
     u(i) = 0.0
     u_new(i) = 0.0
     if (i >= 45 .and. i <= 55) u(i) = 1.0
  end do

  tol = 1.0e-6
  nu = alpha * dt / dx**2
  
  ! 時間発展
  do j = 1, nt
     do i = 1, nx-1
        ! 初期値の設定
        u_new(i) = u(i)
        u_pre(i) = u(i)
     enddo

     ! 反復法による解法
     do k = 1, 500
        error = -1.0
        do i = 1, nx-1
           
           u_new(i) = (u_new(i-1) *nu + u(i+1) *nu + u_pre(i))/(1.0+2.0*nu) 
           error = max(error, abs(u_new(i) - u(i)))
           u(i) = u_new(i)
           
        end do
        write(*,*) k, error
        if (error < tol) exit

     end do

     ! 結果の出力
     write(*,*) j
    
     do i = 0, nx
        write(10*j,*) i*dx, u(i)
     end do
         
  end do
end program diffusion

現状

流体とか差分のプログラムを書いてもらうとそれっぽいプログラムを書いてくれるが，正しいかを割とまじめにチェックする必要がありそう。

補足

英語でお願いしてみました。

Make a Fortran program for solving the one-dimensional diffusion equation. Please use Euler implicit method for time marching algorithm and second accuracy central differencing for spatial discretization. Please solve the linear algebra equation by Gauss-Sidel method.

できたコード

program diffusion
    implicit none

    ! Define parameters
    real(kind=8) :: L = 1.0        ! Length of domain
    integer :: Nx = 100           ! Number of spatial grid points
    real(kind=8) :: T = 1.0       ! Simulation time
    integer :: Nt = 1000          ! Number of time steps
    real(kind=8) :: D = 1.0       ! Diffusion coefficient
    real(kind=8) :: dx, dt       ! Spatial grid spacing, time step size
    real(kind=8), dimension(:), allocatable :: u, unew ! Solution arrays
    real(kind=8), dimension(:,:), allocatable :: A   ! Matrix for implicit method
    real(kind=8) :: r             ! Implicit method coefficient
    integer :: i, n, j            ! Loop variables

    ! Allocate arrays
    allocate(u(Nx), unew(Nx), A(Nx,Nx))

    ! Calculate grid spacing and time step size
    dx = L / real(Nx-1)
    dt = T / real(Nt)

    ! Set initial condition
    do i = 2, Nx-1
        u(i) = sin(real(i-1) * pi * dx)
    end do
    u(1) = 0.0
    u(Nx) = 0.0

    ! Set up matrix for implicit method
    r = D * dt / dx**2
    do i = 1, Nx
        do j = 1, Nx
            if (i == j) then
                A(i,j) = 1.0 + 2.0*r
            else if (i == j+1 .or. i == j-1) then
                A(i,j) = -r
            else
                A(i,j) = 0.0
            end if
        end do
    end do

    ! Time stepping loop
    do n = 1, Nt
        ! Solve linear algebra equation by Gauss-Sidel method
        do i = 1, Nx
            unew(i) = u(i)
        end do
        do i = 2, Nx-1
            u(i) = (unew(i-1) + unew(i+1) + r*u(i-1) + r*u(i+1)) / (1.0 + 2.0*r)
        end do
        u(1) = 0.0
        u(Nx) = 0.0
    end do

    ! Deallocate arrays
    deallocate(u, unew, A)

end program diffusion

係数行列を準備するとかいいじゃないかとおもってみたら，係数行列Aは定義するだけしておいて，使わない。依然として，配列の計算がおかしい。

ちなみにJuliaで作ってっていったら，係数行列をそのまま反転するというプログラムを作ってくれました。

Make a Julia program for solving the one-dimensional diffusion equation. Please use Euler implicit method for time marching algorithm and second accuracy central differencing for spatial discretization.

# Define parameters
L = 1.0        # Length of domain
Nx = 100       # Number of spatial grid points
T = 1.0        # Simulation time
Nt = 1000      # Number of time steps
D = 1.0        # Diffusion coefficient
dx = L / (Nx-1)
dt = T / Nt

# Create arrays for solution
u = zeros(Nx)
u[1] = 0.0
u[Nx] = 0.0
unew = zeros(Nx)

# Set initial condition
for i in 2:Nx-1
    x = (i-1) * dx
    u[i] = sin(pi*x)
end

# Set up matrix for implicit method
r = D * dt / dx^2
A = zeros(Nx,Nx)
A[1,1] = 1.0
A[Nx,Nx] = 1.0
for i in 2:Nx-1
    A[i,i] = 1 + 2*r
    A[i,i-1] = -r
    A[i,i+1] = -r
end

# Time stepping loop
for n in 1:Nt
    b = u
    unew = A \ b
    u = unew
end

# Plot solution
using Plots
x = range(0, L, length=Nx)
plot(x, u, xlabel="x", ylabel="u(x)", label="Numerical solution")
plot!(x, sin.(pi*x), label="Analytical solution")