講義第12回:Pascal プログラミング --- 微分方程式

今回は、計算機の使い道のわりと重要なものの一つである、シミュレーション、 具体的には数値的に常微分方程式を解く方法について学ぶ。説明は、まず一変 数の常微分方程式について数値的に解く原理を説明し、多変数の常微分方程式 に進む。

時間の関係があるので偏微分方程式は扱わない。

常微分方程式 --- 1階線形

program ode1(input, output);
var x, xend, dx,  y, k : real;
begin
    write('Enter k, y0, xend, dx:');
    readln(k, y, xend, dx);
    x := 0;
    while x <= xend do begin
        y := y + dx*k*y;
        x := x + dx;
        writeln('x, y = ', x:20:10,' ', y:20:10);
   end;
end.
上の例は、一階線形常微分方程式の初期値問題

の近似解を求めるプログラムである。計算法は、yのテイラー展開

で、 第3項から後をすべて無視して、単に

としたものである。

練習

  1. 例えば k=1, y0=1 とした時、微分方程式の厳密解は y = exp(x) になる。 x=1 での数値解と厳密解の差を調べて見よう。それは Δx のどのような関数になっているだろうか。また、それはなぜだろ うか。
  2. 上の例で説明した積分法(前進オイラー法)は、プログラムが簡単 であるかわりに誤差が大きい(Δx を小さくしてもなかなか精度が上 がらない)。もう少しましな方法に、2次のルンゲ・クッタ(それぞれ人の名 前)法と呼ばれるものがある。これは、以下のような方法である。

    ここで、 y0 はxでの値、y_1 は x+Δxでの値で、ypは途中 で使うだけである。この方法のプログラムを作って、同じ問題を解いてみよう。 誤差はどうなるだろうか。

  3. 厳密解、近似解をグラフィック表示してみよう。座標軸、目盛など も工夫して書いてみること。

連立常微分方程式

2次ルンゲ・クッタ法

program ode_system(input, output);
const nvar = 2;
type vector = array[1..nvar] of real; 
var i		   : integer;
   y		   : vector;
   x		   : real;
   xinit, xend, dx : real;
   
procedure calculate_dy(x: real; y : vector; var dy : vector);
begin
   dy[1] := y[2];
   dy[2] := -y[1];
end; 

procedure integrate(var y  : vector; var x  : real; dx : real);
var yp,dy : vector;
   i	  : integer;
begin
   calculate_dy(x,y,dy);
   for i:=1 to nvar do yp[i] := y[i]+0.5*dx*dy[i];
   calculate_dy(x,yp,dy);
   for i:=1 to nvar do y[i] := y[i]+dx*dy[i];
   x := x + dx;
end; 

begin
   write('Enter  xend, dx: ');
   readln(xinit, xend, dx);
   y[1]:=0;
   y[2]:=1;
   x := 0;
   while x < xend do begin
      integrate(y, x, dx);
      write(x:16:7);
      for i:= 1 to nvar do write(' ',y[i]:15);
      write(y[1]-sin(x), y[2]-cos(x));
      writeln;
   end;
end.
このプログラムは、(手続き calculate_dy の中身さえ書き換えれば) どのような連立常微分方程式でも解くことができる(はずである)。 今までに使ったことのない Pascal の言語機能としては
const nvar = 2;
type vector = array[1..nvar];
const type の2つがある。 const は、プログラム のあちこちで同じ値になる定数に名前をつける機能である。プログラムが読み やすくなるだけでなく、この例の場合では方程式の本数を変えるのに、ここの 宣言を変えるだけで済むことになる。

type というのは、上のように「大きさ 100 の配列」というものを、実数 とか整数と同じような「型」と宣言する機能である。これにより、いちいち array ... of ... と書かなくて済むし、また配列を手続きの引数にす ることができる。

プログラムでやっていることは、手続き、配列を使っていることの他は上でやっ た1変数常微分方程式とあまり違わない。

台形公式

