/* nag_ode_ivp_adams_roots(d02qfc) Example Program * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 7 revised, 2001. * */ #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) { double y[NEQF], atol[NEQF], rtol[NEQF]; Boolean crit, alter_g, vectol, one_step, sophist; double t, tout, tcrit; Integer i, max_step, neqf, neqg; Nag_Start state; Nag_ODE_Adams opt; Vprintf("d02qfc Example Program Results\n"); neqf = NEQF; neqg = NEQG; tcrit = 10.0; state = Nag_NewStart; vectol = TRUE; one_step = FALSE; crit = TRUE; max_step = 0; sophist = TRUE; for (i = 0; i <= 1; ++i) { rtol[i] = 0.0001; atol[i] = 1e-06; } d02qwc(&state, neqf, vectol, atol, rtol, one_step, crit, tcrit, 0.0, max_step, neqg, &alter_g, sophist, &opt, NAGERR_DEFAULT); t = 0.0; tout = tcrit; y[0] = 0.0; y[1] = 1.0; do { d02qfc(neqf, ftry02, &t, y, tout, gtry02, NAGUSER_DEFAULT, &opt, NAGERR_DEFAULT); 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 * d02qwc to the pointers inside opt. */ d02qyc(&opt); return EXIT_SUCCESS; } 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 */