/* nag_inteq_abel1_weak (d05bec) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static double sol1(double t); static double sol2(double t); static double NAG_CALL ck1(double t, Nag_Comm *comm); static double NAG_CALL cf1(double t, Nag_Comm *comm); static double NAG_CALL cg1(double s, double y, Nag_Comm *comm); static double NAG_CALL ck2(double t, Nag_Comm *comm); static double NAG_CALL cf2(double t, Nag_Comm *comm); static double NAG_CALL cg2(double s, double y, Nag_Comm *comm); #ifdef __cplusplus } #endif int main(void) { /* Scalars */ double err, errmax, h, hi1, soln, t, tlim, tolnl; Integer exit_status = 0; Integer iorder = 4; Integer i, iskip, exno, nmesh, lrwsav; /* Arrays */ double *rwsav = 0, *yn = 0; /* NAG types */ Nag_Comm comm; NagError fail; Nag_WeightMode wtmode; INIT_FAIL(fail); printf("nag_inteq_abel1_weak (d05bec) Example Program Results\n"); nmesh = pow(2, 6) + 7; lrwsav = (2 * iorder + 6) * nmesh + 8 * pow(iorder, 2) - 16 * iorder + 1; if ( !(yn = NAG_ALLOC(nmesh, double)) || !(rwsav = NAG_ALLOC(lrwsav, double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } tolnl = sqrt(nag_machine_precision); for (exno = 1; exno <= 2; exno++) { printf("\nExample %ld\n\n", exno); if (exno==1) { tlim = 7.0; iskip = 5; h = tlim/(double) (nmesh - 1); wtmode = Nag_InitWeights; yn[0] = 0.0; /* nag_inteq_abel1_weak (d05bec). Nonlinear convolution Volterra-Abel equation, first kind, weakly singular. */ nag_inteq_abel1_weak(ck1, cf1, cg1, wtmode, iorder, tlim, tolnl, nmesh, yn, rwsav, lrwsav, &comm, &fail); } else { tlim = 5.0; iskip = 7; h = tlim/(double) (nmesh - 1); wtmode = Nag_ReuseWeights; yn[0] = 1.0; /* nag_inteq_abel1_weak (d05bec) as above. */ nag_inteq_abel1_weak(ck2, cf2, cg2, wtmode, iorder, tlim, tolnl, nmesh, yn, rwsav, lrwsav, &comm, &fail); } if (fail.code != NE_NOERROR) { printf("Error from nag_inteq_abel1_weak (d05bec).\n%s\n", fail.message); exit_status = 1; goto END; } printf("The stepsize h = %8.4f\n\n", h); printf(" t Approximate\n"); printf(" Solution\n\n"); errmax = 0.0; for (i = 0; i < nmesh; i++) { hi1 = (double) (i) * h; err = fabs(yn[i] - ((exno==1)?sol1(hi1):sol2(hi1))); if (err > errmax) { errmax = err; t = hi1; soln = yn[i]; } if (i > 0 && i%iskip == 0) printf("%8.4f%15.4f\n", hi1, yn[i]); } printf("\nThe maximum absolute error, %10.2e, occurred at t = %8.4f\n", errmax, t); printf("with solution %8.4f\n", soln); } END: if (rwsav) NAG_FREE(rwsav); if (yn) NAG_FREE(yn); return exit_status; } static double sol1(double t) { return (1.0 / (sqrt(2.0 * nag_pi) * pow(t, 1.5))) * exp(-pow(1.0 + t, 2) / (2.0 * t)); } static double sol2(double t) { return 1.0/(1.0 + t); } static double NAG_CALL ck1(double t, Nag_Comm *comm) { return exp(-0.5 * t); } static double NAG_CALL cf1(double t, Nag_Comm *comm) { return (-1.0 / sqrt(nag_pi * t)) * exp(-0.5 * pow(1.0 + t, 2) / t); } static double NAG_CALL cg1(double s, double y, Nag_Comm *comm) { return y; } static double NAG_CALL ck2(double t, Nag_Comm *comm) { return sqrt(nag_pi); } static double NAG_CALL cf2(double t, Nag_Comm *comm) { /* Scalars */ double st1; st1 = sqrt(1.0 + t); return -2.0 * log(st1 + sqrt(t)) / st1; } static double NAG_CALL cg2(double s, double y, Nag_Comm *comm) { return y; }