/* pde_nl13.c  solve non linear second order, third degree equation
 * solve  u^2 u'' + 2 u u' + 3u = f(x)
 *        f(x) = (41/3)-(2116/9)x+(4780/3)x^2-(48976/9)x^3+
 *               9984x^4-(88064/9)x^5+(14336/3)x^6-(8192/9)x^7
 * with boundary  u(0)=1, u(1)=1  0 < x < 1 
 * set up system of equations, use specialized Jacobian, simeq_newton3
 * solution u(x) = 1.0 -(20.0/3.0)*x +12.0*x*x -(16.0/3.0)*x*x*x
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "deriv.h"
#include "simeq_newton3.h"
         /* also needs invert.h and invert.c */

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

#define nx 5
/* number of equations and number of linear unknowns */
#define neqn nx
/* number of nonlinear terms */
#define nlin 15
/* total number of terms */
#define ntot neqn+nlin
static double k[neqn][ntot]; /* nonlinear system of equations */
static double xg[nx];    /* xgrid */
static double u[ntot];   /* unknown, except for boundary values */
static double f[neqn];   /* RHS of nonlinear equations */
static double Ua[neqn];  /* only need linear values, other terms computed */
static double cx[nx];    /* numerical first derivative */
static double cxx[nx];   /* numerical second derivative */
static double xmin = 0.0;
static double xmax = 1.0;
static double hx;

/* linear, then nonlinear terms */
static int var1[ntot] =
           { 0, 1, 2, 3, 4, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3};
static int var2[ntot] =
           {-1,-1,-1,-1,-1, 1, 2, 3, 2, 3, 1, 1, 2, 3, 1, 2, 3, 1, 2, 3};
static int var3[ntot] =
           {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 2, 2, 2, 3, 3, 3};


static double F(double x) /* forcing function */
{
  return (41.0/3.0)-(2116.0/9.0)*x +(4780.0/3.0)*x*x -
         (48976.0/9.0)*x*x*x + 9984.0*x*x*x*x -
         (88064.0/9.0)*x*x*x*x*x + (14336.0/3.0)*x*x*x*x*x*x -
         (8192.0/9.0)*x*x*x*x*x*x*x;
} /* end F */

/* exact solution for checking */
static double uana(double x)
{
  return 1.0 -(20.0/3.0)*x +12.0*x*x -(16.0/3.0)*x*x*x;
} /* end uana */

static void check_soln(double U[]) /* check when solution unknown */
{
  /* u^2 u'' + 2 u u' + 3u = f(x) check for each xg inside boundary */
  int i, j;
  double xi, val, err, smaxerr;

  printf("check_soln \n");
  smaxerr = 0.0;
  for(i=1; i<neqn-1; i++)
  {
    xi = xg[i];
    rderiv(1, nx, i, hx, cx);
    rderiv(2, nx, i, hx, cxx);
    err = 3.0*U[i]-F(xi);
    for(j=0; j<nx; j++)
    {
      err += U[i]*U[i]*cxx[j]*U[j];
      err += 2.0*U[i]*cx[j]*U[j];
    }
    smaxerr = max(smaxerr, abs(err));
    printf("i=%2d, err=%g \n", i, err);
  }
  printf("check_soln max error=%g \n\n", smaxerr);
} /* end check_soln */

int main(int argc, char *argv[])
{
  int i, j;
  double xi, xj;
  double val, err, avgerr, maxerr;

  printf("pde_nl13.c running second order, third degree nonlinear \n");
  printf("solve x^2 u''(x) + 2x u'(x) + 3u(x) = f(x)\n");
  printf("  f(x) = (41.0/3.0)-(2116.0/9.0)*x +(4780.0/3.0)*x*x - \n");
  printf("         (48976.0/9.0)*x*x*x + 9984.0*x*x*x*x - \n");
  printf("         (88064.0/9.0)*x*x*x*x*x + (14336.0/3.0)*x*x*x*x*x*x -\n");
  printf("         (8192.0/9.0)*x*x*x*x*x*x*x) \n");
  printf(" \n");
  printf("  0<x<1  Boundary u(0) = 1, u(1) = 1 \n");
  printf("Analytic solution u(x) = 1 -(20/3)*x +12*x^2 -(16/3)*x^3 \n");

  hx = (xmax-xmin)/(double)(nx-1);
  printf("nx=%d, xmin=%f, xmax=%f, hx=%f \n",
         nx, xmin, xmax, hx);
  printf("F forcing function and exact solution \n");
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx; /* does not need to be uniform */
    Ua[i] = uana(xg[i]);
    printf("i=%d, F(%7.3f)=%7.3f, exact=%7.3f \n",
           i, xg[i], F(xg[i]), Ua[i]);
  }  
  printf("\n");
  
  /* initialize k and f */
  for(i=0; i<neqn; i++)
  {
    for(j=0; j<ntot; j++)
    {
      k[i][j] = 0.0;
    }
    f[i] = 0.0;
  }

  /* set boundary conditions */
  k[0][0] = 1.0;
  printf("k[%2d][%2d]=%g \n", 0, 0, k[0][0]);
  u[0] = uana(xmin);
  f[0] = u[0];
  printf("F(%f)=%g, for f at i=%d \n", xmin, f[0], 0);
  k[neqn-1][neqn-1] = 1.0;
  printf("k[%2d][%2d]=%g \n", neqn-1, neqn-1, k[neqn-1][neqn-1]);
  u[neqn-1] = uana(xmax);
  f[neqn-1] = u[neqn-1];
  printf("F(%f)=%g, for f at i=%d \n", xmax, f[neqn-1], neqn-1);

  /* compute entries in stiffness matrix */
  /* method
   *  u^2 u'' + 2 u u' + 3u = f(x)  becomes system of non linear equations
   *  at xg[0], xg[1], ... , xg[nx-1]  denoted x0, x1, ..., x4
   *  special case x0 and x4 are boundary and thus value is known
   *  thus U(x0) and U(x4) are known, U(x1) U(x2) and U(x3) are linear DOF
   *
   *U(x1)*U(x1)*cxx[0]*U(x0)+U(x1)*U(x1)*cxx[1]*U(x1)+U(x1)*U(x1)*cxx[2]*U(x2)+
   *U(x1)*U(x1)*cxx[3]*U(x3)+U(x1)*U(x1)*cxx[4]*U(x4)+
   *2*U(x1)*cx[0]*U(x0)+2*U(x1)*cx[1]*U(x1)+2*U(x1)*cx[2]*U(x2)+
   *2*U(x1)*cx[3]*U(x3)+2*U(x1)*cx[4]*U(x4)+
   *3*U(x1)=f(x1)
   *
   *U(x2)*U(x2)*cxx[0]*U(x0)+U(x2)*U(x2)*cxx[1]*U(x1)+U(x2)*U(x2)*cxx[2]*U(x2)+
   *U(x2)*U(x2)*cxx[3]*U(x3)+U(x2)*U(x2)*cxx[4]*U(x4)+
   *2*U(x2)*cx[0]*U(x0)+2*U(x2)*cx[1]*U(x1)+2*U(x2)*cx[2]*U(x2)+
   *2*U(x2)*cx[3]*U(x3)+2*U(x2)*cx[4]*U(x4)+
   *3*U(x2)=f(x2)
   *
   *U(x3)*U(x3)*cxx[0]*U(x0)+U(x3)*U(x3)*cxx[1]*U(x1)+U(x3)*U(x3)*cxx[2]*U(x2)+
   *U(x3)*U(x3)*cxx[3]*U(x3)+U(x3)*U(x3)*cxx[4]*U(x4)+
   *2*U(x3)*cx[0]*U(x0)+2*U(x3)*cx[1]*U(x1)+2*U(x3)*cx[2]*U(x2)+
   *2*U(x3)*cx[3]*U(x3)+2*U(x3)*cx[4]*U(x4)+
   *3*U(x3)=f(x3)
   *
   * The linear and non linear unknown terms are:
   * U(x1)              U(x2)              U(x3)  
   * U(x1)*U(x1)        U(x2)*U(x2)        U(x3)*U(x3)
   * U(x1)*U(x2)        U(x2)*U(x3)        U(x3)*U(x1)
   * U(x1)*U(x1)*U(x1)  U(x2)*U(x2)*U(x1)  U(x3)*U(x3)*U(x1)
   * U(x1)*U(x1)*U(x2)  U(x2)*U(x2)*U(x2)  U(x3)*U(x3)*U(x2)
   * U(x1)*U(x1)*U(x3)  U(x2)*U(x2)*U(x3)  U(x3)*U(x3)*U(x3)
   *
   * The coefficients of the above terms are:
   * 3+                 3+                 3+
     2*cx[0]*U(x0)+
     2*cx[4]*U(x4)
   * cxx[0]*U(x0)+      cxx[0]*U(x0)+      cxx[0]*U(x0)+
     cxx[4]*U(x4)+      cxx[0]*U(x4)       cxx[0]*U(x4)
     2*cx[1]            2*cx[2]            2*cx[3]
   * 2*cx[2]+2*cx[1]    2*cx[3]+2*cx[2]    2*cx[1]+2*cx[3]
   * cxx[1]             cxx[1]             cxx[1]
   * cxx[2]             cxx[2]             cxx[2]
   * cxx[3]             cxx[3]             cxx[3]
   *
   * var1[] = { 0, 1, 2, 3, 4, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3};
   * var2[] = {-1,-1,-1,-1,-1, 1, 2, 3, 2, 3, 1, 1, 2, 3, 1, 2, 3, 1, 2, 3};
   * var3[] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 2, 2, 2, 3, 3, 3};
   * j index    0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
   */


  printf("compute non linear stiffness matrix \n");
  i=1; /* first=0 and last=neqn-1 are boundary */
  xi = xg[i];
  rderiv(1, nx, i, hx, cx); /* discrete derivative coefficients */
  rderiv(2, nx, i, hx, cxx);
  /* compute val for k(i,j)  */
  j = 1;  /* U(x1) term */
  k[i][j] += 3.0 + 2.0*cx[0]*u[0] + 2.0*cx[neqn-1]*u[neqn-1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 5;  /* U(x1)*U(x1) term */
  k[i][j] += cxx[0]*u[0] + cxx[neqn-1]*u[neqn-1] + 2.0*cx[1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 8;  /* U(x1)*U(x2) term */
  k[i][j] += 2.0*cx[2];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 10; /* U(x3)*U(x1) term */
  k[i][j] += 2.0*cx[3];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 11;  /* U(x1)*U(x1)*U(x1) term */
  k[i][j] += cxx[1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 14; /* U(x1)*U(x1)*U(x2) term */
  k[i][j] += cxx[2];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 17; /* U(x1)*U(x1)*U(x3) term */
  k[i][j] += cxx[3];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  f[i] = F(xi); /* fill in RHS for f(i)  */
  printf("F(%f)=%g, for f at i=%d \n", xi, f[i], i);


  i=2; /* actually the third equation */
  xi = xg[i];
  rderiv(1, nx, i, hx, cx); /* discrete derivative coefficients */
  rderiv(2, nx, i, hx, cxx);
  /* compute val for k(i,j)  */
  j = 2;  /* U(x2) term */
  k[i][j] += 3.0 + 2.0*cx[0]*u[0] + 2.0*cx[neqn-1]*u[neqn-1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 6;  /* U(x2)*U(x2) term */
  k[i][j] += cxx[0]*u[0] + cxx[neqn-1]*u[neqn-1] + 2.0*cx[2];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 8;  /* U(x2)*U(x1) term */
  k[i][j] += 2.0*cx[1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 9;  /* U(x2)*U(x3) term */
  k[i][j] += 2.0*cx[3];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 12; /* U(x2)*U(x2)*U(x1) term */
  k[i][j] += cxx[1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 15; /* U(x2)*U(x2)*U(x2) term */
  k[i][j] += cxx[2]; 
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 18; /* U(x2)*U(x2)*U(x3) term */
  k[i][j] += cxx[3];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  f[i] = F(xi); /* fill in RHS for f(i)  */
  printf("F(%f)=%g, for f at i=%d \n", xi, f[i], i);


  i=3; /* first=0 and last=neqn-1 are boundary */
  xi = xg[i];
  rderiv(1, nx, i, hx, cx); /* discrete derivative coefficients */
  rderiv(2, nx, i, hx, cxx);
  /* compute val for k(i,j)  */
  j = 3;  /* U(x3) term */
  k[i][j] += 3.0 + 2.0*cx[0]*u[0] + 2.0*cx[neqn-1]*u[neqn-1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 7;  /* U(x3)*U(x3) term */
  k[i][j] += cxx[0]*u[0] + cxx[neqn-1]*u[neqn-1] + 2.0*cx[3];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 9;  /* U(x3)*U(x2) term */
  k[i][j] += 2.0*cx[2];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 10; /* U(x3)*U(x1) term */
  k[i][j] += 2.0*cx[1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 13; /* U(x3)*U(x3)*U(x1) term */
  k[i][j] += cxx[1];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 16; /* U(x3)*U(x3)*U(x2) term */
  k[i][j] += cxx[2];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  j = 19; /* U(x3)*U(x3)*U(x3) term */
  k[i][j] += cxx[3];
  printf("k[%2d][%2d]=%g \n", i, j, k[i][j]);
  f[i] = F(xi); /* fill in RHS for f(i)  */
  printf("F(%f)=%g, for f at i=%d \n", xi, f[i], i);

  printf("k computed non linear stiffness matrix, see above \n");
  printf("f computed forcing function, see above \n");
  printf("\n");

  /* guess at initial conditions, boundary known */
  u[0] = uana(xmin); /* also set up previously */
  u[neqn-1] = uana(xmax);
  for(i=1; i<neqn-1; i++) u[i] = 1.0;

  simeq_newton3(neqn, nlin, k, f, var1, var2, var3, u, 6, 0.001, 1);

  avgerr = 0.0;
  maxerr = 0.0;
  printf("u computed solution, Ua analytic, error \n");
  for(i=0; i<nx; i++)
  {
    err = abs(u[i]-Ua[i]);
    if(err>maxerr) maxerr = err;
    avgerr = avgerr + err;
    printf("u[%d]=%7.4f, Ua=%7.4f, err=%g \n",
            i, u[i], Ua[i], u[i]-Ua[i]);
  }
  printf("                                        maxerr=%g, avgerr=%g \n",
         maxerr, avgerr/(double)nx);
  printf("\n");

  printf("check_soln on computed solution against PDE \n");
  check_soln(u);  /* for use when analytic solution is not known */
  printf("just checking check_soln when given correct solution \n");
  check_soln(Ua); /* checking on checker */

  /* try giving the correct solution, to check on solver */
  printf("just checking simeq_newton3 when given correct solution \n");
  u[0] = 1.0;
  u[1] = 0.0;
  u[2] = 0.0;
  u[3] = 0.5;
  u[4] = 1.0;
  simeq_newton3(neqn, nlin, k, f, var1, var2, var3, u, 3, 0.001, 1);
  
  printf("pde_nl13.c finished \n");
  return 0;
} /* end main of pde_nl13.c */