/* nag_inteq_volterra2 (d05bac) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static double NAG_CALL sol(double t); static double NAG_CALL cf(double t, Nag_Comm *comm); static double NAG_CALL ck(double t, Nag_Comm *comm); static double NAG_CALL cg(double s, double y, Nag_Comm *comm); #ifdef __cplusplus } #endif int main(void) { /* Scalars */ double alim = 0.0, tlim = 20.0, tol = 1.e-4; double h, hi, si, thresh; Integer exit_status = 0; Integer iorder = 6, nmesh = 6; Integer i, lwk; /* Arrays */ double *errest = 0, *work = 0, *yn = 0; /* NAG types */ Nag_Comm comm; NagError fail; Nag_ODEMethod method = Nag_Adams; INIT_FAIL(fail); printf("nag_inteq_volterra2 (d05bac) Example Program Results\n"); lwk = 10 * nmesh + 6; if ( !(work = NAG_ALLOC(lwk, double)) || !(yn = NAG_ALLOC(nmesh, double)) || !(errest = NAG_ALLOC(nmesh, double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } h = (tlim - alim)/(double) (nmesh); thresh = nag_machine_precision; /* nag_inteq_volterra2 (d05bac). Nonlinear Volterra convolution equation, second kind. */ nag_inteq_volterra2(ck, cg, cf, method, iorder, alim, tlim, tol, nmesh, thresh, work, lwk, yn, errest, &comm, &fail); /* Loop until the supplied workspace is big enough. */ while (fail.code == NW_OUT_OF_WORKSPACE) { lwk = work[0]; if (work) NAG_FREE(work); if (!(work = NAG_ALLOC(lwk, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } nag_inteq_volterra2(ck, cg, cf, method, iorder, alim, tlim, tol, nmesh, thresh, work, lwk, yn, errest, &comm, &fail); } if (fail.code != NE_NOERROR) { printf("Error from nag_inteq_volterra2 (d05bac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nSize of workspace = %ld\n", lwk); printf("Tolerance = %f\n\n", tol); printf(" t Approx. Sol. True Sol. Est. Error Actual Error\n"); hi = 0.0; for (i = 0; i < nmesh; i++) { hi += h; si = sol(hi); printf("%7.2f%14.5f%14.5f%15.5e%15.5e\n", alim + hi, yn[i], si, errest[i], fabs((yn[i] - si)/si)); } END: if (errest) NAG_FREE(errest); if (yn) NAG_FREE(yn); if (work) NAG_FREE(work); return exit_status; } static double NAG_CALL sol(double t) { return log(t + exp(1.0)); } static double NAG_CALL cf(double t, Nag_Comm *comm) { return exp(-t); } static double NAG_CALL ck(double t, Nag_Comm *comm) { return exp(-t); } static double NAG_CALL cg(double s, double y, Nag_Comm *comm) { return y + exp(-y); }