/* nag_ode_bvp_fd_lin_gen (d02gbc) 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 fcnf(Integer neq, double x, double f[], Nag_User *comm); #ifdef __cplusplus } #endif #define NEQ 2 #define MNP 70 #define C(I, J) c[(I) *tdc + J] #define D(I, J) d[(I) *tdd + J] #define Y(I, J) y[(I) *tdy + J] int main(void) { Integer exit_status = 0, i, j, mnp, neq, np, tdc, tdd, tdy; NagError fail; Nag_User comm; double a, b, *c = 0, *d = 0, eps, *gam = 0, tol, *x = 0, *y = 0; INIT_FAIL(fail); printf("nag_ode_bvp_fd_lin_gen (d02gbc) Example Program Results\n"); /* For communication with function fcnf() * assign address of eps to comm.p. */ comm.p = (Pointer)&eps; neq = NEQ; mnp = MNP; tol = 1.0e-3; np = 0; a = 0.0; b = 1.0; if (mnp >= 32 && neq >= 2) { if (!(c = NAG_ALLOC(NEQ*NEQ, double)) || !(d = NAG_ALLOC(NEQ*NEQ, double)) || !(gam = NAG_ALLOC(NEQ, double)) || !(x = NAG_ALLOC(MNP, double)) || !(y = NAG_ALLOC(NEQ*MNP, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tdc = neq; tdd = neq; tdy = mnp; } else { exit_status = 1; return exit_status; } for (i = 0; i < neq; ++i) { gam[i] = 0.0; for (j = 0; j < neq; ++j) { C(i, j) = 0.0; D(i, j) = 0.0; } } C(0, 0) = 1.0; D(1, 0) = 1.0; gam[1] = 1.0; for (i = 1; i <= 2; ++i) { eps = pow(10.0, (double) -i); printf("\nProblem with epsilon = %7.4f\n", eps); /* nag_ode_bvp_fd_lin_gen (d02gbc). * Ordinary differential equations solver, for general * linear two-point boundary value problems, using a finite * difference technique with deferred correction */ nag_ode_bvp_fd_lin_gen(neq, fcnf, NULLFN, a, b, c, d, gam, mnp, &np, x, y, tol, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_bvp_fd_lin_gen (d02gbc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nApproximate solution on final mesh of %ld points\n", np); printf(" X(I) Y(1,I)\n"); for (j = 0; j < np; ++j) printf("%9.4f %9.4f\n", x[j], Y(0, j)); } END: if (c) NAG_FREE(c); if (d) NAG_FREE(d); if (gam) NAG_FREE(gam); if (x) NAG_FREE(x); if (y) NAG_FREE(y); return exit_status; } static void NAG_CALL fcnf(Integer neq, double x, double f[], Nag_User *comm) { #define F(I, J) f[(I) *neq+J] double *eps = (double *) comm->p; F(0, 0) = 0.0; F(0, 1) = 1.0; F(1, 0) = 0.0; F(1, 1) = -1.0/ *eps; }