/* nag_dae_ivp_dassl_gen (d02nec) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 9, 2009.
*
*/
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd02.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL res(Integer neq, double t, const double y[],
const double ydot[], double r[], Integer *ires,
Nag_Comm *comm);
static void NAG_CALL jac(Integer neq, double t, const double y[],
const double ydot[], double *pd, double cj,
Nag_Comm *comm);
static void NAG_CALL myjac(Integer neq, Integer ml, Integer mu, double t,
const double y[], const double ydot[], double *pd,
double cj);
#ifdef __cplusplus
}
#endif
int main(void)
{
/*Integer scalar and array declarations */
Integer exit_status = 0, maxord = 5;
Nag_Comm comm;
Integer neq, licom, mu, ml, lcom;
Integer i, itask, j;
Nag_Boolean vector_tol;
Integer *icom = 0;
NagError fail;
/*Double scalar and array declarations */
double dt, h0, hmax, t, tout;
double *atol = 0, *com = 0, *rtol = 0, *y = 0, *ydot = 0;
static double ruser[2] = {-1.0, -1.0};
INIT_FAIL(fail);
printf("nag_dae_ivp_dassl_gen (d02nec) Example Program Results\n\n");
/* For communication with user-supplied functions: */
comm.user = ruser;
/* Set problem parameters required to allocate arrays */
neq = 3;
ml = 1;
mu = 2;
licom = 50+neq;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1);
if (
!(atol = NAG_ALLOC(neq, double)) ||
!(com = NAG_ALLOC(lcom, double)) ||
!(rtol = NAG_ALLOC(neq, double)) ||
!(y = NAG_ALLOC(neq, double)) ||
!(ydot = NAG_ALLOC(neq, double)) ||
!(comm.iuser = NAG_ALLOC(2, Integer)) ||
!(icom = NAG_ALLOC(licom, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialize the problem, specifying that the Jacobian is to be */
/* evaluated analytically using the provided routine jac. */
h0 = 0.0;
hmax = 0.0;
vector_tol = Nag_TRUE;
/*
* nag_dae_ivp_dassl_setup (d02mwc)
* Implicit DAE/ODEs, stiff ivp, setup for nag_dae_ivp_dassl_gen (d02nec)
*/
nag_dae_ivp_dassl_setup(neq, maxord, Nag_AnalyticalJacobian, hmax, h0,
vector_tol, icom, licom, com, lcom, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_dae_ivp_dassl_setup (d02mwc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Specify that the Jacobian is banded.
*
* nag_dae_ivp_dassl_linalg (d02npc)
* ODE/DAEs, ivp, linear algebra setup routine for
* nag_dae_ivp_dassl_gen (d02nec)
*/
nag_dae_ivp_dassl_linalg(neq, ml, mu, icom, licom, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_dae_ivp_dassl_linalg (d02npc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Set initial values*/
t = 0.00e0;
tout = 0.00e0;
dt = 0.020e0;
for (i = 0; i < neq; i++)
{
rtol[i] = 1.00e-3;
atol[i] = 1.00e-6;
y[i] = 0.00e0;
ydot[i] = 0.00e0;
}
y[0] = 1.00e0;
/* Use the comm.iuser array to pass the band dimensions through to jac. */
/* An alternative would be to hard code values for ml and mu in jac. */
comm.iuser[0] = ml;
comm.iuser[1] = mu;
printf(" t y(1) y(2) y(3)\n");
printf("%8.4f", t);
for (i = 0; i < neq; i++)
printf("%12.6f%s", y[i], (i+1)%3?" ":"\n");
itask = 0;
/* Obtain the solution at 5 equally spaced values of t.*/
for (j = 0; j < 5; j++)
{
tout = tout + dt;
/*
* nag_dae_ivp_dassl_gen (d02nec)
* dassl integrator
*/
nag_dae_ivp_dassl_gen(neq, &t, tout, y, ydot, rtol, atol, &itask, res,
jac, icom, com, lcom, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_dae_ivp_dassl_gen (d02nec).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("%8.4f", t);
for (i = 0; i < neq; i++)
printf("%12.6f%s", y[i], (i+1)%3?" ":"\n");
/*
* nag_dae_ivp_dassl_cont (d02mcc)
* dassl method continuation resetting function
*/
nag_dae_ivp_dassl_cont(icom);
}
printf("\n");
printf(" The integrator completed task, ITASK = %4ld\n", itask);
END:
NAG_FREE(atol);
NAG_FREE(com);
NAG_FREE(rtol);
NAG_FREE(y);
NAG_FREE(ydot);
NAG_FREE(comm.iuser);
NAG_FREE(icom);
return exit_status;
}
static void NAG_CALL res(Integer neq, double t, const double y[],
const double ydot[], double r[], Integer *ires,
Nag_Comm *comm)
{
if (comm->user[0] == -1.0)
{
printf("(User-supplied callback res, first invocation.)\n");
comm->user[0] = 0.0;
}
r[0] = (-(0.040e0*y[0])) + 1.00e4*y[1]*y[2] - ydot[0];
r[1] = 0.040e0*y[0] - 1.00e4*y[1]*y[2]-3.00e7*y[1]*y[1] - ydot[1];
r[2] = 3.00e7*y[1]*y[1] - ydot[2];
return;
}
static void NAG_CALL jac(Integer neq, double t, const double y[],
const double ydot[], double *pd, double cj,
Nag_Comm *comm)
{
Integer ml, mu;
if (comm->user[1] == -1.0)
{
printf("(User-supplied callback jac, first invocation.)\n");
comm->user[1] = 0.0;
}
ml = comm->iuser[0];
mu = comm->iuser[1];
myjac(neq, ml, mu, t, y, ydot, pd, cj);
return;
}
static void NAG_CALL myjac(Integer neq, Integer ml, Integer mu, double t,
const double y[], const double ydot[], double *pd,
double cj)
{
Integer md, ms;
Integer pdpd;
pdpd = 2*ml+mu+1;
#define PD(I, J) pd[(J-1)*pdpd + I - 1]
/* Main diagonal PDFULL(i,i), i=1,neq*/
md = mu + ml + 1;
PD(md, 1) = -0.040e0 - cj;
PD(md, 2) = -1.00e4*y[2] - 6.00e7*y[1] - cj;
PD(md, 3) = -cj;
/* 1 Sub-diagonal PDFULL(i+1:i), i=1,neq-1*/
ms = md + ml;
PD(ms, 1) = 0.040e0;
PD(ms, 2) = 6.00e7*y[1];
/* First super-diagonal PDFULL(i-1,i), i=2, neq*/
ms = md - 1;
PD(ms, 2) = 1.00e4*y[2];
PD(ms, 3) = -1.00e4*y[1];
/* Second super-diagonal PDFULL(i-2,i), i=3, neq*/
ms = md-2;
PD(ms, 3) = 1.00e4*y[1];
return;
}