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 と表示されます。関数が動作することを意味します。進むをクリックします。