/* beam2.c  ordinary differential equation for a loaded beam */
/*          handle Neumann boundary  dy/dx=0 at x=0 and x=L better */

/*  given a beam of length L, from 0 < x < L
 *  with Young's Modulus of E
 *  with moment of inertia I
 *  with p(x) being the load density  e.g. force per unit length
 *  with both ends fixed, meaning:
 *       y=0 at x=0, y=0 at x=L, dy/dx=0 at x=0, dy/dx=0 at x=L
 *  then the differential equation that defines the y position of
 *  the beam at every x coordinate is
 *
 *      d^4 y
 *   EI ----- = p(x)   with the convention for downward is negative
 *      dx^4
 *
 *  for uniformly distributed force (weight) p(x) becomes  - W/L
 *  This simple case can be integrated and solved analytically:
 *
 *      d^4 y 
 *   EI ----- = -W/L
 *      dx^4
 *
 *     d^3 y
 *  EI ----- = -W/L x + A   ( A is constant of integration, value found later)
 *     dx^3
 *
 *     d^2 y
 *  EI ----- = -W/L x^2/2 + A x + B 
 *     dx^2
 *
 *     dy
 *  EI -- = -W/L x^3/6 + A x^2/2 + B x + C  we see C=0 because dy/dx=0 at x=0
 *     dx
 *
 *  EI y = -W/L x^4/24 + A x^3/6 + B x^2/2 + D  we see D=0 because y=0 at x=0
 *
 *  substituting y=0 at x=L in the equation above gives
 *
 *     0 = -W/L L^4/24 + A L^3/6 + B L^2/2
 *
 *  substituting dy/dx=0 at x=L in the equation above the above gives
 *
 *     0 = -W/L L^3/6 + A L^2/2 + B L
 *
 *  solving two equations in two unknowns  A = W/2  B = - WL/12
 *  then substituting for A and B in  EI y = ... and combining terms
 *
 *         W
 *  y = --------  x^2 (x-L)^2
 *      24 L EI
 *
 *  The known solution for a specific case is valuable to check your
 *  programming of a numerical solution before computing a more
 *  general case of p(x) where no closed form solution may exists.
 */

#include <stdio.h>
#include <math.h>
#undef y1
#undef  abs
#define abs(x) (((x)<0.0)?(-(x)):(x))

/* the right hand side of d^4 y / dx^4  uniform load */
#define F1(x,y) (-W/(L*E*I))

void simeq(int n, double A[], double Y[], double X[]);

static double E = 10.0; /* arbitrary constants, values and units can be found*/
static double I =  9.0;
static double W =  7.0;
static double L =  8.0;

static void difeq();

int main(int argc, char *argv[])
{
  double x, y, h, slope;
  double k1,k2,k3,k4;
  double y4, y3, y2, y1; 
  int i, n;

  printf("beam2.c running \n");
  n = 9;
  h = L/(double)(n-1);
  printf("E=%g, I=%g, W=%g, L=%g, h=%g, n=%d \n",
         E, I, W, L, h, n);
  printf("\nAnalytic Solution.\n");
  for(i=0; i<n; i++)
  {
    x = (double)i*h;
    y = (-W* x*x *(x-L)*(x-L))/(24.0*L*E*I);
    printf("beam x=%g, y=%g \n", x, y);
    slope = ((-W*x*x*x)/(6.0*L) + (W*x*x)/4.0 - (W*L*x)/12.0)/(E*I);
    printf("                       slope=%g \n", slope); 
  }

  printf("\nNumerical Solution.\n");
  /* not, fourth order Runge-Kutta */

  difeq(); /* below */
  printf("\nbeam2.c finished\n");
  return 0;
} /* end main of beam2.c */

static void difeq()
{ /* difference equation solution  */
  double A[10000];
  double U[100];
  double F[100];  /* A U = F  F is p(x) possibly corrected for boundaries */
  double C[5] = {1.0, -4.0, 6.0, -4.0, 1.0};
  double x, y, h;
  int i, j, k, n;

  /*  d^4 y / dx^4 as a difference equation =
   *  1/h^4 ( 1 y(x-2h) - 4 y(x-h) + 6 y(x) - 4 y(x+h) + 1 y(x+2h) )
   *  1/h^4 moved to other side of equations (inner)
   *
   *  dy/dx at x = 1/12h ( -25 y(x) +48 y(x+h) -36 y(x+2h) +16 y(x+3h) -3 y(x+4h) )
   *  multiply both sides by 12h, but constant side is zero for this problem
   *  y(0) = 0 for this problem thus no constant contribution
   *  flip for x=L  -3, +16, -36, +48   y(L)=0 thus no -25
   */
  n = 7; /* number of equations, zeros on front and back */
  h = L/(double)(n+1); /* must count h's for ends */
  printf("Fourth order difference equation solution \n");
  printf("boundary dy/dx=0 fourth order\n");
  printf("n=%d, h=%g, C={%g,%g,%g,%g,%g}\n",
         n, h, C[0], C[1], C[2], C[3], C[4]);
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      A[i*n+j] = 0.0;
  for(i=0; i<n; i++) F[i] = -W*h*h*h*h/(L*E*I);
  /* handle first two and last two rows with boundary conditions */
  i=0;   j=0;   A[i*n+j] =  48.0;
         j=1;   A[i*n+j] = -36.0;
         j=2;   A[i*n+j] =  16.0;
         j=3;   A[i*n+j] =  -3.0;
	        F[i]     =  0.0*12.0*h; /* Neumann */
  i=1;   j=0;   A[i*n+j] = -4.0;
         j=1;   A[i*n+j] =  6.0;
         j=2;   A[i*n+j] = -4.0;
         j=3;   A[i*n+j] =  1.0;    /* Dirichlet */
  i=n-2; j=n-4; A[i*n+j] =  1.0;
         j=n-3; A[i*n+j] = -4.0;
         j=n-2; A[i*n+j] =  6.0;
         j=n-1; A[i*n+j] = -4.0;    /* Dirichlet */
  i=n-1; j=n-4; A[i*n+j] =  -3.0;
         j=n-3; A[i*n+j] =  16.0;
         j=n-2; A[i*n+j] = -36.0;
         j=n-1; A[i*n+j] =  48.0;
	        F[i]     =  0.0*12.0*h; /* Neumann */

  for(i=2; i<n-2; i++) /* all inner rows */
    for(j=i-2; j<i-2+5; j++)
      A[i*n+j] = C[j-i+2];

  /* debug print of simultaneous equations to be solved */
  printf("solve A y = F for y \n");
  printf("                           A                            F \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      printf("%7.3f ", A[i*n+j]);
    printf("%g\n", F[i]);
  }

  simeq(n, A, F, U); /* U has y values */
  printf("\nfourth order difference equations are exact within\n%s\n",
         "numerical accuracy when solution is fourth order.");
  printf("x=%g, y=%g  boundary\n", 0.0, 0.0);
  for(i=0; i<n; i++)
  {
    x = (double)(i+1)*h;
    y = U[i];
    printf("x=%g, y=%g \n", x, y);
  }
  printf("x=%g, y=%g  boundary\n", (double)(n+1)*h, 0.0);

} /* end difeq */
/* end beam2.c */