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 [ ]: