/* nag_ode_ivp_rk_range (d02pcc) Example Program. * * Copyright 1992 Numerical Algorithms Group. * * Mark 3, 1992. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL f(Integer neq, double t1, const double y[], double yp[], Nag_User *comm); #ifdef __cplusplus } #endif #define NEQ 2 #define ZERO 0.0 #define ONE 1.0 #define TWO 2.0 #define FOUR 4.0 int main(void) { Integer exit_status = 0, i, j, neq, nout; NagError fail; Nag_ErrorAssess errass; Nag_ODE_RK opt; Nag_RK_method method; Nag_User comm; double hstart, pi, tend, tgot, *thres = 0, tinc, tol, tstart, twant, *ygot = 0; double *ymax = 0, *ypgot = 0, *ystart = 0; INIT_FAIL(fail); printf("nag_ode_ivp_rk_range (d02pcc) Example Program Results\n"); /* Set initial conditions and input for nag_ode_ivp_rk_setup (d02pvc) */ neq = NEQ; if (neq >= 1) { if (!(thres = NAG_ALLOC(neq, double)) || !(ygot = NAG_ALLOC(neq, double)) || !(ymax = NAG_ALLOC(neq, double)) || !(ypgot = NAG_ALLOC(neq, double)) || !(ystart = NAG_ALLOC(neq, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { exit_status = 1; return exit_status; } /* nag_pi (x01aac). * pi */ pi = nag_pi; tstart = ZERO; ystart[0] = ZERO; ystart[1] = ONE; tend = TWO*pi; for (i = 0; i < neq; i++) thres[i] = 1.0e-8; errass = Nag_ErrorAssess_off; hstart = ZERO; method = Nag_RK_2_3; /* * Set control for output */ nout = 8; tinc = (tend-tstart)/nout; for (i = 1; i <= 2; i++) { if (i == 1) tol = 1.0e-3; else tol = 1.0e-4; /* nag_ode_ivp_rk_setup (d02pvc). * Setup function for use with nag_ode_ivp_rk_range (d02pcc) * and/or nag_ode_ivp_rk_onestep (d02pdc) */ nag_ode_ivp_rk_setup(neq, tstart, ystart, tend, tol, thres, method, Nag_RK_range, errass, hstart, &opt, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_rk_setup (d02pvc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nCalculation with tol = %10.1e\n\n", tol); printf(" t y1 y2\n\n"); printf("%8.3f %8.3f %8.3f\n", tstart, ystart[0], ystart[1]); for (j = nout-1; j >= 0; j--) { twant = tend - j*tinc; /* nag_ode_ivp_rk_range (d02pcc). * Ordinary differential equations solver, initial value * problems over a range using Runge-Kutta methods */ nag_ode_ivp_rk_range(neq, f, twant, &tgot, ygot, ypgot, ymax, &opt, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_rk_range (d02pcc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("%8.3f %8.3f %8.3f\n", tgot, ygot[0], ygot[1]); } printf("\nCost of the integration in evaluations of f is" " %ld\n\n", opt.totfcn); /* nag_ode_ivp_rk_free (d02ppc). * Freeing function for use with the Runge-Kutta suite (d02p * functions) */ nag_ode_ivp_rk_free(&opt); } END: if (thres) NAG_FREE(thres); if (ygot) NAG_FREE(ygot); if (ymax) NAG_FREE(ymax); if (ypgot) NAG_FREE(ypgot); if (ystart) NAG_FREE(ystart); return exit_status; } static void NAG_CALL f(Integer neq, double t, const double y[], double yp[], Nag_User *comm) { yp[0] = y[1]; yp[1] = -y[0]; }