/* 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 out(Integer neq, double *xsol, double y[], Nag_User *comm); static void fcn(Integer neq, double x, double y[], double f[], Nag_User *comm); static double g(Integer neq, double x, double y[], Nag_User *comm); #ifdef __cplusplus } #endif struct user { double xend, h; Integer k; }; int main(void) { Integer i, j, neq; double x, pi, tol; double y[3]; Nag_User comm; struct user s; Vprintf("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; pi = X01AAC; Vprintf("\nCase 1: intermediate output, root-finding\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); Vprintf("\n Calculation with tol = %8.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); Vprintf("\n X Y(1) Y(2) Y(3)\n"); d02cjc(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out, g, &comm, NAGERR_DEFAULT); Vprintf("\n Root of Y(1) = 0.0 at %7.3f\n", x); Vprintf("\n Solution is"); for (i = 0; i < 3; ++i) Vprintf("%10.5f", y[i]); Vprintf("\n"); } Vprintf("\n\nCase 2: no intermediate output, root-finding\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); Vprintf("\n Calculation with tol = %8.1e\n", tol); x = 0.0; y[0] = 0.5; y[1] = 0.5; y[2] = pi / 5.0; d02cjc(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN, g, &comm, NAGERR_DEFAULT); Vprintf("\n Root of Y(1) = 0.0 at %7.3f\n", x); Vprintf("\n Solution is"); for (i = 0; i < 3; ++i) Vprintf("%10.5f", y[i]); Vprintf("\n"); } Vprintf("\n\nCase 3: intermediate output, no root-finding\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); Vprintf("\n Calculation with tol = %8.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); Vprintf("\n X Y(1) Y(2) Y(3)\n"); d02cjc(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out, NULLDFN, &comm, NAGERR_DEFAULT); } Vprintf("\n\nCase 4: no intermediate output, no root-finding"); Vprintf(" ( integrate to xend)\n"); for (j = 4; j <= 5; ++j) { tol = pow(10.0, (double)(-j)); Vprintf("\n Calculation with tol = %8.1e\n", tol); x = 0.0; y[0] = 0.5; y[1] = 0.5; y[2] = pi / 5.0; Vprintf("\n X Y(1) Y(2) Y(3)\n"); Vprintf("%8.2f", x); for (i = 0; i < 3; ++i) Vprintf("%13.5f", y[i]); Vprintf("\n"); d02cjc(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN, NULLDFN, &comm, NAGERR_DEFAULT); Vprintf("%8.2f", x); for (i = 0; i < 3; ++i) Vprintf("%13.5f", y[i]); Vprintf("\n"); } return EXIT_SUCCESS; } static void out(Integer neq, double *xsol, double y[], Nag_User *comm) { Integer i; struct user *s = (struct user *)comm->p; Vprintf("%8.2f", *xsol); for (i = 0; i < 3; ++i) { Vprintf("%13.5f", y[i]); } Vprintf("\n"); *xsol = s->xend - (double)s->k * s->h; s->k--; } /* out */ static void fcn(Integer neq, double x, 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 g(Integer neq, double x, double y[], Nag_User *comm) { return y[0]; } /* g */