NAG Library Manual, Mark 30
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_ode_bvp_coll_nlin_diag (d02tzc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.0, 2024.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

typedef struct {
  double alpha, beta, eps;
  Integer mmax;
} func_data;

#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, neq = 1, mmax = 2, nlbc = 1, nrbc = 1;
  Integer i, iermx, ijermx, j, licomm, lrcomm, mxmesh, ncol, nmesh;
  double alpha, beta, eps, ermx;
  /* Arrays */
  static double ruser[7] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
  double *mesh = 0, *rcomm = 0, *tol = 0, *y = 0;
  double rdum[1];
  Integer *ipmesh = 0, *icomm = 0, *m = 0;
  Integer idum[2];
  /* Nag Types */
  Nag_Boolean failed = Nag_FALSE;
  func_data fd;
  Nag_Comm comm;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_ode_bvp_coll_nlin_diag (d02tzc) Example Program Results\n\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  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;
  }
  /* Set problem orders */
  m[0] = 2;
  scanf("%lf%lf%lf%*[^\n] ", &alpha, &beta, &eps);
  for (i = 0; i < nmesh; i++) {
    scanf("%lf", &mesh[i]);
  }
  scanf("%*[^\n] ");
  for (i = 0; i < nmesh; i++) {
    scanf("%" NAG_IFMT "", &ipmesh[i]);
  }
  scanf("%*[^\n] ");
  for (i = 0; i < neq; i++) {
    scanf("%lf", &tol[i]);
  }
  scanf("%*[^\n] ");

  /* 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;
  }

  eps = 0.1 * eps;
  /* Set data required for the user-supplied functions */
  fd.alpha = alpha;
  fd.beta = beta;
  fd.eps = eps;
  fd.mmax = mmax;
  /* Associate the data structure with comm.p */
  comm.p = (Pointer)&fd;

  for (j = 0; j < 2; j++) {
    printf("\n Tolerance = %8.1e  eps = %10.3e\n", tol[0], eps);
    /* 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;
      /* Continue and print the mesh statistics regardless. */
    }

    /* 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", ijermx);
    if (failed) {
      goto END;
    }

    /* Print solution at every second point on final mesh. */
    printf("\n Solution and derivative at every second point:\n");
    printf("   x           u          u'\n");
    for (i = 0; i < nmesh; i += 2) {

      /* 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(mesh[i], 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.4f %11.5f %11.5f \n", mesh[i], y[0], y[neq]);
    }
    if (j == 0) {
      /* Halve final mesh for new initial mesh and set up for continuation. */
      nmesh = (nmesh + 1) / 2;

      /* nag_ode_bvp_coll_nlin_contin (d02txc).
       * Ordinary differential equations, general nonlinear boundary value
       * problem, continuation facility for
       * nag_ode_bvp_coll_nlin_solve (d02tlc).
       */
      nag_ode_bvp_coll_nlin_contin(mxmesh, nmesh, mesh, ipmesh, rcomm, icomm,
                                   &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_ode_bvp_coll_nlin_contin (d02txc).\n%s\n",
               fail.message);
        exit_status = 4;
        goto END;
      }

      /* Reduce continuation parameter. */
      eps = 0.1 * eps;
      fd.eps = eps;
    }
  }

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) {
  func_data *fd = (func_data *)comm->p;

  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback ffun, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  f[0] = (y[0] - y[0] * y[1]) / fd->eps;
}

static void NAG_CALL fjac(double x, const double y[], Integer neq,
                          const Integer m[], double dfdy[], Nag_Comm *comm) {
  double epsh, fac, ptrb;
  Integer j;
  double f1[1], f2[1], yp[2];

  /* The size of yp here equals m[0]. */

  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback fjac, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  /* nag_machine_precision (x02ajc).
   * The machine precision.
   */
  epsh = 100.0 * nag_machine_precision;
  fac = sqrt(nag_machine_precision);
  for (j = 0; j < 2; j++) {
    yp[j] = y[j];
  }
  for (j = 0; j < 2; j++) {
    ptrb = MAX(epsh, fac * fabs(y[j]));
    yp[j] = y[j] + ptrb;
    ffun(x, yp, neq, m, f1, comm);
    yp[j] = y[j] - ptrb;
    ffun(x, yp, neq, m, f2, comm);
    dfdy[j] = 0.5 * (f1[0] - f2[0]) / ptrb;
    yp[j] = y[j];
  }
}

static void NAG_CALL gafun(const double ya[], Integer neq, const Integer m[],
                           Integer nlbc, double ga[], Nag_Comm *comm) {
  func_data *fd = (func_data *)comm->p;

  if (comm->user[2] == -1.0) {
    printf("(User-supplied callback gafun, first invocation.)\n");
    comm->user[2] = 0.0;
  }
  ga[0] = ya[0] - fd->alpha;
}

static void NAG_CALL gbfun(const double yb[], Integer neq, const Integer m[],
                           Integer nrbc, double gb[], Nag_Comm *comm) {
  func_data *fd = (func_data *)comm->p;

  if (comm->user[3] == -1.0) {
    printf("(User-supplied callback gbfun, first invocation.)\n");
    comm->user[3] = 0.0;
  }
  gb[0] = yb[0] - fd->beta;
}

static void NAG_CALL gajac(const double ya[], Integer neq, const Integer m[],
                           Integer nlbc, double dgady[], Nag_Comm *comm) {
  if (comm->user[4] == -1.0) {
    printf("(User-supplied callback gajac, first invocation.)\n");
    comm->user[4] = 0.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->user[5] == -1.0) {
    printf("(User-supplied callback gbjac, first invocation.)\n");
    comm->user[5] = 0.0;
  }
  dgbdy[0] = 1.0;
}

static void NAG_CALL guess(double x, Integer neq, const Integer m[], double y[],
                           double dym[], Nag_Comm *comm) {
  func_data *fd = (func_data *)comm->p;
  double alpha, beta;

  if (comm->user[6] == -1.0) {
    printf("(User-supplied callback guess, first invocation.)\n");
    comm->user[6] = 0.0;
  }
  alpha = fd->alpha;
  beta = fd->beta;
  y[0] = alpha + (beta - alpha) * x;
  y[1] = beta - alpha;
  dym[0] = 0.0;
}