/* nag_ode_ivp_adams_roots (d02qfc) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void ftry02(Integer neqf, double x, double y[], double yp[], Nag_User *comm); static double gtry02(Integer neqf, double x, double y[], double yp[], Integer k, Nag_User *comm); #ifdef __cplusplus } #endif #define NEQF 2 #define NEQG 2 int main(void) { Nag_Boolean alter_g, crit, one_step, sophist, vectol; Integer exit_status=0, i, max_step, neqf, neqg; NagError fail; Nag_ODE_Adams opt; Nag_Start state; double *atol=0, *rtol=0, t, tcrit, tout, *y=0; INIT_FAIL(fail); Vprintf("nag_ode_ivp_adams_roots (d02qfc) Example Program Results\n"); neqf = NEQF; neqg = NEQG; if (neqf<1) { exit_status = 1; return exit_status; } else { if ( !( y = NAG_ALLOC(neqf, double)) || !( atol = NAG_ALLOC(neqf, double)) || !( rtol = NAG_ALLOC(neqf, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } } tcrit = 10.0; state = Nag_NewStart; vectol = Nag_TRUE; one_step = Nag_FALSE; crit = Nag_TRUE; max_step = 0; sophist = Nag_TRUE; for (i = 0; i <= 1; ++i) { rtol[i] = 0.0001; atol[i] = 1e-06; } /* nag_ode_ivp_adams_setup (d02qwc). * Setup function for nag_ode_ivp_adams_roots (d02qfc) */ nag_ode_ivp_adams_setup(&state, neqf, vectol, atol, rtol, one_step, crit, tcrit, 0.0, max_step, neqg, &alter_g, sophist, &opt, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_ode_ivp_adams_setup (d02qwc).\n%s\n", fail.message); exit_status = 1; goto END; } t = 0.0; tout = tcrit; y[0] = 0.0; y[1] = 1.0; do { /* nag_ode_ivp_adams_roots (d02qfc). * Ordinary differential equation solver using Adams method * (sophisticated use) */ nag_ode_ivp_adams_roots(neqf, ftry02, &t, y, tout, gtry02, NAGUSER_DEFAULT, &opt, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_ode_ivp_adams_roots (d02qfc).\n%s\n", fail.message); exit_status = 1; goto END; } if (opt.root) { Vprintf("\nRoot at %11.5e\n", t); Vprintf("for event equation %1ld", opt.index); Vprintf(" with type %1ld", opt.type); Vprintf(" and residual %11.5e\n", opt.resids[opt.index-1]); Vprintf(" Y(1) = %11.5e Y'(1) = %11.5e\n", y[0], opt.yp[0]); for (i = 1; i <= neqg; ++i) { if (i != opt.index && opt.events[i-1] != 0) { Vprintf("and also for event equation %1ld", i); Vprintf(" with type %1ld", opt.events[i-1]); Vprintf(" and residual %11.5e\n", opt.resids[i-1]); } } } } while (opt.tcurr < tout && opt.root); /* Free the memory which was allocated by * nag_ode_ivp_adams_setup (d02qwc) to the pointers inside opt. */ /* nag_ode_ivp_adams_free (d02qyc). * Freeing function for use with nag_ode_ivp_adams_roots * (d02qfc) */ nag_ode_ivp_adams_free(&opt); END: if (y) NAG_FREE(y); if (atol) NAG_FREE(atol); if (rtol) NAG_FREE(rtol); return exit_status; } static void ftry02(Integer neqf, double x, double y[], double yp[], Nag_User *comm) { yp[0] = y[1]; yp[1] = -y[0]; } /* ftry02 */ static double gtry02(Integer neqf, double x, double y[], double yp[], Integer k, Nag_User *comm) { if (k == 1) return yp[0]; else return y[0]; } /* gtry02 */