課題１

速度に関する微分方程式として考えてみて下さい．x=vと書くべきだったかもしれませんが，例に合わ
せたのと，微分方程式のアルゴリズム上は単なる1階の微分方程式と考えていただければ充分なの
です．

関数は

double fnc(double t, double x){
　double a=0.1;
　double b=0.02;
　double g=9.807;
　return -1.0 * b * x * x - a * x + g;
}

tは出てきません．x'(t)=sin(5.0 * t)のような場合はtを使いますので，このような関数に対応するための
形式を残しておいたのがxとtを引数にした理由です．

ルンゲクッタ法

// 課題1
#include <stdio.h>

const double a = 0.1;
const double b = 0.02;
const double g = 9.807;

// 1階微分方程式
double fnc(double t,double x)
{
　return -b*x*x - a*x +g;
}

// ルンゲクッタ法
double rkm(double(* fnc)(double,double), double t, double x, double h){
double ck0, ck1, ck2, ck3;

　ck0=(*fnc)(t, x)*h;
　ck1=(*fnc)(t+h*0.5, x+ck0*0.5)*h;
　ck2=(*fnc)(t+h*0.5, x+ck1*0.5)*h;
　ck3=(*fnc)(t+h, x+ck2)*h;

　return (x+(ck0+(ck1+ck2)*2.0+ck3)/6.0);
}

int main()
{
　int i;
　double t0 = 0.0; // 始点
　double tn = 5.0; // 終点
　double h = 0.01; // 刻み
　int n = (int)((tn - t0) / h); // 刻み数
　double t;
　double x=0.0;

　for(i=0;i<n;i++){
　　t = t0 + i*h;
　　x = rkm(fnc,t,x,h);
　　if(i%10==0) printf("%lf\n",x);
　}
　return;
}

ルンゲクッタ・ジル法

// 課題1
#include <stdio.h>

const double a = 0.1;
const double b = 0.02;
const double g = 9.807;

// 1階微分方程式
double fnc(double t,double x)
{
　return -b*x*x - a*x +g;
}

// ルンゲクッタ・ジル法
double rkgm(double (*fnc)(double, double), double t, double x, double h, double *q){
　static double p1=0.2928932188134524, p2=1.70710678118547;
　double ck0, ck1, ck2, ck3;
　double q1, q2, q3;
　double r1, r2, r3, r4;
　double tt, x1, x2, x3;

　ck0 = (*fnc)(t,x) * h;
　r1 = ck0 * 0.5 - *q;
　x1 = x + r1;
　q1 = *q + r1 * 3.0 - ck0 * 0.5;
　tt = t + h * 0.5;
　ck1 = (*fnc)(tt,x1) * h;
　r2 = (ck1-q1) * p1;
　x2 = x1 + r2;
　q2 = q1 + r2 *3.0 - ck1 * p1;
　ck2 = (*fnc)(tt,x2) * h;
　r3 = (ck2-q2) * p2;
　x3 = x2 + r3; // x2=x2+r3になっていたのを修正
　q3 = q2 + r3 * 3.0 - ck2 * p2; // q2=q2+r3*3.0-ck2*p2になっていたのを修正
　tt = t + h;
　ck3 = (*fnc)(tt,x3) * h;
　r4 = (ck3-q3-q3) / 6.0;
　*q = q3 + r4 * 3.0 - ck3 * 0.5;

　return (x3+r4);

}

int main()
{
　int i;
　double t0 = 0.0; // 始点
　double tn = 5.0; // 終点
　double h = 0.01; // 刻み
　int n = (int)((tn - t0) / h); // 刻み数
　double t;
　double x=0.0;
　double q = 0;

for(i=0;i<n;i++){
　t = t0 + i*h;
　x = rkgm(fnc,t,x,h,&q);
　if(i%10==0) printf("%lf\n",x);
　}
　return;
}

なお，授業中で省略した，ミルンの方法とアダムスの方法は，授業の始めで説明します．

課題２

細部は異なりますが基本の流れは同様です．

// 課題2
#include <stdio.h>

const double cm = 0.02;
const double km = 0.1;
const double um = 1;

// 1階微分方程式
double f(double t,double x,double z)
{
　return z;
}

// 2階微分方程式
double g(double t,double x,double z)
{
　return -cm*x - km*x + um;
}

// ルンゲクッタ法
double rkm_2(double(* f)(double,double,double),double(* g)(double,double,double), double t, double x, double h, double z){
double k, k0, k1, k2, k3, l, l0, l1, l2, l3;

　k0 = h*(*f)(t, x, z);
　l0 = h*(*g)(t, x, z);
　k1 = h*(*f)(t + h/2, x + k0/2, z + l0/2);
　l1 = h*(*g)(t + h/2, x + k0/2, z + l0/2);
　k2 = h*(*f)(t + h/2, x + k1/2, z + l1/2);
　l2 = h*(*g)(t + h/2, x + k1/2, z + l1/2);
　k3 = h*(*f)(t + h, x + k2, z + l2);
　l3 = h*(*g)(t + h, x + k2, z + l2);

　k = (k0 + 2*k1 + 2*k2 + k3)/6.0;
　l = (l0 + 2*l1 + 2*l2 + l3)/6.0;

　x = x + k;
　z = z + l;

　return x;
}

int main()
{
　int i;
　double t0 = 0.0; // 始点
　double tn = 5.0; // 終点
　double h = 0.01; // 刻み
　int n = (int)((tn - t0) / h); // 刻み数
　double t;
　double z = 1.0;
　double x = 0.05;

　for(i=0;i<n;i++){
　t = t0 + i*h;
　x = rkm_2(f,g,t,x,h,z);
　if(i%10==0) printf("%lf\n",x);
　}
　return;
}