# beam2.py  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.
 
from numpy import array
from numpy.linalg import solve
from numpy.matlib import zeros

# problem definitions
E = 10.0; # constants, values and units can be found
I =  9.0;
W =  7.0;
L =  8.0;
h = 0.0;
# the right hand side of d^4 y / dx^4 
def F1(x, y): # uniform load
    return -W/(L*E*I)

def difeq():   # difference equation solution  
  n = 7 # number of equations, zeros on front and back 
  A = array([[0.0 for j in range(n)] for i in range(n)])
  U = array([0.0 for i in range(n)])
  F = array([0.0 for i in range(n)])  # A U = F  F is p(x) 
                # possibly corrected for boundaries 
  C = array([1.0, -4.0, 6.0, -4.0, 1.0])
  #
  #  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
  # 
  h = L/(n+1) # must count h's for ends 
  print "Fourth order difference equation solution"
  print "boundary dy/dx=0 fourth order"
  print "n=%d, h=%g" %(n,h)
  print "C={%g,%g,%g,%g,%g}" %(C[0],C[1],C[2],C[3],C[4])
  for i in range(n):
    for j in range(n):
      A[i][j] = 0.0
  for i in range(n):
     F[i] = -W*h*h*h*h/(L*E*I)
  # handle first two and last two rows with boundary conditions 
  A[0][0] =  48.0
  A[0][1] = -36.0
  A[0][2] =  16.0
  A[0][3] =  -3.0
  F[0]    =   0.0*12.0*h # Neumann 
  A[1][0] = -4.0
  A[1][1] =  6.0
  A[1][2] = -4.0
  A[1][3] =  1.0    # Dirichlet 
  A[5][3] =  1.0
  A[5][4] = -4.0
  A[5][5] =  6.0
  A[5][6] = -4.0    # Dirichlet 
  A[6][3] =  -3.0
  A[6][4] =  16.0
  A[6][5] = -36.0
  A[6][6] =  48.0
  F[6]    =   0.0*12.0*h # Neumann 
  #
  for i in range(2,5): # all inner rows 
    for j in range(i-2,i-2+5):
      A[i][j] = C[j-i+2]


  # debug print of simultaneous equations to be solved 
  print "solve A U = F for U "
  for i in range(n):
   for j in range(n):
     print "A[%d][%d]=%g" %(i,j,A[i][j])

   print "F[%d]=%g" %(i,F[i])


  U = solve(A, F)
  print "fourth order difference equations are exact within"
  print "numerical accuracy when solution is fourth order."
  print "x=%g, y=%g  boundary" %(0.0,0.0)
  for i in range(n):
    x = (i+1)*h
    y = U[i]
    print "x=%g, y=%g" %(x,y)

  print "x=%g, y=%g  boundary" %((n+1)*h,0.0) 
  # end difeq 


print "beam2.py running"
n = 9
h = L/(n-1)

print "E=%g, I=%g, W=%g, L=%g, h=%g, n=%d" %(E,I,W,L,h,n)
print "Analytic Solution."
for i in range(n):
  x = i*h
  y = (-W* x*x *(x-L)*(x-L))/(24.0*L*E*I)
  print "beam x=%g, y=%g" %(x,y)
  slope = ((-W*x*x*x)/(6.0*L) + (W*x*x)/4.0 - (W*L*x)/12.0)/(E*I)
  print "                       slope=%g" %(slope)

print "Numerical Solution."
# not, fourth order Runge-Kutta 
difeq() 
print "beam2.py finished"
# end  beam2.py