/* nag_ode_ivp_bdf_gen(d02ejc) Example Program * * Copyright 1992 Numerical Algorithms Group. * * Mark 3, 1992. * Mark 7 revised, 2001. * */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void fcn(Integer neq, double x, double y[], double f[], Nag_User *comm); static void pederv(Integer neq, double x, double y[], double pw[], Nag_User *comm); static double g(Integer neq, double x, double y[], Nag_User *comm); static void out(Integer neq, double *tsol, double y[], Nag_User *comm); #ifdef __cplusplus } #endif struct user { double xend, h; Integer k; }; #define NEQ 3 int main(void) { Integer neq; Integer i, j; double x, y[3]; double tol; Nag_User comm; struct user s; Vprintf("d02ejc Example Program Results\n"); /* For communication with function out() * assign address of user defined structure * to comm.p. */ comm.p = (Pointer)&s; neq = NEQ; s.xend = 10.0; Vprintf("\nCase 1: calculating Jacobian internally\n"); Vprintf(" intermediate output, root-finding\n\n"); for (j=3; j<=4; ++j) { tol = pow(10.0, -(double)j); Vprintf("\n Calculation with tol = %3.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; s.k = 4; s.h = (s.xend-x) /(double)(s.k+1); Vprintf( " X Y(1) Y(2) Y(3)\n"); d02ejc(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, out, g, &comm, NAGERR_DEFAULT); Vprintf(" Root of Y(1)-0.9 at %5.3f\n", x); Vprintf(" Solution is "); for (i=0; i<3; ++i) Vprintf("%7.5f ", y[i]); Vprintf("\n"); } Vprintf("\nCase 2: calculating Jacobian by pederv\n"); Vprintf(" intermediate output, root-finding\n\n"); for (j=3; j<=4; ++j) { tol = pow(10.0, -(double)j); Vprintf("\n Calculation with tol = %3.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; s.k = 4; s.h = (s.xend-x) /(double)(s.k+1); Vprintf( " X Y(1) Y(2) Y(3)\n"); d02ejc(neq, fcn, pederv, &x, y, s.xend, tol, Nag_Relative, out, g, &comm, NAGERR_DEFAULT); Vprintf(" Root of Y(1)-0.9 at %5.3f\n", x); Vprintf(" Solution is "); for (i=0; i<3; ++i) Vprintf("%7.5f ", y[i]); Vprintf("\n"); } Vprintf("\nCase 3: calculating Jacobian internally\n"); Vprintf(" no intermediate output, root-finding\n\n"); for (j=3; j<=4; ++j) { tol = pow(10.0, -(double)j); Vprintf("\n Calculation with tol = %3.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; d02ejc(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, NULLFN, g, &comm, NAGERR_DEFAULT); Vprintf(" Root of Y(1)-0.9 at %5.3f\n", x); Vprintf(" Solution is "); for (i=0; i<3; ++i) Vprintf("%7.5f ", y[i]); Vprintf("\n"); } Vprintf("\nCase 4: calculating Jacobian internally\n"); Vprintf(" intermediate output, no root-finding\n\n"); for (j=3; j<=4; ++j) { tol = pow(10.0, -(double)j); Vprintf("\n Calculation with tol = %3.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; s.k = 4; s.h = (s.xend-x) /(double)(s.k+1); Vprintf( " X Y(1) Y(2) Y(3)\n"); d02ejc(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, out, NULLDFN, &comm, NAGERR_DEFAULT); Vprintf("%8.2f", x); for (i=0; i<3; ++i) Vprintf("%13.5f", y[i]); Vprintf("\n"); } Vprintf("\nCase 5: calculating Jacobian internally\n"); Vprintf(" no intermediate output, no root-finding (integrate to xend)\n\n"); for (j=3; j<=4; ++j) { tol = pow(10.0, -(double)j); Vprintf("\n Calculation with tol = %3.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; Vprintf( " 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"); d02ejc(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, 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 fcn(Integer neq, double x, double y[], double f[], Nag_User *comm) { f[0] = y[0] * -0.04 + y[1] * 1e4 * y[2]; f[1] = y[0] * 0.04 - y[1] * 1e4 * y[2] - y[1] * 3e7 * y[1]; f[2] = y[1] * 3e7 * y[1]; } static void pederv(Integer neq, double x, double y[], double pw[], Nag_User *comm) { #define PW(I,J) pw[((I)-1)*neq + (J)-1] PW(1,1) = -0.04; PW(1,2) = y[2] * 1e4; PW(1,3) = y[1] * 1e4; PW(2,1) = 0.04; PW(2,2) = y[2] * -1e4 - y[1] * 6e7; PW(2,3) = y[1] * -1e4; PW(3,1) = 0.0; PW(3,2) = y[1] * 6e7; PW(3,3) = 0.0; } static double g(Integer neq, double x, double y[], Nag_User *comm) { return y[0]-0.9; } static void out(Integer neq, double *xsol, double y[], Nag_User *comm) { Integer j; struct user *s = (struct user *)comm->p; Vprintf("%8.2f", *xsol); for (j=0; j<3; ++j) Vprintf("%13.5f", y[j]); Vprintf("\n"); *xsol = s->xend - (double)s->k * s->h; s->k--; }