/* nag_inteq_volterra2 (d05bac) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 23, 2011.
*/
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd05.h>
#include <nagx02.h>
#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 */
static double ruser[3] = {-1.0, -1.0, -1.0};
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");
/* For communication with user-supplied functions: */
comm.user = ruser;
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];
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:
NAG_FREE(errest);
NAG_FREE(yn);
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)
{
if (comm->user[0] == -1.0)
{
printf("(User-supplied callback cf, first invocation.)\n");
comm->user[0] = 0.0;
}
return exp(-t);
}
static double NAG_CALL ck(double t, Nag_Comm *comm)
{
if (comm->user[1] == -1.0)
{
printf("(User-supplied callback ck, first invocation.)\n");
comm->user[1] = 0.0;
}
return exp(-t);
}
static double NAG_CALL cg(double s, double y, Nag_Comm *comm)
{
if (comm->user[2] == -1.0)
{
printf("(User-supplied callback cg, first invocation.)\n");
comm->user[2] = 0.0;
}
return y + exp(-y);
}