/* nag_ode_ivp_rk_interp (d02pxc) 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 NWANT 1 #define ZERO 0.0 #define ONE 1.0 #define TWO 2.0 #define FOUR 4.0 int main(void) { static Integer use_comm[1] = {1}; Integer exit_status = 0, i, j, neq, nout, nwant; NagError fail; Nag_ErrorAssess errass; Nag_ODE_RK opt; Nag_RK_method method; Nag_User comm; double hstart, pi, tend, *thres = 0, tinc, tnow, tol, tstart, twant, *ynow = 0; double *ypnow = 0, *ypwant = 0, *ystart = 0, *ywant = 0; INIT_FAIL(fail); printf("nag_ode_ivp_rk_interp (d02pxc) Example Program Results\n"); /* For communication with user-supplied functions: */ comm.p = (Pointer)&use_comm; /* Set initial conditions and input for nag_ode_ivp_rk_setup (d02pvc) */ neq = NEQ; nwant = NWANT; if (neq >= 1) { if (!(thres = NAG_ALLOC(neq, double)) || !(ynow = NAG_ALLOC(neq, double)) || !(ypnow = NAG_ALLOC(neq, double)) || !(ystart = NAG_ALLOC(neq, double)) || !(ywant = NAG_ALLOC(nwant, double)) || !(ypwant = NAG_ALLOC(nwant, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { exit_status = 1; return exit_status; } method = Nag_RK_4_5; /* 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; /* * Set control for output */ nwant = NWANT; nout = 16; tinc = tend/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_onestep, 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.4f %8.4f\n", tstart, ystart[0], ystart[1]); j = nout - 1; twant = tend - j*tinc; do { /* nag_ode_ivp_rk_onestep (d02pdc). * Ordinary differential equations solver, initial value * problems, one time step using Runge-Kutta methods */ nag_ode_ivp_rk_onestep(neq, f, &tnow, ynow, ypnow, &opt, &comm, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_ode_ivp_rk_onestep (d02pdc).\n%s\n", fail.message); exit_status = 1; goto END; } while (twant <= tnow) { /* nag_ode_ivp_rk_interp (d02pxc). * Ordinary differential equations solver, computes the * solution by interpolation anywhere on an integration step * taken by nag_ode_ivp_rk_onestep (d02pdc) */ nag_ode_ivp_rk_interp(neq, twant, Nag_SolDer, nwant, ywant, ypwant, f, &opt, &comm, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_ode_ivp_rk_interp (d02pxc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("%8.3f %8.4f %8.4f\n", twant, ywant[0], ypwant[0]); j = j - 1; twant = tend - j*tinc; } } while (tnow < tend); 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: NAG_FREE(thres); NAG_FREE(ynow); NAG_FREE(ypnow); NAG_FREE(ystart); NAG_FREE(ywant); NAG_FREE(ypwant); return exit_status; } static void NAG_CALL f(Integer neq, double t, const double y[], double yp[], Nag_User *comm) { Integer *use_comm = (Integer *)comm->p; if (use_comm[0]) { printf("(User-supplied callback f, first invocation.)\n"); use_comm[0] = 0; } yp[0] = y[1]; yp[1] = -y[0]; }