4. 4章の答え

(以下では、開発環境が手元にない場合は、paiza.io, onlineGDB, programiz, codechef.com, ideone.com, paiza.ioなどを実行してみると良い。)

例えば以下のようなものになる

#include <stdio.h>
#include <math.h>
#define N 1000
#define h 1.0
double f(double t,double x, double y);
double g(double t,double x, double y);
int main(void) {
  int i;
  double x=1000.0;
  double y=100.0;
  double t=0.0;
  double kx1,ky1,kx2,ky2,kx3,ky3,kx4,ky4;
  for (i = 0; i < N; i++) {
    kx1=h*f(t,x,y);ky1=h*g(t,x,y);
    kx2=h*f(t+0.5*h,x+kx1*0.5,y+ky1*0.5);ky2=h*g(t+0.5*h,x+kx1*0.5,y+ky1*0.5);
    kx3=h*f(t+0.5*h,x+kx2*0.5,y+ky2*0.5);ky3=h*g(t+0.5*h,x+kx2*0.5,y+ky2*0.5);
    kx4=h*f(t+h,x+kx3,y+ky3);ky4=h*g(t+h,x+kx3,y+ky3);
    x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0;
    y=y+(ky1+2.0*ky2+2.0*ky3+ky4)/6.0;
    t=t+h;
    printf("%4.4lf\t %4.4lf\t %e\n",x,y,t);
    }
}
double f(double t,double x, double y) {
  return (0.01*x-0.00001*x*y);
}
double g(double t,double x, double y) {
  return (-0.02*y+0.00001*x*y);
}

この結果が

1009.0456 99.0095 0
1018.1831 98.0376 1
1027.4132 97.0843 2
1036.7367 96.1491 3
1046.1546 95.2318 4
1055.6675 94.3323 5
1065.2765 93.4501 6
1074.9824 92.5852 7
1084.7859 91.7372 8
1094.6881 90.9059 9

のようになっていれば、計算はできている。1000回くらい繰り返すと、

のような餌と捕食者の変動の周期的な関係を見ることができる(モデルやパラメタの問題点は議論しない)。

これで常微分方程式に基づいた、さまざまな物理シミュレーションを実施するための基礎的な知識は身に付いた。なお、上記のコードは、ポインタの概念を使えば関数の読み出しの部分で断然簡素化することができるので、興味のある人はやってみよう。