12回目　微分方程式の数値演算(2)

関数ポインタ
変数を関数で渡す方法は既に学んでいると思いますが，実は関数もポインタを利用してやりとりをする
ことが可能なのです．
以下に例を示したいと思います．

#include <stdio.h>

double fnc(double x){
　return x;
}


main(){
　double x = 0.0;
　double r = 0.0;
　double (*pfnc)(double);

　pfnc = fnc;

　x = 10.0;

　r = (*pfnc)(x);

　printf("%f\n",r);
}

通常の関数を宣言します．







メイン関数内で関数ポインタを宣言します．

先の関数のアドレスを関数ポインタ(アドレス)に送ります．



関数ポインタの中身に引数を渡し，先の関数を実行します．

または，

#include <stdio.h>

double fnc(double x){
　return x;
}


main(){
　double x = 0.0;
　double r = 0.0;
　double (*pf)(double);

　pf = *fnc;

　x = 10.0;

　r = pf(x);

　printf("%f\n",r);
}

通常の関数を宣言します．







メイン関数内で関数ポインタを宣言します．

先の関数を関数のポインタの中身に送ります．



関数ポインタに引数を渡し，先の関数を実行します．

のようにしてもやりとりができます．
関数をわざわざポインタにする理由は何かというと，

  1. 関数の再利用性を上げることができる
  2. 関数を動的に何度でも呼び出して利用できる
  3. 関数に対して異なった種類の関数を呼び込むことができる．

既に作っておいた関数を利用する際，ポインタで受け渡しをすることで，例えば関数名が合わなくても
その制約をなくして利用できるので再利用性が向上します．また，関数をモジュール（ユニット）のような
ものとして利用できるので関数内のアルゴリズムの構造を何度も形を変えて利用できるという点があり
ます．何かの数値を算出するための関数ではなく，計算するための仕組みとしての関数にでき，汎用性，
再利用性が格段に上がります．

授業の課題であれば，プログラムが小規模な上に，あまり再利用する機会がないのでメリットがつかみ
にくいと思いますが，自作のプログラムを作る際は大変便利ですので参考にしてみて下さい．

関数ポインタを関数に呼び込む方法が分かりにくいと思いますので，以下の例を見て下さい．
fnc_1の引数として関数ポインタを用いています．
メイン関数内では関数fnc_A, fnc_Bをこのfnc_1に渡しています．
fnc_1の引数(*fnc)(double)は関数ポインタの中身を表しているので，
(*fnc) = fnc_Aとして関数ごと引き渡すことができると考えて下さい．

#include <stdio.h>

double fnc_1(double (*fnc)(double), double a){
　return (*fnc)(a);
}

double fnc_A(double a){
　return a;
}

double fnc_B(double a){
　return a*a;
}

main(){
　double h = 2.0;
　printf("fnc_A = %f\n", fnc_1(fnc_A, h));
　printf("fnc_B = %f\n", fnc_1(fnc_B, h));
}

関数ごとポインタとして引数にする関数



通常の関数fnc_A



通常の関数fnc_B





fnc_1にfnc_Aを読み込んで変数を与える
fnc_1にfnc_Bを読み込んで変数を与える

前回の常微分方程式の解法はシンプルであり，原理も分かりやすいものであると思いますが，計算
精度の面で限界があります．今回は，一般的に使われながら，より精度の高い方法を説明します．
ルンゲクッタ法と呼ばれるもので，1階の常微分方程式dx(t)/dt = f(t, x)を解く際に利用できます．
時間刻みΔtを便宜上hとおくとして，h後の時刻と解は次のようになります．

ti+1 = ti + h
xi+1 = xi + k
k=(1/6)(k0+2k1+2k2+k3)

最後の式がルンゲクッタの解法の特徴ですが，kはそれぞれ，

k0 = h f(ti, xi)
k1 = h f(ti + h/2, xi + k0/2)
k2 = h f(ti + h/2, xi + k1/2)
k3 = h f(ti + h, xi + k2)

となります．複雑な計算をしているように見えますが，離散的ながら，重みを付けて4次の曲線として
数値解を出していると考えていただければよいかと思います．
なお，アルゴリズムは，

  1. 初期値t0, x0を決定し，時間刻みhを適宜決めます．
  2. k0 = h f(ti, xi)を計算します．
  3. k1 = h f(ti + h/2, xi + k0/2)を計算します．
  4. k2 = h f(ti + h/2, xi + k1/2)を計算します．
  5. k3 = h f(ti + h, xi + k2)を計算します．
  6. k=(1/6)(k0+2k1+2k2+k3)を計算し，次のステップ ti+1 = ti + hのkとします．
  7. xi+1 = xi + kでh刻み後の解が算出できます．

あとは2〜7を繰り返します．for文などで繰り返せばよいかと思います．
以下にソースファイルがありますので，適宜利用して下さい．

ルンゲクッタ法→rkm.c
ルンゲクッタ・ジル法→rkgm.c
ミルンの方法→mlnm.c
アダムス法→admsm.c

いずれもルンゲクッタのアルゴリズムを基本としていますが，係数を変えたり，計算を追加して，精度を
上げる工夫をしたものです．但し，計算量は少しとはいえ増える傾向にあります．

2階微分方程式

ルンゲクッタ法は1階の微分方程式を解くための方法です．工夫すれば，連立した1階の微分方程式を
解くことができます．
機械工学で利用する微分方程式は2階もしくはそれ以上の高次のものが大半です．これを解くにはどの
ようにしたらよいかというと，この連立して解く方法を利用するのです．
dx/dt=zとおいて，導関数を新たに変数にしてしまうことにします．x'= zだからx"= z'となるわけです．
従って，新たに変数が発生したため，

dx/dt = f(t, x, z)
dz/dt = g(t, x, z)

の連立方程式を立てて同時に解けばいいのです．

ti+1 = ti + h
xi+1 = xi + k
zi+1 = xi + l
k=(1/6)(k0+2k1+2k2+k3)
l=(1/6)(l0+2l1+2l2+l3)

但し

k0 = h f(ti, xi, zi)
l0 = h g(ti, xi, zi)
k1 = h f(ti + h/2, xi + k0/2, zi + l0/2)
l1 = h g(ti + h/2, xi + k0/2, zi + l0/2)
k2 = h f(ti + h/2, xi + k1/2, zi + l1/2)
l2 = h g(ti + h/2, xi + k1/2, zi + l1/2)
k3 = h f(ti + h, xi + k2, zi + l2)
l3 = h g(ti + h, xi + k2, zi + l2)

のような計算をすればいいのです．式は増えていますがやっていることは変わりません．

微分方程式の数値解法(オイラー法)
以下のような微分方程式（導関数）があるとします．時間を軸に考えてみましょう．

dx/dt=f'(t)

微分方程式を差分にすると以下のように近似できます．微小区間の傾斜が導関数ですので，この
微小区間を小さくすればよいと考えられます．

dx/dt≒Δy/Δt=(f(t+Δt)-f(t))/Δt

Δtは極めて小さいものとすれば即ち微分した値となるので，元の微分方程式は以下のように近似
できます．

(f(t+Δt)-f(t))/Δt≒f'(t)

この式を変形すれば，

(f(t+Δt)≒f(t)+Δt f'(t)

となります．Δt後の値は元の値に幅Δtに変化f'(t)をかけたものを足したものと言えば，離散量らしく
計算できることは想像できると思われます．これがオイラー法の基礎的な考え方なのです．

微分方程式を数値的に解くということは，初期状態（初期値）をもとにある刻み幅で変化させてその
変化を積み重ねていくことです．
アルゴリズムは大変シンプルであり，以下のことをすればいいのです．

(1)f(Δt*0)=初期値：　初期条件を設定する．
(2)f(Δt*1)=f(Δt*0)+Δt*f'(Δt*0)：　(1)の結果を利用
(3)f(Δt*2)=f(Δt*1)+Δt*f'(Δt*1)：　(2)の結果を利用
これを繰り返していきます．

決めるものはf'(t)の式，刻み幅Δtです．
シミュレーションにあたっては，t=0〜何秒行うのか，刻み幅はどのくらいにするのかを決めて計算します．

課題１　上記のアルゴリズムを基にオイラー法の関数を作成せよ．
オイラー法の関数は以下のように書いてください．
double euler(double (*fnc)(double, double), double t, double x, double dt)
引数は
double(* fnc)(double, double)：　やりとりする関数のポインタを示します．
double t：　時刻を示します．
double x：　変数を示します．
double dt：　時間刻みを示します．
となっています．


課題２　以下の関数を作成し，上記のオイラー法によって常微分方程式を解け．
また，hの刻み幅を変化させてどのように解（曲線）が変化するか確かめよ．

x(t)'= -b*x(t)*x(t) - a*x(t)+g

これは速度に比例する抵抗と速度の2乗の抵抗を受ける物体が落下する系に関する運動方程式である．
定数はa=0.1 [1/sec], b=0.02 [1/m], g=9.807 [m/sec/sec]とするがa,bは任意．
t=0〜任意の時間について解けばよい．

xの時間関数であるから以下のような引数として，関数を作ればよい．

double fnc(double t, double x)

よくあるプログラムの説明にはdouble fnc(double x, double y)と書かれていますが，y=f(x)で示すより，
運動を示す時間の関数の方が分かりやすいと思いますので，x=f(t)として扱います．

※　ソースとグラフのファイル必ず添付のこと

課題3　以下の関数について，オイラー法を用いて常微分方程式を解け．

x"(t) = - (c/m) x'(t) - (k/m) x(t) + u/m

ばね-マス-ダンパー系でc/m = 0.02[N sec/m/kg], k/m = 0.1[N/m/kg], u/m = 1[N/kg]位をパラメータとし，
初期値はx'(t)= 1[m/sec], x(t) = 0.05[m]として計算してみよ．
数値は変えてもかまわない．

運動方程式の関数は dx/dt = f(t, x, z)，dz/dt = g(t, x, z)となることに注意せよ．

※　ソースとグラフのファイル必ず添付のこと