/* nag_inteq_abel1_weak (d05bec) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.4, 2022.
*/
#include <math.h>
#include <nag.h>
#ifdef __cplusplus
extern "C" {
#endif
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
static double sol1(double t);
static double sol2(double t);
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 */
static double ruser[6] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
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");
/* For communication with user-supplied functions: */
comm.user = ruser;
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 %" NAG_IFMT "\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;
t = 0.0;
soln = 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:
NAG_FREE(rwsav);
NAG_FREE(yn);
return exit_status;
}
static double sol1(double t) {
if (t == 0.0)
return 0.0;
else
return (1.0 / (sqrt(2.0 * nag_math_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) {
if (comm->user[0] == -1.0) {
printf("(User-supplied callback ck1, first invocation.)\n");
comm->user[0] = 0.0;
}
return exp(-0.5 * t);
}
static double NAG_CALL cf1(double t, Nag_Comm *comm) {
if (comm->user[1] == -1.0) {
printf("(User-supplied callback cf1, first invocation.)\n");
comm->user[1] = 0.0;
}
return (-1.0 / sqrt(nag_math_pi * t)) * exp(-0.5 * pow(1.0 + t, 2) / t);
}
static double NAG_CALL cg1(double s, double y, Nag_Comm *comm) {
if (comm->user[2] == -1.0) {
printf("(User-supplied callback cg1, first invocation.)\n");
comm->user[2] = 0.0;
}
return y;
}
static double NAG_CALL ck2(double t, Nag_Comm *comm) {
if (comm->user[3] == -1.0) {
printf("(User-supplied callback ck2, first invocation.)\n");
comm->user[3] = 0.0;
}
return sqrt(nag_math_pi);
}
static double NAG_CALL cf2(double t, Nag_Comm *comm) {
/* Scalars */
double st1;
if (comm->user[4] == -1.0) {
printf("(User-supplied callback cf2, first invocation.)\n");
comm->user[4] = 0.0;
}
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) {
if (comm->user[5] == -1.0) {
printf("(User-supplied callback cg2, first invocation.)\n");
comm->user[5] = 0.0;
}
return y;
}