/* nag_ode_dae_dassl_setup (d02mwc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*
*/
/* Pre-processor includes */
#include <math.h>
#include <nag.h>
#include <stdio.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);
#ifdef __cplusplus
}
#endif
int main(void) {
/* Scalars */
Integer exit_status = 0;
Integer i, itask, neq, maxord, licom, lcom;
double h0, hmax, g1, g2, t, tout;
/* Arrays */
static double ruser[2] = {-1.0, -1.0};
Integer *icom = 0;
double *atol = 0, *com = 0, *rtol = 0, *y = 0, *ydot = 0;
/* NAG types */
Nag_Boolean vector_tol;
Nag_Comm comm;
NagError fail;
INIT_FAIL(fail);
printf("nag_ode_dae_dassl_setup (d02mwc) Example Program Results\n\n");
/* For communication with user-supplied functions: */
comm.user = ruser;
/* Set problem size and allocate accordingly */
neq = 5;
maxord = 5;
licom = 50 + neq;
lcom = 40 + (maxord + 4 + neq) * neq;
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)) || !(icom = NAG_ALLOC(licom, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Use vector of tolerances */
vector_tol = Nag_TRUE;
for (i = 0; i < neq; i++) {
rtol[i] = 1.00e-8;
atol[i] = 1.00e-8;
}
/* Set up integrator to use supplied Jacobian, default step-sizes and
* vector tolerances using nag_ode_dae_dassl_setup (d02mwc).
*/
h0 = 0.0;
hmax = 0.0;
nag_ode_dae_dassl_setup(neq, maxord, Nag_AnalyticalJacobian, hmax, h0,
vector_tol, icom, licom, com, lcom, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_dae_dassl_setup (d02mwc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Set initial values */
y[0] = 1.0;
y[1] = 0.0;
y[2] = 0.0;
y[3] = 1.0;
y[4] = 1.0;
for (i = 0; i < neq; i++)
ydot[i] = 0.0;
t = 0.0;
tout = 1.0;
/* Print header and initial values */
printf("%7s%12s%12s%12s%12s%12s\n", "t", "y_1", "y_2", "y_3", "y_4", "y_5");
printf(" %6.4f", t);
for (i = 0; i < neq; i++)
printf("%11.6f%s", y[i], (i + 1) % 5 ? " " : "\n");
itask = 0;
while ((itask >= 0) && (itask <= 3) && (t < tout)) {
/* Integrate using nag_ode_dae_dassl_gen (d02nec). */
nag_ode_dae_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_ode_dae_dassl_gen (d02nec).\n%s\n", fail.message);
exit_status = 1;
goto END;
} else {
printf(" %6.4f", t);
for (i = 0; i < neq; i++)
printf("%11.6f%s", y[i], (i + 1) % 5 ? " " : "\n");
printf("\n d02nec returned with ITASK = %4" NAG_IFMT "\n\n", itask);
}
}
if ((itask >= 0) && (itask <= 3)) {
g1 = y[0] * y[0] + y[1] * y[1] - 1.0;
g2 = y[0] * y[2] + y[1] * y[3];
printf(" The position-level constraint G1 = %13.4e\n", g1);
printf(" The velocity-level constraint G2 = %13.4e\n", g2);
}
END:
NAG_FREE(atol);
NAG_FREE(com);
NAG_FREE(rtol);
NAG_FREE(y);
NAG_FREE(ydot);
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] = y[2] - ydot[0];
r[1] = y[3] - ydot[1];
r[2] = -y[4] * y[0] - ydot[2];
r[3] = -y[4] * y[1] - ydot[3] - 1.0;
r[4] = y[2] * y[2] + y[3] * y[3] - y[4] - y[1];
return;
}
static void NAG_CALL jac(Integer neq, double t, const double y[],
const double ydot[], double *pd, double cj,
Nag_Comm *comm) {
Integer pdpd;
if (comm->user[1] == -1.0) {
printf("(User-supplied callback jac, first invocation.)\n");
comm->user[1] = 0.0;
}
pdpd = neq;
#define PD(I, J) pd[(J - 1) * pdpd + I - 1]
PD(1, 1) = -cj;
PD(1, 3) = 1.0;
PD(2, 2) = -cj;
PD(2, 4) = 1.0;
PD(3, 1) = -y[4];
PD(3, 3) = -cj;
PD(3, 5) = -y[0];
PD(4, 2) = -y[4];
PD(4, 4) = -cj;
PD(4, 5) = -y[1];
PD(5, 2) = -1.0;
PD(5, 3) = 2.0 * y[2];
PD(5, 4) = 2.0 * y[3];
PD(5, 5) = -1.0;
return;
}