# Finite Elements (ipynb)¶

## problem 1 (pdf)¶

### 1c. Hermite Polynomials¶

In [22]:
from sympy import *
import numpy as np
from matplotlib import pyplot as plt

In [393]:
x, i, h = symbols('x i h')

expr1 = Piecewise(
(0, (x-i) > h),
(0, (x-i) < -h),
(2*((x-i)/h)**3-3*((x-i)/h)**2+1, (x-i) > 0),
(-2*((x-i)/h)**3-3*((x-i)/h)**2+1, (x-i) <= 0)
)

expr2 = Piecewise(
(0, (x-i) > h),
(0, (x-i) < -h),
(((x-i)/h)**3-2*((x-i)/h)**2+(x-i)/h, (x-i)>0),
(((x-i)/h)**3+2*((x-i)/h)**2+(x-i)/h, (x-i)<=0)
)

L=1
h_val = 1

X = np.linspace(-2,2, 100)

Y1 = [expr1.subs([(i,0),(h,h_val),(x,val)]) for val in X]
Y2 = [expr2.subs([(i,0),(h,h_val),(x,val)]) for val in X]
Y3 = [expr1.subs([(i,1),(h,h_val),(x,val)]) for val in X]
Y4 = [expr2.subs([(i,1),(h,h_val),(x,val)]) for val in X]

plt.plot(X, Y1)
plt.plot(X, Y2)
plt.plot(X, Y3)
plt.plot(X, Y4)

plt.show()


### Matrix Coefficients¶

In [410]:
h_val = 1

exprs = [expr1.subs([(i,0),(h,h_val)]), expr2.subs([(i,0),(h,h_val)]), expr1.subs([(i,1),(h,h_val)]), expr2.subs([(i,1),(h,h_val)])]

rows = np.array([[integrate(a*b, (x,-1,2)) for a in exprs] for b in exprs])

rowsB = np.array([[integrate(diff(a,x)*diff(b,x)-a*diff(b), (x,-1,2)) for a in exprs] for b in exprs])

rows[0,1]=rows[1,0]=rows[2,3]=rows[3,2]=0

print "h = %d" % h_val

print "\n A = \n"

for row in rows:
print row

print "\n B = \n"

for row in rowsB:
print row

h = 1

A =

[26/35 0 9/70 -13/420]
[0 2/105 13/420 -1/140]
[9/70 13/420 26/35 0]
[-13/420 -1/140 0 2/105]

B =

[12/5 1/5 -7/10 0]
[-1/5 4/15 0 -1/20]
[-17/10 -1/5 12/5 1/5]
[1/5 -1/60 -1/5 4/15]


## 2. Beam Bending¶

In [397]:
x, i, h = symbols('x i h')

expr1 = Piecewise(
(0, (x-i) > h),
(0, (x-i) < -h),
(2*((x-i)/h)**3-3*((x-i)/h)**2+1, (x-i) > 0),
(-2*((x-i)/h)**3-3*((x-i)/h)**2+1, (x-i) <= 0)
)

expr2 = Piecewise(
(0, (x-i) > h),
(0, (x-i) < -h),
(((x-i)/h)**3-2*((x-i)/h)**2+(x-i)/h, (x-i)>0),
(((x-i)/h)**3+2*((x-i)/h)**2+(x-i)/h, (x-i)<=0)
)

L=1
h_val = 0.1

exprs2 = [expr1.subs([(i,0),(h, h_val)]), expr2.subs([(i,0),(h, h_val)]), expr1.subs([(i,1),(h, h_val)]), expr2.subs([(i,1),(h, h_val)])]

rows2 = np.array([[integrate(diff(a,x,x)*diff(b,x,x), (x,-1,2)) for a in exprs2] for b in exprs2])

dim = int(2*L/h_val)
A = np.zeros((dim,dim))

A[0,0]=1
A[1,1]=1

for ind in range(2,dim-2,2):
A[ind:ind+4,ind:ind+4] = rows2


In [403]:
def final_exp(f):

B = np.zeros((dim,1))

for ind in range(2,dim-1,2):
B[ind:ind+2,0] = [integrate(a*f,(x,-1,2)) for a in exprs2[:2]]

a_vals = np.linalg.inv(A).dot(B)
bases = []

for ind in np.arange(0,float(L),h_val):
bases.extend([expr1.subs([(i, ind),(h,h_val)]),expr2.subs([(i, ind),(h,h_val)])])

return sum([bases[ind]*a_vals[ind] for ind in range(20)]),bases


In [404]:
Xs = np.arange(0,L,h_val/5)

for f in np.linspace(0,-1,5):
exp,bases = final_exp(f)
Ys = [exp.subs(x,ind) for ind in Xs]
plt.plot(Xs, Ys)
plt.show()

In [ ]:


In [ ]: