/* D02PE_P0W_F C++ Header Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.2, 2021.
*/
#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
using namespace std;
#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL f(void *& ad_handle,
const double & t,
const Integer &n,
const double y[],
double yp[],
Integer iuser[],
double ruser[]);
#ifdef __cplusplus
}
#endif
int main(void)
{
const Integer n = 2, npts = 8;
const Integer liwsav = 130;
const Integer lrwsav = 350 + 32 * n;
Integer exit_status = 0;
double * rwsav = 0, *thresh = 0, *ygot = 0, *ymax = 0;
double * ypgot = 0, *y = 0, ruser[1];
Integer *iwsav = 0, iuser[1];
cout << "D02PE_P0W_F C++ Header Example Program Results\n\n";
thresh = new double[n];
ygot = new double[n];
y = new double[n];
ypgot = new double[n];
ymax = new double[n];
iwsav = new Integer[liwsav];
rwsav = new double[lrwsav];
// Set initial conditions for ODE and parameters for the integrator.
Integer method = 1;
double tol, hstart, tend, tstart;
tstart = 0.0;
tol = 1.0e-4;
tend = 2.0 * nag_math_pi;
y[0] = 0.0;
y[1] = 1.0;
hstart = 0.0;
thresh[0] = 1.0e-8;
thresh[1] = 1.0e-8;
{
cout << "\n Calculation with tol = " << tol << endl;
}
cout.setf(ios::fixed);
cout.setf(ios::right);
cout.precision(3);
cout << "\n t y1 y2" << endl;
cout.width(6);
cout << tstart;
for (int k = 0; k < n; k++)
{
cout.width(10);
cout << y[k];
}
cout << endl;
// Set control for output
double tinc = (tend - tstart) / (double)(npts);
Integer ifail = 0;
void * ad_handle = 0;
// Initialize Runge-Kutta method for integrating ODE
ifail = 0;
nag::ad::d02pq(ad_handle, n, tstart, tend, y, tol, thresh, method, hstart,
iwsav, rwsav, ifail);
double tgot, twant;
twant = tstart;
for (int j = 0; j < npts; ++j)
{
twant = twant + tinc;
ifail = 0;
nag::ad::d02pe(ad_handle, f, n, twant, tgot, ygot, ypgot, ymax, -1, iuser,
-1, ruser, iwsav, rwsav, ifail);
cout.width(6);
cout << tgot;
for (int k = 0; k < n; ++k)
{
cout.width(10);
cout << ygot[k];
}
cout << endl;
}
double hnext, waste;
Integer fevals, stepcost, stepsok;
ifail = 0;
nag::ad::d02pt(ad_handle, fevals, stepcost, waste, stepsok, hnext, iwsav,
rwsav, ifail);
cout << "\n Cost of the integration in evaluations of f is " << fevals;
cout << endl;
delete[] thresh;
delete[] ygot;
delete[] y;
delete[] ypgot;
delete[] ymax;
delete[] iwsav;
delete[] rwsav;
return exit_status;
}
static void NAG_CALL f(void *& ad_handle,
const double & t,
const Integer &n,
const double y[],
double yp[],
Integer iuser[],
double ruser[])
{
yp[0] = y[1];
yp[1] = -y[0];
}