/* nag_dae_ivp_dassl_setup (d02mwc) Example Program. * * Copyright 2009, Numerical Algorithms Group. * * Mark 9, 2009. * */ /* Pre-processor includes */ #include #include #include #include #include #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 */ 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_dae_ivp_dassl_setup (d02mwc) Example Program Results\n\n"); /* 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_dae_ivp_dassl_setup (d02mwc). */ h0 = 0.0; hmax = 0.0; 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; } /* 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_dae_ivp_dassl_gen (d02nec). */ 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; } 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 = %4ld\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: if (atol) NAG_FREE(atol); if (com) NAG_FREE(com); if (rtol) NAG_FREE(rtol); if (y) NAG_FREE(y); if (ydot) NAG_FREE(ydot); if (icom) 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) { 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; 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; }