/* nag::ad::d02pf Tangent over Tangent Example Program.
*/
#include <dco.hpp>
#include <iostream>
#include <nagad.h>
// Function which calls NAG AD routines.
template <typename T> void func(std::vector<T> &r, std::vector<T> &y);
// Driver with the tangent calls.
// Solves the equation y''=-y, y0(=y(0))=0, y0'(=y'(0))=1 by solving the ODE
// system:
// y1'=r1*y2, y2'=-r2*y1
// over the range [0,2*pi] with initial conditions y1=0, y2=1.
// Also computes the sum of all Hessian elements dy2/dr2.
void driver(const std::vector<double> &rv,
std::vector<double> & yv,
double & d2y_dr2);
int main()
{
std::cout
<< "nag::ad::d02pf Tangent over Tangent Example Program Results\n\n";
// Set problem parameters
std::vector<double> rv{1.0, 1.0};
// Solution y
std::vector<double> yv{0.0, 1.0};
// Sum of Hessian elements
double d2y_dr2;
// Call driver
driver(rv, yv, d2y_dr2);
// Print outputs
std::cout << "\n Derivatives calculated: Second order tangents\n";
std::cout << " Computational mode : algorithmic\n";
std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(6);
std::cout << "\n Solution computed with required tolerance " << 1e-4
<< std::endl;
for (std::size_t i = 0; i < yv.size(); i++)
{
std::cout << " y" << i + 1 << " = " << yv[i] << std::endl;
}
// Print derivatives
std::cout
<< "\n Sum of all Hessian elements of solution vector y w.r.t. problem parameters :\n";
std::cout << " sum_ijk [d2y/dr2]_ijk = " << d2y_dr2 << std::endl;
return 0;
}
// Driver with the tangent calls.
// Solves the equation y''=-y, y0(=y(0))=0, y0'(=y'(0))=1 by solving the ODE
// system:
// y1'=r1*y2, y2'=-r2*y1
// over the range [0,2*pi] with initial conditions y1=0, y2=1.
// Also computes the sum of all Hessian elements dy2/dr2.
void driver(const std::vector<double> &rv,
std::vector<double> & yv,
double & d2y_dr2)
{
using mode = dco::gt1s<dco::gt1s<double>::type>;
using T = mode::type;
// Size of system of ODEs
Integer n = rv.size();
// Variables to differentiate w.r.t.
std::vector<T> r(begin(rv), end(rv));
dco::derivative(dco::value(r)) = std::vector<double>(n, 1.0);
dco::value(dco::derivative(r)) = std::vector<double>(n, 1.0);
// Solution y, variable to differentiate
std::vector<T> y(begin(yv), end(yv));
// Call the NAG AD Lib functions
func(r, y);
// Extract the computed solutions
yv = dco::passive_value(y);
// Sum the tangents of r
d2y_dr2 = 0;
for (T const &yi : y)
{
d2y_dr2 += dco::derivative(dco::derivative(yi));
}
}
// function which calls NAG AD Library routines
template <typename T> void func(std::vector<T> &r, std::vector<T> &y)
{
// Active variables
const Integer n = y.size();
const Integer liwsav = 130, lrwsav = 350 + 32 * n;
std::vector<T> thresh(n, 1e-8), ypnow(n), rwsav(lrwsav);
std::vector<Integer> iwsav(liwsav);
// Set parameters for the integrator.
Integer method = 2;
T tol = 1e-4, hstart = 0.0, tend = 2.0 * nag_math_pi, tstart = 0.0;
// Create AD configuration data object
Integer ifail = 0;
nag::ad::handle_t ad_handle;
// Initialize Runge-Kutta method for integrating ODE
ifail = 0;
nag::ad::d02pq(ad_handle, n, tstart, tend, y.data(), tol, thresh.data(),
method, hstart, iwsav.data(), rwsav.data(), ifail);
T tnow = tstart;
auto f = [&](nag::ad::handle_t &ad_handle,
const T & t,
const Integer & n,
const T y[],
T yp[])
{
yp[0] = r[0] * y[1];
yp[1] = -r[1] * y[0];
};
do
{
ifail = 0;
// Solve an initial value problem for a 1st order system of ODEs
nag::ad::d02pf(ad_handle, f, n, tnow, y.data(), ypnow.data(), iwsav.data(), rwsav.data(), ifail);
} while (tnow < tend);
}