4.2.2.26 Fitting with an Ordinary Differential Equation

1 Summary
2 What you will learn
3 Example and Steps

User-defined Fitting Functions

Summary

In this tutorial we will show you how to define an ordinary differential equation (ODE) in the Fitting function Builder dialog and perform a fit of the data using this fitting function.

Minimum Origin Version Required: Origin 9.1 SR0

What you will learn

This tutorial will show you how to:

Define an ODE fitting function.
Call NAG functions using Origin C code.
Recalculate the ODE result only when parameters are updated.
Interpolate on the ODE result.

Example and Steps

In this tutorial, we will use a first order ordinary differential equation as an example:

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

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

where a is a parameter in the ordinary differential equation and y0 is the initial value for the ODE. NAG functions d02pvc and d02pcc are called using the Runge–Kutta method to solve the ODE problem.

Import Data

Start with an empty workbook. Select Help: Open Folder: Sample Folder... to open the "Samples" folder. In this folder, open the Curve Fitting subfolder and find the file Exponential Growth.dat. Drag-and-drop this file into the empty worksheet to import it.
Highlight B column, and select Plot: Symbol: Scatter from Origin menu. The graph should look like:

Define Fitting Function

The fitting function can be defined using the Fitting Function Builder tool.

Select Tools: Fitting Function Builder from Origin menu.
On the Fitting Function Builder dialog's Goal page, ensure the Create New Function radio button is selected. Click the Next button.
On the Name and Type page, From the Select or create a Category drop-down list select User Defined. For Function Name type FitODE. In the Function Type section select the Origin C radio button. Click Next.
On the Variables and Parameters page, in the Parameters field type a, y0. Click Next.

On the Origin C Fitting Function page, click the

button on the right of the Function Body edit box and define the function in Code Builder as follows. Include header files for NAG and fitting:

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

Define a static function to solve ODE by calling NAG functions. Call NAG function d02pvc to establish ODE model and d02pcc to solve the model:

struct user   // parameter in the ODE
{
	double a;
};

//Define the differentiate equation: y'=a*y
static void NAG_CALL f(double t, Integer neq, const double *y, double *yp, Nag_Comm *comm)
{
	neq; //Number of ordinary differential equations
	t; //Independent variable
	y; //Dependent variables
	yp; //First derivative

	struct user *sp = (struct user *)(comm->p);
	double a;
	a = sp->a;

	yp[0] = a*y[0];
}


//Use Runge-Kutta ODE23 to solve ODE
static bool nag_ode_fit( const double a, const double y0, const double tstart, 
  const double tend, const int nout, vector &vP )
{
	//nout: Number of points to output	
	if( nout < 2 )
		return false;

	vP.SetSize( nout );
	vP[0] = y0;

	int neq = 1; //Number of ordinary differential equations
	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;
	//Setup 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;
		//Solve 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;
}

Define the fitting function body _nlsfFitODE:

NLFitContext *pCtxt = Project.GetNLFitContext();
if ( pCtxt )
{
  static vector vX, vY;
  static int nSize;

  BOOL bIsNewParamValues = pCtxt->IsNewParamValues();
  // If parameters were updated, we will recalculate the ODE result.
  if ( bIsNewParamValues )
  {
    //Initial and final values of the independent variable
    double tstart = 0.02, tend = 2, tinc;
    int nout = 100; //Number of points
		
    tinc = (tend-tstart)/(nout-1);
    vX.Data( tstart, tend, tinc );
    nSize = vX.GetSize();
		
    if( !nag_ode_fit( a, y0, tstart, tend, nout, vY) )
      return;
  }
	
  //Interpolate y from fitting data's x on the ODE result.
  ocmath_interpolate( &x, &y, 1, vX, vY, nSize );
}

Click the Compile button to compile the function body. To return, click Return to Dialog button. The Evaluate button on the dialog has a black figure of a running person. Clicking on it will show y=2.6627270424371. This implies that the function works. Click Next.

In the Parameter Initialization Code page, select Use Custom Code radio button, and type the code in the Initialization Code box:
```
//Set the start y value in fitting data as the initial value of y0
y0 = y_data[0];
a = 1;
```
and click Finish.

Note: In order to monitor the the fitted parameters, NLFitContext class was introduced in defining fitting function to achieve some key information within the fitter

Fit the Curve

Highlight B column, then select Analysis: Fitting: Nonlinear Curve Fit from Origin menu. In the NLFit dialog, select Settings: Function Selection, in the page select User Defined from the Category drop-down list and FitODE function from the Function drop-down list.
Select Parameters tab, and fix y0 as shown in the dialog.
Click Fit button to fit the curve.