/* nag_dae_ivp_dassl_gen (d02nec) Example Program. * * Copyright 2009, Numerical Algorithms Group. * * Mark 9, 2009. * */ /* Pre-processor includes */ #include #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); 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; INIT_FAIL(fail); printf("nag_dae_ivp_dassl_gen (d02nec) Example Program Results\n\n"); /* 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: 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 (comm.iuser) NAG_FREE(comm.iuser); 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] = (-(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; 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; }