11回目

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

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

n=5 比率 n=50 比率 n=500 比率
*****fnc1****
台形則 3.84 1 3.8334 1 3.833334 1
ガウス則 3.833333 0.998263802 3.833333 0.999982522 3.833333 0.999999739
ニュートン則 1.95 0.5078125 3.5988 0.938801064 3.809388 0.993753218
シンプソン則 2.794667 0.727777865 3.755979 0.979803569 3.82536 0.997919826
*****fnc2****
台形則 9.561919 1 9.56761 1 9.567667 1
ガウス則 9.567668 1.000601239 9.567668 1.000006062 9.567668 1.000000105
ニュートン則 5.662929 0.592237709 9.173388 0.958796188 9.528212 0.995876215
シンプソン則 7.60533 0.795376953 9.436289 0.98627442 9.554516 0.998625475
*****fnc3****
台形則 1.128858 1 1.39955 1 1.401666 1
ガウス則 1.622695 1.437466005 1.622695 1.159440534 1.622695 1.157690206
ニュートン則 4.199242 3.719902769 1.524332 1.089158658 1.41481 1.009377412
シンプソン則 2.970974 2.63184032 1.443382 1.031318638 1.406026 1.003110584

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

n=50000 比率
*****fnc1****
台形則 3.833333 1 もともと誤差は少なかったが
殆ど誤差はなくなっている．
アルゴリズムの差があまり出ないことが分かる．
ガウス則 3.833333 1
ニュートン則 3.833309 0.999993739
シンプソン則 3.833325 0.999997913
*****fnc2****
台形則 9.567668 1 回数を繰り返すことで
殆ど誤差はなくなっている．
回数を繰り返すことで解消できることが分かる．
ガウス則 9.567668 1
ニュートン則 9.567628 0.999995819
シンプソン則 9.567654 0.999998537
*****fnc3****
台形則 1.401687 1 回数を増やしても，ガウス則の誤差はある．
fnc3に対してガウス則は不適合であると分かる．
ガウス則 1.622695 1.157672861
ニュートン則 1.401701 1.000009988
シンプソン則 1.401692 1.000003567

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

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

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〜何秒行うのか，刻み幅はどのくらいにするのかを決めて計算します．

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

#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;

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

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


課題２　以下の関数を作成し，上記のオイラー法によって常微分方程式を解け．
また，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)として扱います．

提出するデータは
課題１のオイラー法のアルゴリズムの関数を収めたソースファイル
課題２のシミュレーション結果と考察を記したデータ
である．

※fncのtは要らないか？
この式であれば要りません．double fnc(double ｘ)でも構いません．右辺に例えば，tの関数となるもの
例えばsin(5.0*t)があるときには，double fnc(double t, double x)とする必要があります．