/* nag_ode_ivp_rk_step_revcomm (d02pgc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.6, 2022.
*/
#include <math.h>
#include <nag.h>
int main(void) {
/* Constants */
double const tol = 1.0e-5;
Integer const n = 2;
Integer const liwsav = 130;
/* Scalars */
double hnext, hstart, t, t1, t2, tend, tnow, tout, tprev, waste;
Integer ind, irevcm, j, k, nchange, stepcost, stepsok, totf, lrwsav, lwcomm,
exit_status = 0;
/* Arrays */
double c[17];
double *rwsav = 0, *thresh = 0, *troot = 0, *wcomm = 0, *y = 0, *ynow = 0,
*yout = 0, *yp = 0, *ypnow = 0, *yprev = 0;
Integer *iroot = 0, *iwsav = 0;
char nag_enum_arg[40];
/* Nag Types */
Nag_Boolean icheck;
NagError fail, fail2;
Nag_RK_method method;
Nag_ErrorAssess errass;
INIT_FAIL(fail);
INIT_FAIL(fail2);
printf("nag_ode_ivp_rk_step_revcomm (d02pgc) Example Program Results\n\n");
lrwsav = 350 + 32 * n;
lwcomm = 6 * n;
if (!(thresh = NAG_ALLOC((n), double)) ||
!(iwsav = NAG_ALLOC((liwsav), Integer)) ||
!(rwsav = NAG_ALLOC((lrwsav), double)) ||
!(ynow = NAG_ALLOC((n), double)) || !(ypnow = NAG_ALLOC((n), double)) ||
!(yprev = NAG_ALLOC((n), double)) ||
!(wcomm = NAG_ALLOC((lwcomm), double)) ||
!(yout = NAG_ALLOC((n), double)) || !(iroot = NAG_ALLOC((n), Integer)) ||
!(y = NAG_ALLOC((n), double)) || !(yp = NAG_ALLOC((n), double)) ||
!(troot = NAG_ALLOC((n), double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Set initial conditions for ODE and parameters for the integrator. */
scanf(" %39s%*[^\n] ", nag_enum_arg);
/* nag_enum_name_to_value (x04nac) Converts NAG enum member name to value. */
method = (Nag_RK_method)nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s%*[^\n] ", nag_enum_arg);
errass = (Nag_ErrorAssess)nag_enum_name_to_value(nag_enum_arg);
scanf("%lf%lf%*[^\n] ", &t, &tend);
for (j = 0; j < n; j++)
scanf("%lf", &ynow[j]);
scanf("%*[^\n] ");
scanf("%lf%*[^\n] ", &hstart);
for (j = 0; j < n; j++)
scanf("%lf", &thresh[j]);
scanf("%*[^\n] ");
/* Initialize Runge-Kutta method for integrating ODE using
* nag_ode_ivp_rkts_setup (d02pqc).
*/
nag_ode_ivp_rkts_setup(n, t, tend, ynow, tol, thresh, method, errass, hstart,
iwsav, rwsav, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_ivp_rkts_setup (d02pqc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf(" Calculation with tol = %8.1e\n", tol);
printf(" t y1 y2\n");
printf("%7.3f", t);
for (k = 0; k < n; k++)
printf("%11.4f", ynow[k]);
printf("\n");
tout = 0.1;
tnow = t;
while (tnow < tend) {
tprev = tnow;
for (k = 0; k < n; ++k)
yprev[k] = ynow[k];
/* Solve ODE by Runge-Kutta method by a sequence of single steps.
* Each step requires a reverse communication loop around
* nag_ode_ivp_rk_step_revcomm (d02pgc).
*/
irevcm = 0;
while (irevcm >= 0) {
nag_ode_ivp_rk_step_revcomm(&irevcm, n, &tnow, ynow, ypnow, iwsav, rwsav,
&fail);
if (irevcm > 0) {
ypnow[0] = ynow[1];
ypnow[1] = -ynow[0];
}
}
if (irevcm == -2) {
if (fail.code != NE_RK_POINTS && fail.code != NE_STIFF_PROBLEM &&
fail.code != NW_RK_TOO_MANY) {
printf("Error from nag_ode_ivp_rk_step (d02pgc).\n%s\n", fail.message);
exit_status = 2;
goto END;
}
}
/* Detect sign changes in last step */
for (k = 0; k < n; ++k)
iroot[k] = 0;
nchange = 0;
for (k = 0; k < n; ++k) {
if (ynow[k] * yprev[k] < 0.0) {
iroot[nchange] = k + 1;
nchange++;
}
}
if (tnow >= tout || nchange > 0) {
/* nag_ode_ivp_rk_interp_setup (d02phc).
* Compute interpolant for the last step taken by the Runge-Kutta
* integrator nag_ode_ivp_rk_step_revcomm (d02pgc).
*/
irevcm = 0;
while (irevcm >= 0) {
nag_ode_ivp_rk_interp_setup(&irevcm, n, n, &t, y, yp, wcomm, lwcomm,
iwsav, rwsav, &fail);
if (irevcm > 0) {
yp[0] = y[1];
yp[1] = -y[0];
}
}
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_ivp_rk_interp_setup (d02phc).\n%s\n",
fail.message);
exit_status = 3;
goto END;
}
icheck = Nag_TRUE;
for (k = 0; k < nchange; ++k) {
j = iroot[k] - 1;
t1 = tprev;
t2 = tnow;
ind = 1;
/* nag_roots_contfn_brent_rcomm (c05azc).
* Locates a simple zero of a continuous function.
* Reverse communication.
*/
while (ind != 0) {
nag_roots_contfn_brent_rcomm(&t1, &t2, y[j], tol, Nag_Mixed, c, &ind,
&fail);
if (ind > 1) {
/* nag_ode_ivp_rk_interp_eval (d02pjc).
* Evaluate interpolant at a point in the last integrated step
* as computed by nag_ode_ivp_rk_interp_setup (d02phc).
*/
nag_ode_ivp_rk_interp_eval(icheck, n, n, t1, 0, y, wcomm, lwcomm,
iwsav, rwsav, &fail2);
icheck = Nag_FALSE;
}
}
if (fail.code != NE_NOERROR) {
printf("Error from nag_roots_contfn_brent_rcomm (c05azc).\n%s\n",
fail.message);
exit_status = 4;
goto END;
}
troot[k] = t1;
}
while (tnow >= tout) {
for (k = 0; k < nchange; k++) {
if (troot[k] < tout && iroot[k] > 0) {
printf("Component %2" NAG_IFMT " has a root at t = %7.4f\n",
iroot[k], troot[k]);
iroot[k] = -iroot[k];
}
}
nag_ode_ivp_rk_interp_eval(icheck, n, n, tout, 0, yout, wcomm, lwcomm,
iwsav, rwsav, &fail2);
icheck = Nag_FALSE;
printf("%7.3f", tout);
for (k = 0; k < n; k++) {
printf("%11.4f", yout[k]);
}
printf("\n");
tout = tout + 0.1;
}
for (k = 0; k < nchange; k++) {
if (iroot[k] > 0) {
printf("Component %2" NAG_IFMT " has a root at t = %7.4f\n", iroot[k],
troot[k]);
}
}
}
}
/* Get diagnostics on whole integration using
* nag_ode_ivp_rkts_diag (d02ptc).
*/
nag_ode_ivp_rkts_diag(&totf, &stepcost, &waste, &stepsok, &hnext, iwsav,
rwsav, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_ivp_rkts_diag (d02ptc).\n%s\n", fail.message);
exit_status = 5;
goto END;
}
printf("\nCost of the integration in evaluations of f is %6" NAG_IFMT "\n\n",
totf);
END:
NAG_FREE(thresh);
NAG_FREE(ynow);
NAG_FREE(ypnow);
NAG_FREE(yprev);
NAG_FREE(yout);
NAG_FREE(y);
NAG_FREE(yp);
NAG_FREE(wcomm);
NAG_FREE(rwsav);
NAG_FREE(iwsav);
NAG_FREE(iroot);
NAG_FREE(troot);
return exit_status;
}