（以下では例によって、開発環境が手元にない場合は、paiza.io, onlineGDB, programiz, codechef.com, ideone.com, paiza.ioなどを実行してみると良い）
連立常微分方程式も、同様の手法で解くことができる。たとえば生物の捕食―被食関係における個体数の変動を表現する数理モデルに、ロトカ・ヴォルテラの方程式がある。捕食者（たとえばヤマネコ）の個体数$x$が、被食者（たとえばウサギ）の個体数$y$とが、以下のような2元連立非線形常微分方程式の関係であらわされるとする考え方である：
$$\frac{dx}{dt}=ax – cxy$$
$$\frac{dy}{dt}=-by+cxy$$
こうした連立微分方程式であっても、同様にルンゲ・クッタ法を適用することで、現実的に解くことができる。やり方としては同じだが
$$\frac{dx}{dt}=f(t,x,y)$$
$$\frac{dy}{dt}=g(t,x,y)$$
と整理して、各々に1段ずつルンゲクッタ法を適用していけばよい。すなわち、
$$\begin{matrix}k_{x1}=hf(t_i,x_i,y_i) & k_{y1}=hg(t_i,x_i,y_i)\\k_{x2}=hf(t_i+\frac{h}{2},x_i+\frac{k_{x1}}{2},y_i+\frac{k_{y1}}{2})&k_{y2}=hg(t_i+\frac{h}{2},x_i+\frac{k_{x1}}{2},y_i+\frac{k_{y1}}{2}) \\k_{x3}=hf(t_i+\frac{h}{2},x_i+\frac{k_{x2}}{2},y_i+\frac{k_{y2}}{2})&k_{y3}=hg(t_i+\frac{h}{2},x_i+\frac{k_{x2}}{2},y_i+\frac{k_{y2}}{2})\\k_{x4}=hf(t_i+h,x_i+k_{x3},y_i+k_{y3}) &k_{y4}=hg(t_i+h,x_i+k_{x3},y_i+k_{y3})\end{matrix}$$
として順次計算し、
$$x_{i+1}=x_i+\frac{1}{6}(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})$$
$$y_{i+1}=y_i+\frac{1}{6}(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})$$
として次のステップの $x, y$を計算していけばよい。一見複雑に見えるが、やっていることは極めて単純で、これをプログラムにすると
#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=（ここに記述する）; ky1=（ここに記述する）;
kx2=（ここに記述する）; ky2=（ここに記述する）;
kx3=（ここに記述する）; ky3=（ここに記述する）;
kx4=（ここに記述する）; ky4=（ここに記述する）;
x=（ここに記述する）
y=（ここに記述する）
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回くらい繰り返すと、
のような餌と捕食者の変動の周期的な関係を見ることができる（モデルやパラメタの問題点は議論しない）。
これで常微分方程式に基づいた、さまざまな物理シミュレーションを実施するための基礎的な知識は身に付いた。なお、上記のコードは、ポインタの概念を使えば関数の読み出しの部分で断然簡素化することができるので、興味のある人はやってみよう。
今週の課題
課題は講義中にやってもらうことを想定しています、以下について回答してみてください。オリジナリティが高い回答ほど高評価で、おなじオリジナリティであれば提出時間が早いほど高評価。提出は10/10 深夜まで。提出が遅くてもオリジナリティが非常に高いものが最高評価、とします。提出はITCではなくて、https://forms.gle/5SAxcgPD1kkkKZ7w9 にお願いします。
ルンゲクッタの４次公式は、以下のようなものであった：
微分方程式$\frac{dy}{dx}=f(x,y)$ に対し、初期条件$y(x_0)=y_0$としたときの数値解は、刻み幅を$h$として、$x_n=x_0+nh$における$y$の値を$y_n$から$y_{n+1}=y_n+k$と書く場合、$k$は以下で近似できる
$$\begin{array}{ccc}k_1 & = & hf(x_n, y_n)\\k_2 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_1}{2}\right )\\k_3 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_2}{2}\right )\\k_4&=&hf(x_n+h,y_n+k_3)\\k&=&\frac{1}{6}(k_1+2k_2+2k_3+k_4) \end{array}$$
連立微分方程式
$$\left \{ \begin{array}{ccc}\frac{dy}{dx} &=&f(x,y,z),&初期条件が& y(x_0)=y_0\\\frac{dz}{dx} &=&g(x,y,z),&初期条件が& z(x_0)=z_0 \end{array} \right .$$
の数値解は、刻み幅を$h$、$x_n=x_0+nh$として以下で与えられる：
$$\begin{array}{ccc}k_1 & = & hf(x_n, y_n, z_n)\\l_1 & = & hg(x_n, y_n, z_n)\\k_2 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_1}{2},z_n+\frac{l_1}{2}\right )\\l_2 & = & hg \left ( x_n+\frac{h}{2},y_n+\frac{k_1}{2},z_n+\frac{l_1}{2}\right )\\k_3 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_2}{2},z_n+\frac{l_2}{2}\right )\\l_3 & = & hg \left ( x_n+\frac{h}{2},y_n+\frac{k_2}{2},z_n+\frac{l_2}{2}\right )\\k_4&=&hf(x_n+h,y_n+k_3, z_n+l_3)\\l_4&=&hg(x_n+h,y_n+k_3, z_n+l_3)\\k&=&\frac{1}{6}(k_1+2k_2+2k_3+k_4)\\l&=&\frac{1}{6}(l_1+2l_2+2l_3+l_4) \\y_{n+1}&=&y_n+k\\z_{n+1}&=&z_n+l\end{array} $$
上の考え方を使って、以下の問題を解いてみよう
課題1：$y^\prime =y-10x +1$を刻み幅0.1で$x=0$から$x=1.0$まで解け。ただし初期条件を$y(0)=1$とする。
課題２：2階微分方程式$\frac{d^2 y}{dx^2}-2x\frac{dy}{dx}+xy=0$を初期条件$y(0)=1, y^\prime(0)=1$で解け。
課題3: 野球ボールがどこまで飛ぶか考えよ。ただし空気抵抗が速度の2乗に比例する以下の形で書けるとする。
$$\begin{array}{ccc}a_x&=&-\frac{3}{4D}\frac{\rho_{air}}{\rho_{ball}} C_d V V_x\\a_y&=&-g-\frac{3}{4D}\frac{\rho_{air}}{\rho_{ball}}C_d V V_y \\V&=&\sqrt{V_x^2+V_y^2}\end{array}$$
ここで$a_x$を水平方向の加速度、$a_y$を垂直方向の加速度、$D$をボールの直径=7[cm]、$g$を重力加速度=9.8[m/s$^2$]、$\rho_{air}$を空気の密度=1.2[kg/m$^3$]、$\rho_{ball}$をボールの密度=650[kg/m$^3$]、$C_d$を抵抗係数=0.2とする。初速を50[m/s]、バットで打った場所が高さ1.0[m]である場合、45度で飛び出す時が最大距離ではないことを何らかの形で示せ。なおボールの回転による揚力や、風など空気の揺らぎは考えない。
なお、以下のコードを参考にしても良い。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 1000
#define h 0.01
double G=-9.8;
double k=-0.00396;
double f(double t,double x, double y);
double g(double t,double x, double y);
int main(void) {
int i;
double v=50.;
double degree=38.;
double radian=degree*3.141592/180.0;
double x=v*cos(radian);
double y=v*sin(radian);
double t=0.0;
double xx=0.,yy=1.0;
double vx, vy,x1,x2,x3,x4,y1,y2,y3,y4;
double kx1,ky1,kx2,ky2,kx3,ky3,kx4,ky4;
for (i = 0; i < N; i++) {
x1=ここに記述;y1=ここに記述;
kx1=ここに記述;ここに記述;
vx=ここに記述;ここに記述;
x2=ここに記述;y2=ここに記述;
kx2=ここに記述;ky2=ここに記述;
vx=ここに記述;vy=ここに記述;
x3=ここに記述;y3=ここに記述;
kx3=ここに記述;ky3=ここに記述;
vx=ここに記述;vy=ここに記述;
x4=ここに記述;y4=ここに記述;
kx4=ここに記述;ky4=ここに記述;
xx=ここに記述;
yy=ここに記述;
x=ここに記述;
y=ここに記述;
t=t+h;
if (yy<0.0){
printf("max=%4.1lfm\t%4.4lf\t %4.4lf\t %e\t %4.4e\t %4.4e\n",xx,x,y,t,xx,yy);
exit(0);
}
}
}
double f(double t,double x, double y) {
return (k*sqrt(x*x+y*y)*x);
}
double g(double t,double x, double y) {
return (G+k*sqrt(x*x+y*y)*y);
}