/* nag_ode_ivp_bdf_gen (d02ejc) Example Program. * * Copyright 1992 Numerical Algorithms Group. * * Mark 3, 1992. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[], Nag_User *comm); static void NAG_CALL pederv(Integer neq, double x, const double y[], double pw[], Nag_User *comm); static double NAG_CALL g(Integer neq, double x, const double y[], Nag_User *comm); static void NAG_CALL out(Integer neq, double *tsol, const double y[], Nag_User *comm); #ifdef __cplusplus } #endif struct user { double xend, h; Integer k; }; #define NEQ 3 int main(void) { Integer exit_status = 0, i, j, neq; NagError fail; Nag_User comm; double tol, x, *y = 0; struct user s; INIT_FAIL(fail); printf("nag_ode_ivp_bdf_gen (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; if (neq >= 1) { if (!(y = NAG_ALLOC(neq, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } } else { exit_status = 1; return exit_status; } s.xend = 10.0; printf("\nCase 1: calculating Jacobian internally\n"); printf(" intermediate output, root-finding\n\n"); for (j = 3; j <= 4; ++j) { tol = pow(10.0, -(double) j); printf("\n Calculation with tol = %10.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); printf(" X Y(1) Y(2) Y(3)\n"); /* nag_ode_ivp_bdf_gen (d02ejc). * Ordinary differential equations solver, stiff, initial * value problems using the Backward Differentiation * Formulae */ nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, out, g, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" Root of Y(1)-0.9 at %5.3f\n", x); printf(" Solution is "); for (i = 0; i < 3; ++i) printf("%7.5f ", y[i]); printf("\n"); } printf("\nCase 2: calculating Jacobian by pederv\n"); printf(" intermediate output, root-finding\n\n"); for (j = 3; j <= 4; ++j) { tol = pow(10.0, -(double) j); printf("\n Calculation with tol = %10.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); printf(" X Y(1) Y(2) Y(3)\n"); /* nag_ode_ivp_bdf_gen (d02ejc), see above. */ nag_ode_ivp_bdf_gen(neq, fcn, pederv, &x, y, s.xend, tol, Nag_Relative, out, g, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" Root of Y(1)-0.9 at %5.3f\n", x); printf(" Solution is "); for (i = 0; i < 3; ++i) printf("%7.5f ", y[i]); printf("\n"); } printf("\nCase 3: calculating Jacobian internally\n"); printf(" no intermediate output, root-finding\n\n"); for (j = 3; j <= 4; ++j) { tol = pow(10.0, -(double) j); printf("\n Calculation with tol = %10.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; /* nag_ode_ivp_bdf_gen (d02ejc), see above. */ nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, NULLFN, g, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" Root of Y(1)-0.9 at %5.3f\n", x); printf(" Solution is "); for (i = 0; i < 3; ++i) printf("%7.5f ", y[i]); printf("\n"); } printf("\nCase 4: calculating Jacobian internally\n"); printf(" intermediate output, no root-finding\n\n"); for (j = 3; j <= 4; ++j) { tol = pow(10.0, -(double) j); printf("\n Calculation with tol = %10.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); printf(" X Y(1) Y(2) Y(3)\n"); /* nag_ode_ivp_bdf_gen (d02ejc), see above. */ nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, out, NULLDFN, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\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"); } printf("\nCase 5: calculating Jacobian internally\n"); printf( " no intermediate output, no root-finding (integrate to xend)\n\n"); for (j = 3; j <= 4; ++j) { tol = pow(10.0, -(double) j); printf("\n Calculation with tol = %10.1e\n", tol); x = 0.0; y[0] = 1.0; y[1] = 0.0; y[2] = 0.0; printf(" 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_bdf_gen (d02ejc), see above. */ nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative, NULLFN, NULLDFN, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\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: if (y) NAG_FREE(y); return exit_status; } static void NAG_CALL fcn(Integer neq, double x, const 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 NAG_CALL pederv(Integer neq, double x, const 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 NAG_CALL g(Integer neq, double x, const double y[], Nag_User *comm) { return y[0]-0.9; } static void NAG_CALL out(Integer neq, double *xsol, const double y[], Nag_User *comm) { Integer j; struct user *s = (struct user *) comm->p; printf("%8.2f", *xsol); for (j = 0; j < 3; ++j) printf("%13.5f", y[j]); printf("\n"); *xsol = s->xend - (double) s->k * s->h; s->k--; }