next up previous contents index
: 不安定現象 : 常微分方程式の数値解 : 常微分方程式の数値解   目次   索引

差分法

$f'(t)$$h$ が十分小さいとき,

\begin{displaymath}\frac{f(t+h)-f(t)}{h}\end{displaymath}

で近似できる. 同じように $f''(t)$$h$ が十分小さいとき,
\begin{displaymath}
\frac{f(t+h)-2f(t)+f(t-h)}{h^2}
\end{displaymath} (9.1)

で近似できることをたとえば Taylor 展開を用いて証明できる. 微分をこのように近似することにより, 微分方程式の近似解を, 数列の漸化式を解くことにより もとめることが可能である.

たとえば, 単振動の方程式

\begin{displaymath}{{d^2} \over {dt^2}} y + y = 0\end{displaymath}

を数値的に解くことを考えてみよう. 解く前に, この方程式の物理的意味を復習しておこう. 古典物理で習うように,
質量 × 加速度 = 力
なる関係式がなりたつ. この関係式を Newton の運動方程式という. 物体の時刻 $t$ における位置を $q(t)$ とおくと, 速度は $q'(t)$, 加速度は $q''(t)$ である.

1 次元的に単振動をするバネについた質量 $1$ の物体 W を考えよう. 時刻 $t$ における, $W$ の位置を $y(t)$ とすることにする. ただし, バネが自然な長さにあるとき $y=0$ とする. フックの法則によると, 自然な長さから $y$ だけのびた(負のときはちじんだとみなす)とき物体 $W$ にかかる力は $-k y $ である. ここで $k$ はバネできまる定数. ここで $k=1$ と仮定して Newton の運動方程式を適用すると, 単振動の方程式

\begin{displaymath}y'' + y = 0\end{displaymath}

を得る.

$A$, $B$ を定数として一般解は $A \cos(t) + B \sin(t)$ であるが, この方程式を数値解法で解こう. 数値解法の利点は, $\cos$$\sin$ で解を書けないときでも, 微分方程式の近似解がわかることである. 式 (9.1) より, $h$ が十分小さいとき,

\begin{displaymath}\frac{y(t+h)-2y(t)+y(t-h)}{h^2} + y(t) \end{displaymath}

は大体 $0$ に等しい. したがって,

\begin{displaymath}y(t+h) = 2y(t)-y(t-h) - h^2 y(t) \end{displaymath}

が近似的になりたつとしてよいであろう. $Y_k = y(kh)$, $k=0, 1, 2, \ldots$ で数列 $Y_k$ を定める. このとき,

\begin{displaymath}Y_{k+2} = 2 Y_{k+1} - Y_k - h^2 Y_{k} \end{displaymath}

なる漸化式が近似的になりたつとしてよい. したがって, 漸化式
\begin{displaymath}
y_{k+2} = 2 y_{k+1} - y_k - h^2 h_{k+1}
\end{displaymath} (9.2)

を満たす数列を決めることにより, 解の近似をもとめることが可能であると予想できる. このように解の近似を求める方法を差分法とよぶ. 漸化式 (9.2) を元の微分方程式の 差分化, または差分スキームとよぶ.

数列の値を決めるには初期条件が必要であるが, それは, 微分方程式の初期条件 $y(0) = a, y'(0) = b$ を用いて 関係式

\begin{displaymath}y(0) = a = y_0 , \
y'(0)= b \simeq \frac{y_1-y_0}{h}
\end{displaymath}

できめてやればよい.

例題 9.1   強制振動の方程式

\begin{displaymath}{{d^2} \over {dt^2}} y + y = \sin(at) \end{displaymath}

を数値的に解いてみよう. $\sin(at)$ が外部から単振動する物体に与えられる力である.

load("glib")$

def osci() {
  glib_window(0,-5,50,5);
  glib_clear();
  glib_line(0,0,50,0);
  glib_line(0,-10,0,10);

  X1 = 0.5; X2 = 0.501; A=0.5;

  Dt = 0.07; T = 0;
  while (T<50) {
     X3 = 2*X2-X1+Dt*Dt*(eval(sin(A*T))-X2);
     glib_putpixel(T,X1);
    /*     print([T,X1]); */
     T=T+Dt;
     X1=X2; X2=X3;
  }
}

print("Type in osci()")$
end$

これで, 微分方程式を近似的に解く問題が, 漸化式をみたす数列を求める 問題になったのであるが, このような近似解が本当の解に収束するかとか, 全くことなる解しかとらえられない不安定現象が起こる場合が あるとかの議論をやらないといけない.


next up previous contents index
: 不安定現象 : 常微分方程式の数値解 : 常微分方程式の数値解   目次   索引
Nobuki Takayama 平成15年9月12日