3.9.5.2 Ordinary Differential Equations

User can access NAG Ordinary Differential Equation (ODE) solvers to evaluate the intial value problem. NAG supplies flexible ODE solvers and methods to solve this kind of problem.

Version Info

Minimum Origin Version Required: Origin 8 SR1

A simple second order ODE

#include <Origin.h>
////////////////////////////////////////////////////////////////////////////////////

#include <OC_nag.h>
////////////////////////////////////////////////////////////////////////////////////
// Start your functions here.

#define NEQ 2
#define ZERO 0.0
#define ONE 1.0
#define TWO 2.0
#define FOUR 4.0


// Supply your own differential equation system
// Here is an example of second order oscillator system
// For higher order differential equation, you need to supply its equivalent first-order differential equation system
// i.e.  y'' = -y with initial value y(0)=0, y'(0)=1 can be rewriten as y1'= y2; y2'= -y1
// In addition, by Nag_User struct you can also pass in the parameters for the system
static void NAG_CALL f(double t, Integer neq, const double *y, double *yp, Nag_Comm *comm)
{
	yp[0] = y[1];
	yp[1] = -y[0];
}


void test_diff_eqn()
{
	Integer neq;
	Nag_RK_method method;
	double hstart, pii, tgot, tend, tinc;
	double  tol, tstart, twant;
	Integer i, j, nout;
	double thres[NEQ], ygot[NEQ], ymax[NEQ], ypgot[NEQ], ystart[NEQ];
	Nag_ErrorAssess errass;
	Nag_Comm comm;

	/* Set initial conditions and input for d02pvc */
	neq = NEQ;  // number of differential equations

	pii = 3.1415926;
	tstart = ZERO; //  span of the time interval (the independant variable) for the solution
	tend = TWO*pii;

	ystart[0] = ZERO; // initial values
	ystart[1] = ONE;

	for (i=0; i<neq; i++)
		thres[i] = 1.0e-8;
	errass = Nag_ErrorAssess_off;
	hstart = ZERO;
	method  = Nag_RK_2_3;

	/* Set control for output */
	nout = 8;
	tinc = (tend-tstart)/nout;

	for (i=1; i<=2; i++)
	{
		if (i==1)    // solve for two kinds of accuracy for the solution
			tol = 1.0e-3;
		else
			tol = 1.0e-4;
		// setup function need to be called to call the solver nag_ode_ivp_rkts_range
		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, NAGERR_DEFAULT);

		printf("\nCalculation with tol = %8.1e\n\n",tol);
		printf("%8.3f   %8.3f   %8.3f\n", tstart, ystart[0], ystart[1]);

		for (j=nout-1; j>=0; j--)
		{
			twant = tend - j*tinc; // the solution at twant is wanted

			nag_ode_ivp_rkts_range(f, neq, twant, &tgot, ygot, ypgot, ymax, &comm, iwsav, rwsav, NAGERR_DEFAULT);
			// obtain the solution at tgot, and when no error exit tgot should equal to twant
			// and ygot[0] will then be the solution we want at time twant
			printf("%8.3f   %8.3f   %8.3f\n", tgot, ygot[0], ygot[1]);   
		}

		free(rwsav);
	}
}

Second order ODE with Parameters

The following example defines a differential equation system with a damped parameter.

#include <Origin.h>
////////////////////////////////////////////////////////////////////////////////////

#include <OC_nag.h>
////////////////////////////////////////////////////////////////////////////////////
// Start your functions here.

#define NEQ 2
#define ZERO 0.0
#define ONE 1.0
#define TWO 2.0
#define FOUR 4.0



struct user   // parameters in the differential equations
{
	double beta;
 
};

// Supply your own differential equation system here
// Here is an example of damped oscillator system
// The equivalent first-order differential equation system for this damped system
// y'' = -y-beta*y' with initial value y(0)=0, y'(0)=1 is y1'= y2; y2'= -y1-beta*y2
// By Nag_User struct we can passed in the parameter beta

static void NAG_CALL f(double t, Integer neq, const double *y, double *yp, Nag_Comm *comm)
{
	struct user *sp = (struct user *)(comm->p);

	double beta; // temp variable
	beta = sp->beta;

	yp[0] = y[1];
	yp[1] = -y[0]-beta*y[1];
}


int test_diff_eqn_damp_osc()
{
	Integer neq;
	Nag_RK_method method;
	double hstart, pii, tgot, tend, tinc;
	double  tol, tstart, twant;
	Integer i, j, nout;
	double thres[NEQ], ygot[NEQ], ymax[NEQ], ypgot[NEQ], ystart[NEQ];
	Nag_ErrorAssess errass;

	Nag_Comm comm;  
	struct user s;
	double beta = 0.02;
	s.beta = beta; // the parameters of the differential equation systems is passed through this Nag_User communication struct
	comm.p = (Pointer)&s;  // you can define your own parametrized differential equation system through this way

	/* Set initial conditions and input for d02pvc */
	neq = NEQ;  // number of differential equations

	pii = 3.1415926;
	tstart = ZERO; //  span of the time interval (the independant variable) for the solution
	tend = TWO*pii;

	ystart[0] = ZERO; // initial values
	ystart[1] = ONE;


	for (i=0; i<neq; i++)
		thres[i] = 1.0e-8;
	errass = Nag_ErrorAssess_off;
	hstart = ZERO;
	method  = Nag_RK_2_3;

	/* Set control for output */
	nout = 8;
	tinc = (tend-tstart)/nout;


	for (i=1; i<=2; i++)
	{
		if (i==1)    // solve for two kinds of accuracy for the solution
			tol = 1.0e-3;
		else
			tol = 1.0e-4;
		// setup function need to be called to call the solver nag_ode_ivp_rkts_range
		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, NAGERR_DEFAULT);

		printf("\nCalculation with tol = %8.1e\n\n",tol);
		printf("%8.3f   %8.3f   %8.3f\n", tstart, ystart[0], ystart[1]);

		for (j=nout-1; j>=0; j--)
		{
			twant = tend - j*tinc; // the solution at twant is wanted

			nag_ode_ivp_rkts_range(f, neq, twant, &tgot, ygot, ypgot, ymax, &comm, iwsav, rwsav, NAGERR_DEFAULT);
			// obtain the solution at tgot, and when no error exit tgot should equal to twant
			// and ygot[0] will then be the solution we want at time twant
			printf("%8.3f   %8.3f   %8.3f\n", tgot, ygot[0], ygot[1]);   
		  
		}

		free(rwsav);
	}
	
	return 1;
}