講義第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項から後をすべて無視して、単に

としたものである。

この、線形の方程式の場合にはもちろん解析的に解けるので、数値解を計算したりする必要はない。しかし、非線形で解析解がない問題というのは無数にあり、そのときにこのような数値的な近似は強力な手段となる。

練習

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

ここで、 y0 はxでの値、y_1 は x+Δxでの値で、ypは途中で使うだけである。この方法のプログラムを作って、同じ問題を解いてみよう。誤差はどうなるだろうか。
厳密解、近似解をグラフィック表示してみよう。座標軸、目盛なども工夫して書いてみること。

連立常微分方程式

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 の２つがある。 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.

練習

Van der Polの方程式

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

をプログラムしてみよう。例えば y[1]} を x, y[2]} を v_x, y[3] を y, y[4]} を v_y とおき、 calculate_dy を書き直す。また、 nvar=4 とする。
Lotka-Volterra の方程式

をプログラムして、数値解を求めてみよう。また、定常状態はどこか？
2次元ケプラー問題の場合に、ルンゲ・クッタと台形公式のそれぞれで、 Δxを変えたときに誤差がどう変わるかを調べてみよう。また、 1周期積分した時と、100 周期、 10000 周期と伸ばしていった時について、公式による違いを調べて見よう。
最初の例の一階線形常微分方程式で、系数 k の値を大きな負の数にしたとき、オイラー法、ルンゲ・クッタ、台形公式のそれぞれでどのような解が求まるかを調べて見よう。うまく解が求まらないのはどういう場合か、またそれは何故か。