procedure integrate(var y  : vector;
		    var x  : real;
			dx : real);
var y1,dy0,dy1 : vector;
   i,k	  : integer;
begin
   calculate_dy(x,y,dy0);
   for i:=1 to nvar do y1[i] := y[i]+dx*dy0[i];
   for k:=1 to 3 do begin 
      calculate_dy(x,y1,dy1);
      for i:=1 to nvar do y1[i] := y[i]+0.5*dx*(dy0[i] + dy1[i]);
   end;
   for i:=1 to nvar do y[i] := y1[i];
   x := x + dx;
end; 
台形公式は、ルンゲ・クッタ公式と同様に高い精度が得られる。この方式の特徴は
y1 = y + 0.5*dx*(dy0 + dy1)
としているところで、 dy1 が実は y1 の関数になっているところである。こ のために、ここでは前回やった「直接代入法」で、 y1 についての方程式を解 いている。ただし、 答が正しく求まったかどうかの判定はしないで、最初の 値(オイラー法でもとめた)から数回(この例では3回)繰り返して打ちきっ ている。

台形公式は、計算量が増えるという欠点はあるが、重要な特長がいくつかあり、 この方法を発展させた様々な方法が広く用いられている。

グラフィック

begin
   write('Enter xend, dx: ');
   readln(xend, dx);
   y[1]:=0;
   y[2]:=1;
   x := 0;
   initgraph(gdriver, gmode, '');
   moveto(400+round(200*y[1]),400-round(200*y[2]));
   while x < xend do begin
      integrate(y, x, dx);
      lineto(400+round(200*y[1]),400-round(200*y[2]));
   end;
   readln;
end.
メインプログラムを上のように書き換えると、数値積分した答をグラフィック 表示することができる。(gdriver, gmode の宣言を忘れないように)

練習

独立変数を横軸、y[1] を縦軸にするグラフを書いてみよう。また、縦 軸、横軸を自由に指定できるようにするにはどうすればいいか考えて見よう。

一変数高階微分方程式

一般に、高階微分方程式 f(x,y,y',y'',y''' ,...) = 0数値解は、最高次の次数を nとして以下のように計算できる:

例:強制振動

例えば強制振動の方程式

は、以下のようにcalculate_dy を書き換えれば計算できる。
(省略)
 a, b, c,omega	   : real;
procedure calculate_dy(x : real;y : vector; var dy : vector);
begin
   dy[1] := y[2];
   dy[2] := -2*b*y[2] -c*y[1] + a* sin(omega*x);
end;
(省略)
begin
   write('Enter xinit, xend, dx: ');
   readln(xinit, xend, dx);
   write('Enter a, b, c,  omega ');
   readln(a, b,c, omega);
   y[1]:=1;
   y[2]:=0;
   x := xinit;
   initgraph(gdriver, gmode, '');
   moveto(400+round(200*y[1]),400-round(200*y[2]));
   while x < xend do begin
      integrate(y, x, dx);
      lineto(400+round(200*y[1]),400-round(200*y[2]));
   end;
   readln;
end.

練習

  1. Van der Polの方程式

    をプログラムしてみよう。これは非線形振動(自励振動)のもっとも簡単なモ デルの一つである。
  2. 2次元ケプラー問題

    をプログラムしてみよう。例えば y[1]} を x, y[2]} を v_x, y[3] を y, y[4]} を v_y とおき、 calculate_dy を書き 直す。また、 nvar=4 とする。
  3. Lotka-Volterra の方程式
    をプログラムして、数値解を求めてみよう。また、定常状態はどこか?
  4. 2次元ケプラー問題の場合に、ルンゲ・クッタと台形公式のそれぞれで、 Δxを変えたときに誤差がどう変わるかを調べてみよう。また、 1周期積 分した時と、100 周期、 10000 周期と伸ばしていった時について、公式 による違いを調べて見よう。

次回予告

次回は、計算機を使ってはじめて実用になるような、「モンテカルロ法」と呼 ばれる問題の解き方について学ぶ。