数値積分を用いた未知係数の最適化を行うプログ

数値積分を用いた未知係数の最適化を行うプログラムを教えてください。数値積分を用いた未知係数の最適化を行うプログラムを教えてください。 先日、ルンゲクッタ法による微分方程式の数値解析の方法を教えていただいたのですが、今度は実測値とルンゲクッタ法で計算する曲線がフィットするように、未知の係数を最小二乗法で決定したいのですが、どのようなアルゴリズムを書けばいいのですが。Excel2003のマクロで書きたいと思いますので、お願いします。質問の内容が不明瞭であったことを、まずお詫びさせていただきます。ある時刻における、ある2種類の物質の濃度C1,C2を求める反応速度の微分方程式は解くことができずdC1/dt=F(C1,C2),dC2/dt=G(C1,C2)の形で与えられています。ここからRunge-Kutta法により数値的に解くことには成功しましたが、ほかの方法で任意の時刻における実測の濃度を決め、それに合うように速度定数k1,k2の値を求める方法を探しています。速度定数とは化学物質の濃度C1,C2の時間変化について

__dC1/dt_=_v1_=_k1*(C1^a11)*(C2^a12)______①

__dC2/dt_=_v2_=_k2*(C1^a21)*(C2^a22)______②

における、k1,_k2_を指すものと考えます。これらを

求めるためのデータとしてC1,_C2が共存する状態で

時間tに対するC1,_C2_の測定値C1(t),C2(t)_が得られ

ている場合に、最小二乗法を利用するものと理解します

 

結論に近い方から述べます。

式①、②において、k1,a11,a12,_k2,a21,a22が未知数で

v1(t),_C1(t),_v2(t),_C2(t)が時刻t_について既知とし

t_=_t0,_t1,_..._,_tnについて記録されていると

①、②のlogをとって

__log(v1)_=_log(k1)+a11*log(C1)+a12*log(C2)___③

__log(v2)_=_log(k2)+a21*log(C1)+a22*log(C2)___④

を得、K1=log(k1),_K2=log(k2),_a11,_a12,_a21,_a22

について、最小二乗解を求めることが出来ます。

ここで、v1,v2_については、時刻歴C1(t),C2(t)から

求めることになりますが、C1(t),_C2(t)は測定値で、

時刻に対して独立なランダム誤差によるバラツキを

含んでいるので

__v1(t)_=_{C1(t+dt)-C1(t)}/dt

で直接求めた場合、誤差の影響が拡大されて、適当で

ありません。平滑化された変化速度dC1/dtやdC2/dtを

得るためには、C1(t),_C2(t)の時間的平滑化が必要で

す。それには、これらの測定値に対して

多項式

__C1(t)_=_c10_+_c11*t_+_c12*t^2_+_・・・__⑤

__C2(t)_=_c20_+_c21*t_+_c22*t^2_+_・・・__⑥

で、代表させる方法が考えられます。そのためには

c10,_c11,_c12,_..._および、c20,_c21,_c22,..

を最小二乗法で求めることができます。ここで多項式

の次数はどうするかですが、これは低すぎれば充分に

傾向を表示できません。逆に高すぎれば測定値の細動

に敏感に対応して、曲線に振動傾向を持ち込む可能性

があります。

 

以上の考察から、この問題については2段階の計算を

要し、最初の計算で多項式⑤、⑥を定めるには

指定次数の代数多項式のための最小二乗法が必要で、

この結果を用いてv1(t),_v2(t)の時刻歴を求めた後

③、④の最小二乗法により、k1,_k2を求めることに

なると思われます。2種類の最小二乗法はかなり様子

が異なるので、共通の部分、連立方程式の解法とか

最小二乗法の正規方程式(連立方程式)組立てとか

は同じソースを用いるが、実行ファイルは別物として

扱う方が楽勝と考えられます。

 

以上の理由で、まず可変次数の代数多項式を求める

最小二乗法を提示しようと考えます。

 

ご意見があれば、別途、質問の形でお出し下さい。

 

(以上)市販の微分方程式の数値解析の本の中でも、実際にどの方法を適用すればいいか分からなかったので、非常に助かりました。