/* nag_ode_bvp_fd_nonlin_gen (d02rac) Example Program. * * Copyright 1992 Numerical Algorithms Group. * * Mark 3, 1992. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL fcn(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm); static void NAG_CALL g(Integer neq, double eps, const double ya[], const double yb[], double bc[], Nag_User *comm); static void NAG_CALL jaceps(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm); static void NAG_CALL jacgep(Integer neq, double eps, const double ya[], const double yb[], double bcep[], Nag_User *comm); static void NAG_CALL jacobf(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm); static void NAG_CALL jacobg(Integer neq, double eps, const double ya[], const double yb[], double aj[], double bj[], Nag_User *comm); #ifdef __cplusplus } #endif #define NEQ 3 #define MNP 40 #define Y(I, J) y[(I) *tdy + J] int main(void) { static Integer use_comm[6] = {1, 1, 1, 1, 1, 1}; double *abt = 0; double deleps; double tol; double *x = 0, *y = 0; Integer exit_status = 0; Integer i, j; Integer np; Integer numbeg, nummix; Integer neq, mnp, tdy; Nag_User comm; NagError fail; INIT_FAIL(fail); printf( "nag_ode_bvp_fd_nonlin_gen (d02rac) Example Program Results\n"); /* For communication with user-supplied functions: */ comm.p = (Pointer)&use_comm; printf("\nCalculation using analytic Jacobians\n\n"); neq = NEQ; mnp = MNP; if (neq >= 1) { if (!(abt = NAG_ALLOC(neq, double)) || !(x = NAG_ALLOC(mnp, double)) || !(y = NAG_ALLOC(neq*mnp, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tdy = mnp; } else { exit_status = 1; return exit_status; } tol = 1.0e-4; np = 17; numbeg = 2; nummix = 0; x[0] = 0.0; x[np-1] = 10.0; deleps = 0.1; /* nag_ode_bvp_fd_nonlin_gen (d02rac). * Ordinary differential equations solver, for general * nonlinear two-point boundary value problems, using a * finite difference technique with deferred correction */ nag_ode_bvp_fd_nonlin_gen(neq, &deleps, fcn, numbeg, nummix, g, Nag_DefInitMesh, mnp, &np, x, y, tol, abt, jacobf, jacobg, jaceps, jacgep, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_bvp_fd_nonlin_gen (d02rac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("Solution on final mesh of %ld points \n", np); printf(" X Y(1) Y(2) Y(3)\n"); for (j = 0; j < np; ++j) { printf(" %9.3f ", x[j]); for (i = 0; i < neq; ++i) printf(" %9.4f", Y(i, j)); printf("\n"); } printf("\n\nMaximum estimated error by components \n"); for (i = 1; i <= 3; ++i) printf(" %11.2e", abt[i-1]); printf(" \n"); END: NAG_FREE(abt); NAG_FREE(x); NAG_FREE(y); return exit_status; } #undef Y static void NAG_CALL fcn(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm) { Integer *use_comm = (Integer *)comm->p; if (use_comm[0]) { printf("(User-supplied callback fcn, first invocation.)\n"); use_comm[0] = 0; } f[0] = y[1]; f[1] = y[2]; f[2] = -y[0] * y[2] - (1.0 - y[1]*y[1])*2.0*eps; } static void NAG_CALL g(Integer neq, double eps, const double ya[], const double yb[], double bc[], Nag_User *comm) { Integer *use_comm = (Integer *)comm->p; if (use_comm[1]) { printf("(User-supplied callback g, first invocation.)\n"); use_comm[1] = 0; } bc[0] = ya[0]; bc[1] = ya[1]; bc[2] = yb[1] - 1.0; } /* g */ static void NAG_CALL jaceps(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm) { Integer *use_comm = (Integer *)comm->p; if (use_comm[2]) { printf("(User-supplied callback jaceps, first invocation.)\n"); use_comm[2] = 0; } f[0] = 0.0; f[1] = 0.0; f[2] = (1.0 - y[1]*y[1]) * -2.0; } static void NAG_CALL jacgep(Integer neq, double eps, const double ya[], const double yb[], double bcep[], Nag_User *comm) { Integer i; Integer *use_comm = (Integer *)comm->p; if (use_comm[3]) { printf("(User-supplied callback jacgep, first invocation.)\n"); use_comm[3] = 0; } for (i = 0; i < neq; ++i) bcep[i] = 0.0; } static void NAG_CALL jacobf(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm) { Integer i, j; Integer *use_comm = (Integer *)comm->p; #define Y(I) y[(I) -1] #define F(I, J) f[((I) -1)*neq+(J) -1] if (use_comm[4]) { printf("(User-supplied callback jacobf, first invocation.)\n"); use_comm[4] = 0; } for (i = 1; i <= neq; ++i) { for (j = 1; j <= neq; ++j) F(i, j) = 0.0; } F(1, 2) = 1.0; F(2, 3) = 1.0; F(3, 1) = -Y(3); F(3, 2) = Y(2) * 4.0 * eps; F(3, 3) = -Y(1); } static void NAG_CALL jacobg(Integer neq, double eps, const double ya[], const double yb[], double aj[], double bj[], Nag_User *comm) { Integer i, j; Integer *use_comm = (Integer *)comm->p; #define YA(I) ya[(I) -1] #define YB(I) yb[(I) -1] #define AJ(I, J) aj[((I) -1)*neq+(J) -1] #define BJ(I, J) bj[((I) -1)*neq+(J) -1] if (use_comm[5]) { printf("(User-supplied callback jacobg, first invocation.)\n"); use_comm[5] = 0; } for (i = 1; i <= neq; ++i) for (j = 1; j <= neq; ++j) { AJ(i, j) = 0.0; BJ(i, j) = 0.0; } AJ(1, 1) = 1.0; AJ(2, 2) = 1.0; BJ(3, 2) = 1.0; }