/* 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 #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL fcn(Integer neq, double x, double eps, double y[], double f[], Nag_User *comm); static void NAG_CALL g(Integer neq, double eps, double ya[], double yb[], double bc[], Nag_User *comm); static void NAG_CALL jaceps(Integer neq, double x, double eps, double y[], double f[], Nag_User *comm); static void NAG_CALL jacgep(Integer neq, double eps, double ya[], double yb[], double bcep[], Nag_User *comm); static void NAG_CALL jacobf(Integer neq, double x, double eps, double y[], double f[], Nag_User *comm); static void NAG_CALL jacobg(Integer neq, double eps, double ya[], 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(int argc, char *argv[]) { FILE *fpout; 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); /* Check for command-line IO options */ fpout = nag_example_file_io(argc, argv, "-results", NULL); fprintf(fpout, "nag_ode_bvp_fd_nonlin_gen (d02rac) Example Program Results\n"); fprintf(fpout, "\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))) { fprintf(fpout, "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) { fprintf(fpout, "Error from nag_ode_bvp_fd_nonlin_gen (d02rac).\n%s\n", fail.message); exit_status = 1; goto END; } fprintf(fpout, "Solution on final mesh of %ld points \n", np); fprintf(fpout, " X Y(1) Y(2) Y(3)\n"); for (j = 0; j < np; ++j) { fprintf(fpout, " %9.3f ", x[j]); for (i = 0; i < neq; ++i) fprintf(fpout, " %9.4f", Y(i, j)); fprintf(fpout, "\n"); } fprintf(fpout, "\n\nMaximum estimated error by components \n"); for (i = 1; i <= 3; ++i) fprintf(fpout, " %11.2e", abt[i-1]); fprintf(fpout, " \n"); END: if (fpout != stdout) fclose(fpout); if (abt) NAG_FREE(abt); if (x) NAG_FREE(x); if (y) NAG_FREE(y); return exit_status; } #undef Y static void NAG_CALL fcn(Integer neq, double x, double eps, double y[], double f[], Nag_User *comm) { 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, double ya[], double yb[], double bc[], Nag_User *comm) { 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, double y[], double f[], Nag_User *comm) { 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, double ya[], double yb[], double bcep[], Nag_User *comm) { Integer i; for (i = 0; i < neq; ++i) bcep[i] = 0.0; } static void NAG_CALL jacobf(Integer neq, double x, double eps, double y[], double f[], Nag_User *comm) { Integer i, j; #define Y(I) y[(I) -1] #define F(I, J) f[((I) -1)*neq+(J) -1] 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, double ya[], double yb[], double aj[], double bj[], Nag_User *comm) { Integer i, j; #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] 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; }