Problem Set 1

Home Page

Problem 9.1:

(a)

Using the Galerkin Method to solve : ###\({d^2u \over {dt^2}} = v^2{d^2u \over {dx^2}} + \gamma {d \over dt}{d^2u \over {dx^2}}\)

Taking the Residual ###\(R= {d^2u \over {dt^2}} - v^2{d^2u \over {dx^2}} - \gamma {d \over dt}{d^2u \over {dx^2}}\)

The Galerkin method: ###\(\int _0 ^1 {{R \phi_j}}\,dx = 0\) where \(u(x,t) = \sum_{i} a_i(t) \phi_i(x_i)\)

The integral then becomes :

\(\sum_{i} \int _0 ^1 {({d^2u \over {dt^2}} - v^2{d^2u \over {dx^2}} - \gamma {d \over dt}{d^2u \over {dx^2}})\phi_j}\,dx = \sum_{i} \int _0 ^1 {({d^2a \over {dt^2}}\phi_i - v^2a{d^2\phi_i \over {dx^2}} - \gamma {da \over dt}{d^2\phi_i \over {dx^2}})\phi_j}\,dx\)

Carrying out Integration By Parts:

\(\sum_{i}[ {d^2a \over {dt^2}} \underbrace{\int _0 ^1 {\phi_i \phi_j }\,dx}_{A_{i j}} + av^2\underbrace{(\int _0 ^1 {{{d\phi_i \over dx} {d\phi_j \over dx} }}\,dx-[{d\phi_i \over dx} \phi_j]_0^1)}_{B_{i j}}+{da \over dt}\gamma\underbrace{(\int _0 ^1 {{{d\phi_i \over dx} {d\phi_j \over dx} }}\,dx-[{d\phi_i \over dx} \phi_j]_0^1)}_{B_{i j}} ]\)

Which is a system of second order ODE's of the form:

\(\bf{A}\) \({d^2a \over dt^2}\) + \({av^2}\) \(\bf{B}\) + \(\gamma {da \over dt}\) \(\bf{B}\) = 0

(b)

Linear hat expansion where the basis functions are hat functions \(\phi_i\) :

