11回目　(2015年度版)

今回から数値計算を扱っていきます．
前回は，入り口として数値微分と数値積分をやってみました．
手で計算しても分かるものでしたが，細かい誤差が生じたことを確認していただけたと思います．

前回の計算例を示します．
意図的に誤差が大きくなるようなシミュレーションを行っていますが，以下のような比となります．
5%以内の誤差であれば，一応参考になると考えても良いのですが，10%ともなると実用的とは
言えず，誤差なのか，条件が不適切なのか，もともとのモデルが不充分なのかなどが分からなく
なってしまいます．

n=5では大きくばらついていますし，n=50でもまだ違いがあります．n=500にしてもfnc3などは
まだ誤差があります．

従って，数値演算を行う際には，以下の点について留意しておく必要があります．

• 数値演算のアルゴリズムにより扱う式や条件に向き不向きがあることを認識する．
• 算出された数値が正しいかどうかについて検証する．
• 近似式や簡素化したモデルの解析解と比較してみる．
• いくつかのアルゴリズムによる解を比較してみる．

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

または，

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

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

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

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

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

常微分方程式の数値解法
今回は，一般的に使われながら，より精度の高いルンゲクッタ法と呼ばれるもので，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)

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

補足説明
ミルンの方法とアダムスの方法について，補足説明をします．
ルンゲクッタ法は基本になりますが，4次までの近似で，なおかつ多角形的なものであることはオイラー法と
変わりがありません．
その欠点を経験則的に解決したのがルンゲクッタ・ジル法であり，係数が細かく調整してあります．それでも，
微分方程式の曲線のカーブの変化が大きい時などは誤差が蓄積します．
ミルンの方法，アダムスの方法はこの問題点を解決する方法であり，関数中に設けたepsというパラメータが
あることに着目してみて下さい．
これは設定した誤差であり，一ステップ当たりの誤差がeps以下になるまでルンゲクッタによる解を修正すると
いうものです．
ルンゲクッタによる４ステップ分の計算をまず行い，これらの解を基に修正を１ステップごとに繰り返していくと
いうものです．ですから計算量は多くなります．epsを小さくすればするほど計算量が増える可能性が大きく
なるといえます．
ミルンの方法→comp_mlnm.c
アダムスの方法→comp_admsm.c

課題１　以下の関数を作成しし，上記の4つの方法によって常微分方程式を解け．
関数に関するヒントはここにある．
4つのファイルをダウンロードし，ソースファイルを完成せよ．また解いた結果を比較し，考察せよ．

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)

課題２　以下の関数について，ルンゲクッタ法を用いて常微分方程式を解け．

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]として計算してみよ．
数値は変えてもかまわない．

ルンゲクッタ法→rkm.c　をダウンロードし，名前を変えて先のアルゴリズムに合わせて変更すればよい．
要は，上記のアルゴリズムをそのまま並べて書けば充分である．rkm_2.cのような名前にしておけばよい．
運動方程式の関数は dx/dt = f(t, x, z)，dz/dt = g(t, x, z)となることに注意せよ．

課題に取り組むにあたり，ソースファイルは以下のものを参考にしてもよい．
ミルンの方法　→comp_mlnm.c
アダムスの方法　→comp_admsm.c