---
robots: noindex, nofollow
lang: ja
dir: ltr
breaks: true
---
$\newcommand{\bv}[1]{\boldsymbol{ #1 }}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\grad}[0]{\text{Grad}}
\newcommand{\w}[0]{\hat{\mathcal W}}
\DeclareMathOperator{\det}{det}
\DeclareMathOperator{\trace}{tr}
\DeclareMathOperator{\lim}{lim}$

# 反応拡散方程式を差分法(陽的解法)で解く
JuliaかC++で1次元での反応拡散方程式を解いていきます.差分法で解いていきますが,陰的解法だと非線形の連立方程式を解いていく必要があり,それは面倒なので,陽的解法で解いていきます.

## 1.反応拡散方程式
1次元の反応拡散方程式は以下のようになる.
\begin{align}
\frac{\partial u}{\partial t}=\epsilon^2 \frac{\partial^2 u}{\partial x^2}+u(1-u^2)
\end{align}
また,今回の解析領域は$x=[0,2\pi]$とする.

## 2.差分近似
差分近似について述べる.2変数関数$u(x,t)$に関して，テイラー展開を考える．
まず，$x$に関してテイラー展開すると
\begin{align}
u(x+\Delta x,t) = u(x,t) +\Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} + \cdots
\end{align}
\begin{align}
u(x-\Delta x,t) = u(x,t) -\Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2!} \frac{\partial^2 u}{\partial x^2} + \cdots
\end{align}

となる．上2式の辺々を加えることにより,
\begin{align}
\frac{\partial^2 u}{\partial x^2} = \frac{u(x+\Delta x,t) -2 u(x,t) + u(x - \Delta x,t)}{\Delta x^2}
\end{align}
が得られる.

## 3.陽的解法

$x$方向の刻み幅を$h$，$t$方向の刻み幅を$k$とする.　$i$を$x$方向のステップ数, $j$を$t$方向のステップ数とし,以降は, $u(x,t)=u_{i, j}$と表す．このとき，$\frac{\partial^2 u}{\partial x^2}$は中心差分を用いると，以下のように差分化される．
\begin{align}
\frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1, j} -2 u_{i, j} + u_{i-1, j}}{h^2}
\end{align}

また$\frac{\partial u}{\partial t}$は前進差分により
\begin{align}
\frac{\partial u}{\partial t}=\frac{u_{i, j+1}-u_{i, j}}{k}
\end{align}
となる.

これらを1節の反応拡散方程式に代入すると以下のようになる.
\begin{align}
\frac{u_{i, j+1}-u_{i, j}}{k}=\epsilon^2 \frac{u_{i+1, j}-2u_{i, j}+u_{i-1, j}}{h^2}+u_{i, j}\left(1-u_{i, j}^2\right)
\end{align}
これを$u_{i, j+1}$について解くと,
\begin{align}
u_{i, j+1}=\frac{k \epsilon^2}{h^2}u_{i+1, j}-ku_{i, j}^3+\left(1+k-\frac{2k\epsilon^2}{h^2}\right)u_{i, j}+\frac{k\epsilon^2}{h^2}u_{i-1, j}
\end{align}

## 4.初期条件と境界条件
### 4.1 初期条件($t=0$)
初期条件は以下のようにする.
\begin{align}
u_{i, 0}=\cos(x)=\cos(ih)
\end{align}
### 4.2 境界条件($x=0, 2\pi$)

<!-- Newman境界条件として,以下の条件を課す.
\begin{align}
\left.\frac{\partial u}{\partial x}\right|_{x=0}=\left.\frac{\partial u}{\partial x}\right|_{x=1}=0
\end{align}

この式に差分を適用すると以下のようになる.
\begin{align}
\frac{u_{1, j}-u_{0, j}}{h}=\frac{u_{N, j}-u_{N-1, j}}{h}=0
\end{align}

よって,以下の条件が導かれる.
\begin{align}
u_{1, j}=u_{0, j},\hspace{10mm}u_{N, j}=u_{N-1, j}
\end{align} -->
Dirichlet境界条件として以下を課す.
\begin{align}
&u_{0, j}=1\\
&u_{N, j}=1
\end{align}
<!-- ## 5.計算結果(例) -->
#### サンプルプログラム

using Plots
using FileIO

#位置xは0から2πまで
n=100
x=zeros(Float64,n)
h=2π/(n-1)
for i=1:n-1
    x[i+1]=i*h
end
#時刻は0から200000まで(十分大きく、反応が拡散しきった定常状態まで至ればOK)
m=2000
t=zeros(Float64,m)
k=0.01
#t[1]は0
for j=1:m-1
    t[j+1]=t[j]+k
end

#濃度u
u=zeros(Float64,(n,m))

#拡散係数
ϵ=0.01
#初期条件
for i=1:n
    u[i,1]=cos((i-1)*h)#h=2π/(n-1)なのでi=nでちょうどcos2πになる
end

#境界条件
u[1,:]=ones(Float64,m)
u[n,:]=ones(Float64,m)
for j=1:m-1
    for i=2:n-1
        u[i,j+1]=k*ϵ^2/h^2*u[i+1,j]-k*u[i,j]^3+(1+k-2*k*ϵ^2/h^2)u[i,j]+k*ϵ^2/h^2*u[i-1,j]
    end
end
#@animate for 時刻plot()で時刻ごとのグラフをつなげてアニメにできる
#plot(x,f(x))でf(x)のグラフが書ける。オプションで範囲指定できる。leg=falseで凡#例を消す。
anim=@animate for j=1:m
    plot(x,u[:,j],ylim=(-1.0,1.0),leg=false)
end
gif(anim,"hannoukakusan.gif",fps=50)