\[ \ \phi_i = \left\{% \begin{array}{11} {x-x_{i-1} \over {x_i-x_{i-1}}} &\textrm{if }x_{i-1} \le x < x_i\\ {x_{i+1}-x \over {x_{i+1}-x_{i}}} &\textrm{if }x_i \le x < x_{i+1}\\ 0 &\textrm{if } x<x_{i-1} \textrm{ or } x \ge x_{i+1} \end{array}% \right. \ \]

Considering the step size is h, then \(x_{i+1} = x_i + h\) and \(x_{i-1} = x_i - h\)

Calculating the coefficient terms:

\[ A_{i,i} = \int_{x_i-h}^{x_i}[{{x-x_i+h} \over h}]^2\,dx + \int_{x_i}^{x_{i+h}}[{{x_i+h-x} \over h}]^2\,dx \]

\[ A_{i,i-1} = \int_{x_i-h}^{x_i}[{{x_i-x}\over h}][{{x-(x_i-h)}\over h}]\,dx \]

Since the spatial distribution is equidistant and uniform of a step size h, then

\[ A_{i,i-1} = A_ {i,i+1}= A_{i-1,i} = A_ {i+1,i} \]

\[ B_{i,i} = \int_{x_i-h}^{x_i}[{{1} \over h}]^2\,dx + \int_{x_i}^{x_{i+h}}[{1 \over h}][{{x_{i}+h-1-0-x_{i}-h} \over h}]\,dx = \int_{x_i-h}^{x_i}[{{1} \over h}]^2\,dx + \int_{x_i}^{x_{i+h}}[{{1} \over h}]^2\,dx \]

\[ B_{i,i-1} = \int_{x_i-h}^{x_i}[{{1} \over h}][{{-1} \over h}]\,dx \]

In [1]: 
from sage.all import *
x= var('x')
xi=var('xi')
h=var('h')
Aii= integral(((x-xi+h)/h)**2,x, xi-h,xi)+integral(((xi+h-x)/h)**2,x, xi,xi+h)
print "Aii = " + str(Aii)

Aii1= integral(((xi-x)/h)*((x-xi+h)/h),x, xi-h,xi)
print "Ai,i-1 = " + str(Aii1)

Bii= integral((1/h**2),x,xi-h,xi)+ integral((1/h**2),x,xi,xi+h)
print "Bi,i = " + str(Bii)

Bii1= integral((-1/h**2),x,xi-h,xi)
print "Bi,i-1 = " + str(Bii1)

Aii = 2/3*h
Ai,i-1 = 1/6*h
Bi,i = 2/h
Bi,i-1 = -1/h

(c)

In [2]: 
# Inverting the matrix to get a value for a's in terms of u
K= matrix([[1, 0, 0, 0],[0, 1, 0, 0],[1,h,h**2,h**3],[0,1,2*h,3*h**2]])
invK=~K

poly=matrix([1, x, x**2, x**3])

u=poly*invK #this returns the u vector without being multpilied by the a's which is basically the transpose of the phi matrix

phi=transpose(u)
dphi=diff(phi,x)

# A total of 4 degrees of freedom per element, since we are spanning on two elements 
# then we expect a total of 8 different degrees of freedom per element, which could be reduced by specifying the u's
#A[i,j]= integral(phi(i)*phi(j))(ui-1,ui)+integral(phi(i)*phi(j))(ui,ui+1)

#1
Aubub= integral(phi[0]*phi[0],x,0,h)+integral(phi[2]*phi[2],x,0,h)
print "Aui,ui = " + str(Aubub)+'\n'

#2
Adubdub= integral(phi[1]*phi[1],x,0,h)+integral(phi[3]*phi[3],x,0,h)
print "Au'i,u'i = " + str(Adubdub)+'\n'

#3-4
Aubdub= integral(phi[0]*phi[1],x,0,h)+integral(phi[2]*phi[3],x,0,h) #=Adubub
print "Au'i,ui = Aui,u'i =" + str(Aubdub)+'\n'

#5
Auaub= integral(phi[0]*phi[2],x,0,h)
print "Aui-1,ui = " + str(Auaub)+'\n'

#6
Aduadub= integral(phi[1]*phi[3],x,0,h)
print "Au'i-1,u'i = " + str(Aduadub)+'\n'

#7
Aduaub= integral(phi[1]*phi[2],x,0,h) #=Auadub
print "Au'i-1,ui = " + str(Aduaub)+'\n'

#8
Auadub= integral(phi[0]*phi[3],x,0,h) #=Auadub
print "Aui-1,u'i= " + str(Auadub) +'\n'+'\n'


#1
Bubub= integral(dphi[0]*dphi[0],x,0,h)+integral(dphi[2]*dphi[2],x,0,h)
print "Bui,ui = " + str(Bubub)+'\n'

#2
Bdubdub= integral(dphi[1]*dphi[1],x,0,h)+integral(dphi[3]*dphi[3],x,0,h)
print "Bu'i,u'i = " + str(Bdubdub)+'\n'

#3-4
Bubdub= integral(dphi[0]*dphi[1],x,0,h)+integral(dphi[2]*dphi[3],x,0,h) #=Adubub
print "Bu'i,ui = Bui,u'i = " + str(Bubdub)+'\n'

#5
Buaub= integral(dphi[0]*dphi[2],x,0,h)
print "Bui-1,ui = " + str(Buaub)+'\n'

#6
Bduadub= integral(dphi[1]*dphi[3],x,0,h)
print "Bu'i-1,u'i = " + str(Bduadub)+'\n'

#7
Bduaub= integral(dphi[1]*dphi[2],x,0,h) #=Auadub
print "Bu'i-1,ui = " + str(Bduaub)+'\n'

#8
Bduaub= integral(dphi[0]*dphi[3],x,0,h) #=Auadub
print "Bui-1,u'i = " + str(Bduaub)+'\n'

#MS=MatrixSpace(ComplexField(),2,2)
#A=MS(0)
#A[1,1]=5 

# A1= integral(phi(0)*phi(1))  -coressponds to u0-u1
# A2= integral(phi(1)*phi(1))  - 

Aui,ui = 26/35*h

Au'i,u'i = 2/105*h^3

Au'i,ui = Aui,u'i =0

Aui-1,ui = 9/70*h

Au'i-1,u'i = -1/140*h^3

9.2

The Potential Energy of a beam V is given by the following expression: \(V=\int_0^L{{1 \over 2}EI {d^2u \over dx^2}-u(x)f(x)}\,dx\)

where u is approximated as \(u(x)= \sum(a_i,\phi_i)\) where the basis (shape functions) are the hermite interpolation coefficients which have been found in 9.1 above

In [61]: 
E=29000.0 #ksi
I=300.0 #in^4


K= matrix([[1, 0, 0, 0],[0, 1, 0, 0],[1,h,h**2,h**3],[0,1,2*h,3*h**2]])

d2phi=diff(dphi,x) # since A is formed of terms that consist of d^2phi(i)/dx^2
"""
for i in range(4):
    print integral(phi[i],x,0,h)
"""    

M=d2phi*transpose(d2phi)# since A[i,j]= integral(d^2phi(i)/dx^2*d^2phi(j)/dx^2)(ui-1,ui) = integral (M)
a=matrix([[h,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]])
k=copy(a)

for i in range(4):
    for j in range(4):
        k[i,j]=integral(M[i][j],x,0,h)
K=copy(k)     


# Boundary conditions---
for i in range(4):
    for j in range(4):
        if i ==0 or j==0:
            K[i,j]=0
        elif i==1 or j==1:
            K[i,j]=0
K[1,1]=1
K[0,0]=1

# Solving the expression

var('P')        
Force= matrix([[0],[0],[P],[0]])  


u= invK*(K.inverse()*Force)  

u=poly*u  # expression of u(x)

# for a particular case
up=u.subs(P=-100,h=10)
up=(1/E*I)*up

m=[]
uplot=[]

for l in range(0,1000):
    uplot.append(up.subs(x=l/100.0)[0]);
    m.append(l)

plt.plot(m,uplot)









In [56]:





from matplotlib.patches import Rectangle 
import matplotlib.tri as tri

#plt.tripcolor(m, uplot, global_triang)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
Edge= Rectangle([0,0],0.1,2,color="red")
ax.add_patch(Edge)
deltax=0.5
for x in range(len(m)):
    Beam= Rectangle([x/100.0,uplot[x][0]],0.2,0.2,color="red")
    ax.add_patch(Beam)
ax.set_xlim(0,10)
ax.set_ylim(-300,0)
















    Out[56]:






(-300, 0)


































In []: