Published on

1次元の拡散方程式シミュレーション

Authors

1次元の拡散方程式

差分近似方程式の特解

ujk+1=ujk+A[uj+1k+uj1k2ujk]u_{j}^{k+1}=u_{j}^{k}+A\left[u_{j+1}^{k}+u_{j-1}^{k}-2 u_{j}^{k}\right]

この式の特解を ujk=(g)keiqxju_{j}^{k}=(g)^{k} e^{i q x_{j}} とおくと (g=ujk+1/ujk)\left(g=u_{j}^{k+1} / u_{j}^{k}\right)

gujk=[1+A(eiqΔx+eiqΔx2)]ujkg=1+2A(cosqΔx1)\begin{array}{r} g u_{j}^{k} &=& \left[1+A\left(e^{i q \Delta x}+e^{-i q \Delta x}-2\right)\right] u_{j}^{k} \\ g &=& 1+2 A(\cos q \Delta x-1) \end{array}

解が有限にとどまる条件: g1|g| \leq 1

1cosqΔx1-1 \leq \cos q \Delta x \leq 1より A12A \leq \frac{1}{2}が安定の条件.

DΔtΔx212ΔtΔx22D\frac{D \Delta t}{\Delta x^{2}} \leq \frac{1}{2} \rightarrow \Delta t \leq \frac{\Delta x^{2}}{2 D}

問題設定

1次元の拡散方程式:

u(x,t)t=D2u(x,t)x2(0x1,t0)\frac{\partial u(x, t)}{\partial t}=D \frac{\partial^{2} u(x, t)}{\partial x^{2}} \quad(0 \leq x \leq 1, \quad t \geq 0)

境界条件:

u(8t)=u(1,t)=0u(8 t)=u(1, t)=0

初期条件:

u(x,0)={2x(0x12)2(1x)(12x1)u(x, 0)=\left\{\begin{aligned} 2 x &\left(0 \leq x \leq \frac{1}{2}\right) \\ 2(1-x) &\left(\frac{1}{2} \leq x \leq 1\right) \end{aligned}\right.

計算条件:

  • 前進オイラー+2 次精度中心差分
    uik+1=uik+A[ui+1k+ui1k2uik],A=DΔtΔx2u_{i}^{k+1}=u_{i}^{k}+A\left[u_{i+1}^{k}+u_{i-1}^{k}-2 u_{i}^{k}\right], \quad A=\frac{D \Delta t}{\Delta x^{2}} xi=iΔx(i=0,1,2,N)tk=kΔt(k=0,1,2,,M)x_{i}=i \Delta x(i=0,1,2 \cdots, N) \quad t_{k}=k \Delta t(k=0,1,2, \cdots, M)

  • 計算パラメータ
    D=1,N=10,M=100,Δx=0.1,Δt=0.002D=1, N=10, M=100, \Delta x=0.1, \Delta t=0.002

プログラム

fortran で実装した

            integer::i, k, tstp
            integer, parameter::M=100,N=10
            real*8::u(0:N),utmp(0:N),uout(0:N,0:10)
            real*8,parameter::L=1.d0,D=1.d0,dt=2.d-3
            real*8,parameter::dx=L/dble(N),A=D*dt/dx/dx

            tstp=0

            do i=0, N/2
                u(i)=2*i*dx
            enddo
            do i=N/2+1,N
                u(i)=2*(1.d0-i*dx)
            enddo

            do k=0,M
                if (mod(k,20)==0) then
                    do i=0,N
                        uout(i,tstp)=u(i)
                    enddo
                    tstp=tstp+1
                endif

                u(0)=0d0
                u(N)=0d0
                do i=1,N-1
                    utmp(i)=u(i)+A*(u(i+1)+u(i-1)-2*u(i))
                enddo

                do i=1,N-1
                    u(i)=utmp(i)
                enddo
            enddo

            open(10,file='u.dat', status='unknown')
            do i=0,N
                write(10,'(10f9.5)')i*dx,(uout(i,k),k=0,tstp-1)
            enddo
            close(10)

            end

実行結果

null

安定していることが分かる