常微分方程式によるフィット

概要

このチュートリアルでは、フィット関数ビルダーを使用して常微分方程式（ODE）を作成し、この関数を使用してデータのフィットを実行する方法をご紹介します。

必要なOriginのバージョン:Origin 9.1 SR0以降

学習する項目

このチュートリアルでは、以下の項目について解説します。

ODFフィット関数を定義する
OriginCコードを使用してNAG関数を呼び出す
パラメータが更新されたときのみODE結果を再計算する
ODF結果の補間

3 サンプルと操作

このチュートリアルでは、以下のサンプルのような、1階常微分方程式を使用します。

$\frac{\mathrm{d} y}{\mathrm{d} x}=ay$

$y|_{x=x0}=y0$

ここで、aは常微分方程式のパラメータで、y0はODEの初期値です。このODEの問題を解くために、Runge–Kuttaメソッドを使用して、NAG関数 d02pvc と d02pcc が呼び出されます。

データのインポート

新しいワークブックを用意します。ヘルプ: フォルダを開く: サンプルフォルダを選択して、サンプルフォルダを開きます。このフォルダ内のCurve FittingサブフォルダにあるExponential Growth.dat ファイルを探します。空のワークシートにファイルをドラッグアンドドロップしてインポートします。
B列を選択して、Originのメニューから作図：シンボル図：散布図を選択します。グラフは次のようになります。

フィット関数を定義する

フィット関数は、フィット関数ビルダーを使用して定義できます。

メニューから、ツール：フィット関数ビルダーを選択します。
開いたフィット関数ビルダーの処理のゴールページで、新しい関数の作成が選択されているの確認します。進むボタンをクリックします。
関数名と関数形式のページでは、関数カテゴリーの選択ドロップダウンリストでUser Definedを選択し、関数名をFitODEとします。関数形式は、Origin Cにします。進むボタンをクリックします。進むをクリックします。
変数とパラメータページでは、パラメータとして、a, y0 と入力します。進むをクリックします。

Origin Cフィット関数のページで、関数内容編集ボックスの右上にある、

ボタンをクリックしてコードビルダを開き、関数を以下のように定義します。 NAGとフィットのためにヘッダファイルを含めます。

#include <oc_nag.h>
#include <ONLSF.H>

NAG関数を呼び出して、ODEを解くための静的な関数を定義します。NAG関数d02pvcを呼び出して、ODEモデルを確立し、d02pccでモデルを解きます。

struct user   //ODEのパラメータ
{
	double a;
};
//微分方程式を定義: y'=a*y
static void NAG_CALL f(double t, Integer neq, const double *y, double *yp, Nag_Comm *comm)
{
	neq; //常微分方程式の数
	t; //独立変数
	y; //従属変数
	yp; //一階微分
	struct user *sp = (struct user *)(comm->p);
	double a;
	a = sp->a;
	yp[0] = a*y[0];
}

//ODEを解くためにRunge-Kutta ODE23を使用
static bool nag_ode_fit( const double a, const double y0, const double tstart, 
  const double tend, const int nout, vector &vP )
{
	//nout: 出力ポイント数	
	if( nout < 2 )
		return false;
	vP.SetSize( nout );
	vP[0] = y0;
	int neq = 1; //常微分方程式の数
	Nag_RK_method method;
	double hstart, tgot, tinc;
	double tol, twant;
	int i, j;
	vector thres(neq), ygot(neq), ymax(neq), ypgot(neq), ystart(neq);
		
	Nag_ErrorAssess errass;
	Nag_Comm comm;
	struct user s;
	s.a = a;
	comm.p = (Pointer)&s;
	ystart[0] = y0;
	for (i=0; i<neq; i++)
		thres[i] = 1.0e-8;
	errass = Nag_ErrorAssess_off;
	hstart = 0;
	method = Nag_RK_2_3;
	tinc = (tend-tstart)/(nout-1);
	tol = 1.0e-3;
	NagError nagErr1;
	//ODEセットアップ
	int iwsav[130];
	double* rwsav = (double*)malloc((350 + 32*neq)*sizeof(double));
	nag_ode_ivp_rkts_setup(neq, tstart, tend, ystart, tol, thres, method, errass, hstart, iwsav, rwsav, &nagErr1);
	if( nagErr1.code != NE_NOERROR )
	{
		free(rwsav);
		return false;
	}
	for (j=1; j<nout; j++)
	{
		twant = tstart + j*tinc;
		NagError nagErr2;
		//ODEを解く
		nag_ode_ivp_rkts_range(f, neq, twant, &tgot, ygot, ypgot, ymax, &comm, iwsav, rwsav, &nagErr2);
			
		if( nagErr2.code != NE_NOERROR )
		{
			free(rwsav);
			return false;
		}
			
		vP[j] = ygot[0];
	}
	free(rwsav);
	return true;
}

フィット関数内容_nlsfFitODEを定義

NLFitContext *pCtxt = Project.GetNLFitContext();
if ( pCtxt )
{
  static vector vX, vY;
  static int nSize;
  BOOL bIsNewParamValues = pCtxt->IsNewParamValues();
  // パラメータを更新したら、ODEの結果を再計算
  if ( bIsNewParamValues )
  {
    //独立変数の初期および最終の値
    double tstart = 0.02, tend = 2, tinc;
    int nout = 100; //ポイント数
		
    tinc = (tend-tstart)/(nout-1);
    vX.Data( tstart, tend, tinc );
    nSize = vX.GetSize();
		
    if( !nag_ode_fit( a, y0, tstart, tend, nout, vY) )
      return;
  }
	
  //ODE結果の近似データのxからyを内挿
  ocmath_interpolate( &x, &y, 1, vX, vY, nSize );
  }

コンパイルボタンをクリックして関数内容をコンパイルします。NLSFに戻るボタンをクリックし、戻ります。ダイアログの人が走っているマークのボタンは評価ボタンです。これをクリックすると、y=2.6627270424371 と表示されます。関数が動作することを意味します。進むをクリックします。

パラメータの初期化コードのページで、カスタムコードを使用するのラジオボタンを選択して、初期化ボックスのコードを入力します。
```
// y0の初期値として、フィットデータのyの開始値をセット
y0 = y_data[0];
a = 1;
```
完了ボタンをクリックします。

Note:NLFitContext classでフィットしたパラメータを確認できます。

曲線をフィットする

B列を選択して、Originメニューから、解析：フィット：非線形曲線フィットを選択します。開いたNLFitダイアログの、設定：関数選択のページで、カテゴリドロップダウンリストからUser Defined を選択し、関数ドロップダウンリストからFitODEを選択します。
パラメータタブを開き、下図のようにy0を固定します。
フィットボタンをクリックして、曲線をフィットします。

フィット結果

グラフは下図のようになります。

フィットパラメータは以下のようになります

パラメータ	値	標準誤差
a	1.30272	0.00858
y0	0.77038	0

Note：より複雑なODEフィット関数を使用したフィットも同様の方法で実行できます。