In [7]:
from sympy import *
from numpy import *
from sympy.abc import x
from matplotlib import pyplot as plt
%matplotlib inline

9.1 (b)

In [9]:
x, h = symbols('x h')
p_0 = x/h;
p_1 = 1-x/h;

print 2*integrate(p_0 * p_0, (x,0,h))
print integrate(p_0 * p_1, (x,0,h))

2*h/3
h/6

9.1 (c)

In [10]:
# x, h = symbols('x h', postive=True, real=True)
h = symbols('h')
M = Matrix([[1,0,0,0], [0,1,0,0], [1,h,h**2,h**3], [0,1,2*h,3*h**2]])
M_i = M.inv()
x_vec = Matrix([[1,x,x**2,x**3]])


p_0 = sum(x_vec * M_i[:,0])
p_1 = sum(x_vec * M_i[:,1])
p_2 = sum(x_vec * M_i[:,2])
p_3 = sum(x_vec * M_i[:,3])

dp_0 = diff(p_0, x)
dp_1 = diff(p_1, x)
dp_2 = diff(p_2, x)
dp_3 = diff(p_3, x)

# A(u_i, u_i)
print "A(u_i,u_i): ", simplify(integrate(p_2*p_2, (x, 0, h)) + integrate(p_0*p_0, (x, 0, h)))
print "A(du_i,du_i): ", simplify(integrate(p_3*p_3, (x, 0, h)) + integrate(p_1*p_1, (x, 0, h)))
print "A(u_i,du_i): ", simplify(integrate(p_2*p_3, (x, 0, h)) + integrate(p_0*p_1, (x, 0, h)))
print "A(u_i-1,u_i): ", simplify(integrate(p_0*p_2, (x, 0, h)))
print "A(du_i-1,du_i): ", simplify(integrate(p_1*p_3, (x, 0, h)))
print "A(du_i-1,u_i): ", simplify(integrate(p_1*p_2, (x, 0, h)))
                               

print "B_(u_i,u_i): ", simplify(integrate(dp_2*dp_2, (x, 0, h)) + integrate(dp_0*dp_0, (x, 0, h)))
print "B_(du_i,du_i): ", simplify(integrate(dp_3*dp_3, (x, 0, h)) + integrate(dp_1*dp_1, (x, 0, h)))
print "B_(u_i,du_i): ", simplify(integrate(dp_2*dp_3, (x, 0, h)) + integrate(dp_0*dp_1, (x, 0, h)))
print "B_(u_i-1,u_i): ", simplify(integrate(dp_0*dp_2, (x, 0, h)))
print "B_(du_i-1,du_i): ", simplify(integrate(dp_1*dp_3, (x, 0, h)))
print "B_(du_i-1,u_i): ", simplify(integrate(dp_1*dp_2, (x, 0, h)))

A(u_i,u_i):  26*h/35
A(du_i,du_i):  2*h**3/105
A(u_i,du_i):  0
A(u_i-1,u_i):  9*h/70
A(du_i-1,du_i):  -h**3/140
A(du_i-1,u_i):  13*h**2/420
B_(u_i,u_i):  12/(5*h)
B_(du_i,du_i):  4*h/15
B_(u_i,du_i):  0
B_(u_i-1,u_i):  -6/(5*h)
B_(du_i-1,du_i):  -h/30
B_(du_i-1,u_i):  -1/10

9.2

In [ ]:
# second derivatives in space
ddp_0 = diff(p_0,x,2)
ddp_1 = diff(p_1,x,2)
ddp_2 = diff(p_2,x,2)
ddp_3 = diff(p_3,x,2)

E = 0.5
I = 0.6
nodes = 10
f = x**2 #what do these normally look like?
n = np.arange(nodes)
n_1 = np.arange(nodes - 1)
A = np.zeros((nodes*2, nodes*2))
#Fill A matrix postive cases
A[n*2,n*2] = 2*integrate(ddp_0*ddp_0, (x,0,h)) 
A[(n*2)+1,(n*2)+1] = 2*integrate(ddp_1*ddp_1, (x,0,h))
A[n*2,(n*2)+1] = integrate(ddp_0*ddp_1, (x,0,h)) + integrate(ddp_2*ddp_3, (x,0,h))
A[(nodes*2)-2,(nodes-1)*2] = integrate(ddp_2*ddp_3, (x,0,h))
# use n_1 so we dont hit edge case
A[(n_1*2)+1,(n_1*2)+2] = integrate(ddp_1*ddp_2, (x,0,h))
A[n_1*2,(n_1*2)+2] = integrate(ddp_0*ddp_2)
A[(n_1*2)+1,(n_1*2)+3] = integrate(ddp_1*ddp_3, (x,0,h))
A[(n_1*2),(n_1*2)+3] = integrate(ddp_0*ddp_3, (x,0,h))
# Fill A matrix negative cases
#A[(n_1*2)-1,(n_1*2)-2] = integrate(ddp_1*ddp_2, (x,0,h))
# we can use this trick instead!
A = A + A.T - diag(A.diagonal)
A *= 0.5*E*I

#Fill B matrix
B = np.zeros(n*2)
# edge values
B[(nodes*2)-2] = integrate(ddp_2*f, (x,0,h))
B[(nodes*2)-1] = integrate(ddp_3*f, (x,0,h))
B[n*2] =   integrate(p_2*f, (x,0,h)) + integrate(p_0*f, (x,0,h))
B[(n*2) + 1] =   integrate(p_3*f, (x,0,h)) + integrate(p_2*f, (x,0,h))

#boundary conditions? 
# not sure how to fix edges?
A[:,0] = np.zeros(nodes*2)
B[0] = 0

a = np.linalg.solve(A,B)
print a
In [ ]:
In [ ]: