// beam2.java  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.
//
//  needs simeq.java 
 
class beam2
{
  // problem definitions
  double E = 10.0; // constants, values and units can be found
  double I =  9.0;
  double W =  7.0;
  double L =  8.0;
  double h;
  // the right hand side of d^4 y / dx^4 
  double F1(double x, double y) // uniform load
  {
    return -W/(L*E*I);
  } // end F1

  beam2()
  {
    double x, y, slope;
    double k1,k2,k3,k4;
    double y4, y3, y2, y1; 
    int i, n;

    System.out.println("beam2.java running ");
    n = 9;
    h = L/(double)(n-1);
    System.out.println("E="+E+", I="+I+", W="+W+", L="+L+", h="+h+", n="+n);
    System.out.println("Analytic Solution.");
    for(i=0; i<n; i++)
    {
      x = (double)i*h;
      y = (-W* x*x *(x-L)*(x-L))/(24.0*L*E*I);
      System.out.println("beam x="+x+", y="+y);
      slope = ((-W*x*x*x)/(6.0*L) + (W*x*x)/4.0 - (W*L*x)/12.0)/(E*I);
      System.out.println("                       slope"+slope); 
    }

   System.out.println("Numerical Solution.");
   // not, fourth order Runge-Kutta 

   difeq(); // below 
    System.out.println("beam2.java finished");
  } // end main of beam2

  void difeq()
  { // difference equation solution  
    double A[] = new double[10000];
    double U[] = new double[100];
    double F[] = new double[100];  // A U = F  F is p(x) 
                                   // possibly corrected for boundaries 
    double C[] = {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 
    System.out.println("Fourth order difference equation solution ");
    System.out.println("boundary dy/dx=0 fourth order");
    System.out.println("n="+n+", h="+h);
    System.out.println("C={"+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 
    System.out.println("solve A y = F for y ");
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
        System.out.println("A["+i+","+j+"]="+A[i*n+j]);
      System.out.println("F["+i+"]="+F[i]);
    }


    new simeq(n, A, F, U); // U have y values
    System.out.println("fourth order difference equations are exact within");
    System.out.println("numerical accuracy when solution is fourth order.");
    System.out.println("x="+0.0+", y="+0.0+"  boundary");
    for(i=0; i<n; i++)
    {
      x = (double)(i+1)*h;
      y = U[i];
      System.out.println("x="+x+", y="+y);
    }
    System.out.println("x="+((double)(n+1)*h)+", y="+0.0+"  boundary"); 
  } // end difeq 

  public static void main (String[] args)
  {
    new beam2();
  } // end main
} // end beam2.java