/* nag_ode_bvp_fd_nonlin_fixedbc (d02gac) 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, const double y[], double f[], Nag_User *comm); #ifdef __cplusplus } #endif #define NEQ 3 #define MNP 40 #define U(I, J) u[(I) *tdu + J] #define Y(I, J) y[(I) *tdy + J] #define V(I, J) v[(I) *tdv + J] int main(void) { Integer exit_status = 0, i, j, k, mnp, neq, np, tdu, tdv, tdy, *v = 0; NagError fail; Nag_User comm; double a, b, beta, tol, *u = 0, *x = 0, *y = 0; INIT_FAIL(fail); printf( "nag_ode_bvp_fd_nonlin_fixedbc (d02gac) Example Program Results\n"); /* For communication with function fcn() * assign address of beta to comm.p. */ comm.p = (Pointer)β neq = NEQ; mnp = MNP; tol = 0.001; np = 26; a = 0.0; b = 10.0; beta = 0.0; if (mnp >= 32 && neq >= 2) { if (!(u = NAG_ALLOC(neq*2, double)) || !(x = NAG_ALLOC(mnp, double)) || !(y = NAG_ALLOC(neq*mnp, double)) || !(v = NAG_ALLOC(neq*2, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tdu = 2; tdy = mnp; tdv = 2; } else { exit_status = 1; return exit_status; } for (i = 0; i < neq; ++i) for (j = 0; j < 2; ++j) { U(i, j) = 0.0; V(i, j) = 0; } V(2, 0) = 1; V(0, 1) = 1; V(2, 1) = 1; U(1, 1) = 1.0; U(0, 1) = b; x[0] = a; for (i = 2; i <= np-1; ++i) x[i-1] = ((double)(np-i)*a + (double)(i-1)*b)/ (double)(np-1); x[np-1] = b; for (k = 1; k <= 2; ++k) { printf("\nProblem with beta = %7.4f\n", beta); /* nag_ode_bvp_fd_nonlin_fixedbc (d02gac). * Ordinary differential equations solver, for simple * nonlinear two-point boundary value problems, using a * finite difference technique with deferred correction */ nag_ode_bvp_fd_nonlin_fixedbc(neq, fcn, a, b, u, v, mnp, &np, x, y, tol, &comm, &fail); if (fail.code == NE_NOERROR || fail.code == NE_CONV_ROUNDOFF) { if (fail.code == NE_CONV_ROUNDOFF) { printf("Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac)" ".\n%s\n", fail.message); exit_status = 2; } printf("\nSolution on final mesh of %ld points\n", np); printf(" X Y(1) Y(2) Y(3)\n"); for (i = 0; i <= np-1; ++i) { printf(" %9.4f ", x[i]); for (j = 0; j < neq; ++j) printf(" %9.4f ", Y(j, i)); printf("\n"); } beta += 0.2; } else { printf( "Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac).\n%s\n", fail.message); exit_status = 1; } } END: if (u) NAG_FREE(u); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (v) NAG_FREE(v); return exit_status; } static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[], Nag_User *comm) { double *beta = (double *) comm->p; f[0] = y[1]; f[1] = y[2]; f[2] = -y[0] * y[2] - *beta * (1.0-y[1]*y[1]); }