# 3.9.5.2 Ordinary Differential Equations

## Contents

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>
////////////////////////////////////////////////////////////////////////////////////

#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 = y;
yp = -y;
}

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 = ZERO; // initial values
ystart = 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;
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, ystart);

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 will then be the solution we want at time twant
printf("%8.3f   %8.3f   %8.3f\n", tgot, ygot, ygot);
}

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>
////////////////////////////////////////////////////////////////////////////////////

#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 = y;
yp = -y-beta*y;
}

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 = ZERO; // initial values
ystart = 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;
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, ystart);

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 will then be the solution we want at time twant
printf("%8.3f   %8.3f   %8.3f\n", tgot, ygot, ygot);

}

free(rwsav);
}

return 1;
}