/* nag_ode_bvp_coll_nlin_interp (d02tyc) Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.3, 2021.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL ffun(double x, const double y[], Integer neq,
const Integer m[], double f[], Nag_Comm *comm);
static void NAG_CALL fjac(double x, const double y[], Integer neq,
const Integer m[], double dfdy[], Nag_Comm *comm);
static void NAG_CALL gafun(const double ya[], Integer neq, const Integer m[],
Integer nlbc, double ga[], Nag_Comm *comm);
static void NAG_CALL gbfun(const double yb[], Integer neq, const Integer m[],
Integer nrbc, double gb[], Nag_Comm *comm);
static void NAG_CALL gajac(const double ya[], Integer neq, const Integer m[],
Integer nlbc, double dgady[], Nag_Comm *comm);
static void NAG_CALL gbjac(const double yb[], Integer neq, const Integer m[],
Integer nrbc, double dgbdy[], Nag_Comm *comm);
static void NAG_CALL guess(double x, Integer neq, const Integer m[], double y[],
double dym[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
int main(void) {
/* Scalars */
Integer exit_status = 0, nmesh_out = 11;
Integer neq, mmax, nlbc, nrbc;
Integer i, iermx, ijermx, licomm, lrcomm, mxmesh, ncol, nmesh;
double a, ainc, ermx, x;
/* Arrays */
static Integer iuser[7] = {-1, -1, -1, -1, -1, -1, -1};
double *mesh = 0, *rcomm = 0, *tol = 0, *y = 0;
double rdum[1], ruser[1];
Integer *ipmesh = 0, *icomm = 0, *m = 0;
Integer idum[2];
/* Nag Types */
Nag_Boolean failed = Nag_FALSE;
Nag_Comm comm;
NagError fail;
INIT_FAIL(fail);
printf("nag_ode_bvp_coll_nlin_interp (d02tyc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "%" NAG_IFMT "", &neq, &mmax);
scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nlbc, &nrbc);
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &ncol, &nmesh,
&mxmesh);
if (!(mesh = NAG_ALLOC(mxmesh, double)) || !(m = NAG_ALLOC(neq, Integer)) ||
!(tol = NAG_ALLOC(neq, double)) || !(y = NAG_ALLOC(neq * mmax, double)) ||
!(ipmesh = NAG_ALLOC(mxmesh, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
comm.user = ruser;
comm.iuser = iuser;
for (i = 0; i < neq; i++) {
scanf("%" NAG_IFMT "", &m[i]);
}
scanf("%*[^\n] ");
scanf("%lf%*[^\n] ", &a);
for (i = 0; i < neq; i++) {
scanf("%lf", &tol[i]);
}
scanf("%*[^\n] ");
ainc = a / (double)(nmesh - 1);
mesh[0] = 0.0;
for (i = 1; i < nmesh - 1; i++) {
mesh[i] = mesh[i - 1] + ainc;
}
mesh[nmesh - 1] = a;
ipmesh[0] = 1;
for (i = 1; i < nmesh - 1; i++) {
ipmesh[i] = 2;
}
ipmesh[nmesh - 1] = 1;
/* For communication with function guess store a in comm.user[0]. */
ruser[0] = a;
/* Communication space query to get size of rcomm and icomm
* by setting lrcomm=0 in call to
* nag_ode_bvp_coll_nlin_setup (d02tvc):
* Ordinary differential equations, general nonlinear boundary value problem,
* setup for nag_ode_bvp_coll_nlin_solve (d02tlc).
*/
nag_ode_bvp_coll_nlin_setup(neq, m, nlbc, nrbc, ncol, tol, mxmesh, nmesh,
mesh, ipmesh, rdum, 0, idum, 2, &fail);
if (fail.code == NE_NOERROR) {
lrcomm = idum[0];
licomm = idum[1];
if (!(rcomm = NAG_ALLOC(lrcomm, double)) ||
!(icomm = NAG_ALLOC(licomm, Integer))) {
printf("Allocation failure\n");
exit_status = -2;
goto END;
}
/* Initialize, again using nag_ode_bvp_coll_nlin_setup (d02tvc). */
nag_ode_bvp_coll_nlin_setup(neq, m, nlbc, nrbc, ncol, tol, mxmesh, nmesh,
mesh, ipmesh, rcomm, lrcomm, icomm, licomm,
&fail);
}
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_coll_nlin_setup (d02tvc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf(" Tolerance = %8.1e, a = %8.2f\n", tol[0], a);
/* Solve */
/* nag_ode_bvp_coll_nlin_solve (d02tlc).
* Ordinary differential equations, general nonlinear boundary value
* problem, collocation technique.
*/
nag_ode_bvp_coll_nlin_solve(ffun, fjac, gafun, gbfun, gajac, gbjac, guess,
rcomm, icomm, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_coll_nlin_solve (d02tlc).\n%s\n",
fail.message);
failed = Nag_TRUE;
}
/* Extract mesh. */
/* nag_ode_bvp_coll_nlin_diag (d02tzc).
* Ordinary differential equations, general nonlinear boundary value
* problem, diagnostics for nag_ode_bvp_coll_nlin_solve (d02tlc).
*/
nag_ode_bvp_coll_nlin_diag(mxmesh, &nmesh, mesh, ipmesh, &ermx, &iermx,
&ijermx, rcomm, icomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_coll_nlin_diag (d02tzc).\n%s\n",
fail.message);
exit_status = 2;
goto END;
}
/* Print mesh statistics */
printf("\n Used a mesh of %4" NAG_IFMT " points\n", nmesh);
printf(" Maximum error = %10.2e in interval %4" NAG_IFMT "", ermx, iermx);
printf(" for component %4" NAG_IFMT " \n\n\n", ijermx);
printf(" Mesh points:\n");
for (i = 0; i < nmesh; i++) {
printf("%4" NAG_IFMT "(%1" NAG_IFMT ")%13.4e%s", i + 1, ipmesh[i], mesh[i],
(i + 1) % 4 ? "" : "\n");
}
printf("\n");
if (!failed) {
/* Print solution on output mesh. */
printf("\n\n Computed solution\n x solution derivative\n");
x = 0.0;
ainc = a / (double)(nmesh_out - 1);
for (i = 0; i < nmesh_out; i++) {
/* nag_ode_bvp_coll_nlin_interp (d02tyc).
* Ordinary differential equations, general nonlinear boundary value
* problem, interpolation for nag_ode_bvp_coll_nlin_solve (d02tlc).
*/
nag_ode_bvp_coll_nlin_interp(x, y, neq, mmax, rcomm, icomm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_coll_nlin_interp (d02tyc).\n%s\n",
fail.message);
exit_status = 3;
goto END;
}
printf("%8.2f %11.5f %11.5f\n", x, y[0], y[neq]);
x = x + ainc;
}
}
END:
NAG_FREE(mesh);
NAG_FREE(m);
NAG_FREE(tol);
NAG_FREE(rcomm);
NAG_FREE(y);
NAG_FREE(ipmesh);
NAG_FREE(icomm);
return exit_status;
}
static void NAG_CALL ffun(double x, const double y[], Integer neq,
const Integer m[], double f[], Nag_Comm *comm) {
if (comm->iuser[0] == -1) {
printf("(User-supplied callback ffun, first invocation.)\n");
comm->iuser[0] = 0;
}
if (y[0] <= 0.0) {
f[0] = 0.0;
} else {
f[0] = pow(y[0], 1.5) / sqrt(x);
}
}
static void NAG_CALL fjac(double x, const double y[], Integer neq,
const Integer m[], double dfdy[], Nag_Comm *comm) {
if (comm->iuser[1] == -1) {
printf("(User-supplied callback fjac, first invocation.)\n");
comm->iuser[1] = 0;
}
if (y[0] <= 0.0) {
dfdy[0] = 0.0;
} else {
dfdy[0] = 1.5 * sqrt(y[0]) / sqrt(x);
}
}
static void NAG_CALL gafun(const double ya[], Integer neq, const Integer m[],
Integer nlbc, double ga[], Nag_Comm *comm) {
if (comm->iuser[2] == -1) {
printf("(User-supplied callback gafun, first invocation.)\n");
comm->iuser[2] = 0;
}
ga[0] = ya[0] - 1.0;
}
static void NAG_CALL gbfun(const double yb[], Integer neq, const Integer m[],
Integer nrbc, double gb[], Nag_Comm *comm) {
if (comm->iuser[3] == -1) {
printf("(User-supplied callback gbfun, first invocation.)\n");
comm->iuser[3] = 0;
}
gb[0] = yb[0];
}
static void NAG_CALL gajac(const double ya[], Integer neq, const Integer m[],
Integer nlbc, double dgady[], Nag_Comm *comm) {
if (comm->iuser[4] == -1) {
printf("(User-supplied callback gajac, first invocation.)\n");
comm->iuser[4] = 0;
}
dgady[0] = 1.0;
}
static void NAG_CALL gbjac(const double yb[], Integer neq, const Integer m[],
Integer nrbc, double dgbdy[], Nag_Comm *comm) {
if (comm->iuser[5] == -1) {
printf("(User-supplied callback gbjac, first invocation.)\n");
comm->iuser[5] = 0;
}
dgbdy[0] = 1.0;
}
static void NAG_CALL guess(double x, Integer neq, const Integer m[], double y[],
double dym[], Nag_Comm *comm) {
double a = comm->user[0];
if (comm->iuser[6] == -1) {
printf("(User-supplied callback guess, first invocation.)\n");
comm->iuser[6] = 0;
}
y[0] = 1.0 - x / (a);
y[neq] = -1.0 / (a);
dym[0] = 0.0;
}