多分これが原因

スプライン補間の近似すげーと思った。
元の離散的なデータから滑らかな曲線で近似できる。

でも、この「滑らか」という部分で勘違いをしてしまっていたようだ。
具体的には、移動平均みたいな、滑らかな「トレンド」を抽出してくれると勘違いしてしまっていた。なんともお恥ずかしい。でもまあ、スプライン近似の定義からありえないのだが（笑）

以下の図は元のデータがギザギザな部分を含む場合。

グリッドの数は400とまあまあ細かい。
これに対してspline.f90の中の3次スプライン補間法を適用。
するとどうだろう。まあ当たり前だけど、ギザギザの部分も見事に再現してくれた。
でも、自分としてはこのギザギザを削ぎ落とした滑らかな曲線がほしいわけだ。まあ考えられる方法はあるのでいいのだ。とりあえず問題発見はできた。これがプログラムの結果が上手くいかない１つの原因と思われる。

あと、もう１つは、外挿の精度の低さ。これも図の右端でわかる。これもプログラムが上手く結果を出さない１つの原因だと思われる。

さて、これからシャワー。今日はテニスだ！

【Fortran90】 3次スプライン補間法

関数を何らかの方法で補間する必要があるので、その方法として3次スプライン補間法を考える。よく使われているらしいから。

でも、これがよく使われる理由とは？データがたくさんあるときでも、ラグランジュ補間のように当てはまりの精度が問題にならないこと、なのかな。なぜ3次かはわからない。

B-スプラインというのも耳にしたことあるけど、これはデータの点を通過しないという点で、個人的には、問題がある。

スプラインのサブルーチンはFortran Source CodesのSPLINEがある。

3次スプライン補間法（Piecewise Cubic Spline）を使うには、次のサブルーチンを使う。

• spline_cubic_val(n,t,y,ypp,tval,yval,ypval,yppval)


• spline_cubic_set(n,t,y,ibcbeg,ybcbeg,ibcend,ybcend,ypp)
• 2階微分の情報を与えるサブルーチン。
• スプライン補間の境界条件はここで設定する。自然境界(natural boundary。2つの境界点における近似関数の2階微分が0となる)を使う場合、ibcbeg=2, ybcbeg=0.0D0, ibcend=2, ybcend=0.0D0とする必要がある。

Codingの手順（自然境界条件の場合）

1. n, t, yにそれぞれデータ・ポイントの数（整数）、argumentの配列（倍精度、n個の要素）、関数値の配列（倍精度、n個の要素）を設定する。

2. ``spline_cubic_set''サブルーチンをコールして2階微分の配列yppを得る。

3. ``spline_cubic_val''サブルーチンをコールして、argumentがtval（倍精度、input）の時の関数の近似値yval（倍精度、output）、1階微分ypval（倍精度、output）、2階微分yppval（倍精度、output）を得る。


Note 1: 微分可能でない関数の近似精度は低い。特に、微分可能でない点の近傍においてそれは顕著。

非線形方程式の解法：特殊なケース

　今、次のような2変数の非線形方程式を数値計算で解くことを考えているとする：

　f_1(x_1, x_2) = 0
　f_2(x_1, x_2) = 0.

さらに、変数がとりうる値の範囲が次のようになっているとする：

　l_1 ≦ x_1 ≦ u_1
　l_2 ≦ x_2 ≦ u_2.

　内点解が保証され、かつf_1とf_2が連続微分可能であれば、Newton-Raphson法を修正して問題を解くことができるはず。

　では、もし、微分可能性が満たされない場合はどうだろうか。この場合はSecant methodを使うことが考えられる。さらに、特殊な場合として、区間[l_2,u_2]が非常に小さいとする。
　このとき、2変数のシステムにsecant methodを適用するのではなく、x_1についてだけSecant methodを適用し、x_2については単純に加重平均でアップデートするという方法が考えられる。f_1とf_2が連続であれば、根の付近においてはそのような近似が利用できると考えられる。

非線形方程式を解く場合は次元を下げて

　''f(x) = 0''を満たすようなn次元ユークリッド空間内の点xを求める方法について、Newton-Raphson methodやMultivariate Secant methodなどの方法がある。
　次元nが大きくなればなるほど、計算時間が掛かるだけでなく、iterationが収束しにくいということが予想される。後者について、システムの均衡条件を使えば実は次元が2になるのに、それに気が付かずに次元が4のままで計算してしまうということもあったりする。特に、次元が大きいままでは収束しなかったものが、次元を低くすれば収束したというような場合ではかなり大きな違いが出てくる。均衡の条件を使うことでより効率的に数値計算を行うことができる可能性が示唆される。

連立方程式を解く

　連立方程式を数値計算で解く方法にはいくつかある。写像が（連続）微分可能か否かで分類することができる。

　○微分可能な場合
　・Newton-Raphson method

　○微分可能でない場合（※もちろん、微分可能な場合にも適用できる）
　・Bisection search
　・Secant method

いずれも一変数である場合は利用は単純。ただ、変数が複数ある場合、Bisection searchはよくわからん。Secant methodは複数あるときでも何らかの仮定の下では比較的簡単に実践できると思われる。Newton-Raphson methodについてとにかく偏微分してしまえばいいので、特に複雑なことをするわけでもない。