/* nag_ode_ivp_rk_errass (d02pzc) 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 4 #define ZERO 0.0 #define ONE 1.0 #define THREE 3.0 #define ECC 0.7 int main(void) { Integer exit_status = 0, i, neq; NagError fail; Nag_ErrorAssess errass; Nag_ODE_RK opt; Nag_RK_method method; Nag_User comm; double errmax, hstart, pi, *rmserr = 0, tend, terrmx, tgot, *thres = 0, tol; double tstart, twant, *ygot = 0, *ymax = 0, *ypgot = 0, *ystart = 0; INIT_FAIL(fail); printf("nag_ode_ivp_rk_errass (d02pzc) 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)) || !(ypgot = NAG_ALLOC(neq, double)) || !(ystart = NAG_ALLOC(neq, double)) || !(ymax = NAG_ALLOC(NEQ, double)) || !(rmserr = 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] = ONE - ECC; ystart[1] = ZERO; ystart[2] = ZERO; ystart[3] = sqrt((ONE+ECC)/(ONE-ECC)); tend = THREE*pi; for (i = 0; i < neq; i++) thres[i] = 1.0e-10; errass = Nag_ErrorAssess_on; hstart = ZERO; tol = 1.0e-6; method = Nag_RK_7_8; /* 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 y3 y4\n\n"); printf("%8.3f %8.4f %8.4f %8.4f %8.4f\n", tstart, ystart[0], ystart[1], ystart[2], ystart[3]); twant = tend; do { /* 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); } while (fail.code == NE_RK_PDC_POINTS || fail.code == NE_STIFF_PROBLEM); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_rk_range (d02pcc).\n%s\n", fail.message); exit_status = 1; goto END; } else { printf("%8.3f %8.4f %8.4f %8.4f %8.4f\n\n", tgot, ygot[0], ygot[1], ygot[2], ygot[3]); /* nag_ode_ivp_rk_errass (d02pzc). * A function to provide global error assessment during an * integration with either nag_ode_ivp_rk_range (d02pcc) or * nag_ode_ivp_rk_onestep (d02pdc) */ nag_ode_ivp_rk_errass(neq, rmserr, &errmax, &terrmx, &opt, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_rk_errass (d02pzc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("Componentwise error assessment\n "); for (i = 0; i < neq; i++) printf("%11.2e ", rmserr[i]); printf("\n\n"); printf("Worst global error observed was %11.2e - " "it occurred at t = %6.3f\n\n", errmax, terrmx); printf("Cost of the integration in evaluations of f is %ld\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 (ypgot) NAG_FREE(ypgot); if (ystart) NAG_FREE(ystart); if (ymax) NAG_FREE(ymax); if (rmserr) NAG_FREE(rmserr); return exit_status; } static void NAG_CALL f(Integer neq, double t, const double y[], double yp[], Nag_User *comm) { double r, rp3; r = sqrt(y[0]*y[0]+y[1]*y[1]); rp3 = pow(r, 3.0); yp[0] = y[2]; yp[1] = y[3]; yp[2] = -y[0]/rp3; yp[3] = -y[1]/rp3; }