/* nag_ode_ivp_adams_gen (d02cjc) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 3 revised, 1994. * Mark 7 revised, 2001. * */ #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL out(Integer neq, double *xsol, const double y[], Nag_User *comm); static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[], Nag_User *comm); static double NAG_CALL g(Integer neq, double x, const double y[], Nag_User *comm); #ifdef __cplusplus } #endif struct user { double xend, h; Integer k; }; int main(void) { Integer exit_status = 0, i, j, neq; Nag_User comm; double pi, tol, x, y[3]; struct user s; NagError fail; INIT_FAIL(fail); printf("nag_ode_ivp_adams_gen (d02cjc) Example Program Results\n"); /* For communication with function out() * assign address of user defined structure * to Nag pointer. */ comm.p = (Pointer)&s; neq = 3; s.xend = 10.0; /* nag_pi (x01aac). * pi */ pi = nag_pi; printf("\nCase 1: intermediate output, root-finding\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); printf("\n Calculation with tol = %10.1e\n", tol); x = 0.0; y[0] = 0.5; y[1] = 0.5; y[2] = pi / 5.0; s.k = 4; s.h = (s.xend - x) / (double)(s.k + 1); printf("\n X Y(1) Y(2) Y(3)\n"); /* nag_ode_ivp_adams_gen (d02cjc). * Ordinary differential equation solver using a * variable-order variable-step Adams method (Black Box) */ nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out, g, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n Root of Y(1) = 0.0 at %7.3f\n", x); printf("\n Solution is"); for (i = 0; i < 3; ++i) printf("%10.5f", y[i]); printf("\n"); } printf("\n\nCase 2: no intermediate output, root-finding\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); printf("\n Calculation with tol = %10.1e\n", tol); x = 0.0; y[0] = 0.5; y[1] = 0.5; y[2] = pi / 5.0; /* nag_ode_ivp_adams_gen (d02cjc), see above. */ nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN, g, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n Root of Y(1) = 0.0 at %7.3f\n", x); printf("\n Solution is"); for (i = 0; i < 3; ++i) printf("%10.5f", y[i]); printf("\n"); } printf("\n\nCase 3: intermediate output, no root-finding\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); printf("\n Calculation with tol = %10.1e\n", tol); x = 0.0; y[0] = 0.5; y[1] = 0.5; y[2] = pi / 5.0; s.k = 4; s.h = (s.xend - x) / (double)(s.k + 1); printf("\n X Y(1) Y(2) Y(3)\n"); /* nag_ode_ivp_adams_gen (d02cjc), see above. */ nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out, NULLDFN, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n", fail.message); exit_status = 1; goto END; } } printf("\n\nCase 4: no intermediate output, no root-finding"); printf(" ( integrate to xend)\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); printf("\n Calculation with tol = %10.1e\n", tol); x = 0.0; y[0] = 0.5; y[1] = 0.5; y[2] = pi / 5.0; printf("\n X Y(1) Y(2) Y(3)\n"); printf("%8.2f", x); for (i = 0; i < 3; ++i) printf("%13.5f", y[i]); printf("\n"); /* nag_ode_ivp_adams_gen (d02cjc), see above. */ nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN, NULLDFN, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("%8.2f", x); for (i = 0; i < 3; ++i) printf("%13.5f", y[i]); printf("\n"); } END: return exit_status; } static void NAG_CALL out(Integer neq, double *xsol, const double y[], Nag_User *comm) { Integer i; struct user *s = (struct user *) comm->p; printf("%8.2f", *xsol); for (i = 0; i < 3; ++i) { printf("%13.5f", y[i]); } printf("\n"); *xsol = s->xend - (double) s->k * s->h; s->k--; } /* out */ static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[], Nag_User *comm) { double pwr; f[0] = tan(y[2]); f[1] = -0.032*tan(y[2])/y[1] - 0.02*y[1]/cos(y[2]); pwr = y[1]; f[2] = -0.032/(pwr*pwr); } /* fcn */ static double NAG_CALL g(Integer neq, double x, const double y[], Nag_User *comm) { return y[0]; } /* g */