In [1]:
import numpy
import matplotlib.pyplot as plt
from scipy import linalg
%matplotlib inline
In [2]:
number_of_points = 1000
number_of_elements = number_of_points - 1
In [3]:
x = numpy.linspace(0, 1, number_of_points)
In [4]:
h = 1 / (number_of_points - 1)
In [5]:
def f(x):
return x * numpy.sin(3 * numpy.pi * x / 2) + (9/4) * (numpy.pi ** 2) * numpy.sin(3 * numpy.pi * x /2)
In [6]:
def expected_solution_calc(x):
return numpy.sin(3 * numpy.pi * x / 2)
In [7]:
K = numpy.zeros((3, 3))
K[0][0] = K[2][2]= 19 / 6
K[1][0] = K[0][1] = - 16 / 3
K[2][0] = K[0][2] = 13 / 6
K[1][1] = 32 / 3
K[1][2] = K[2][1] = - 16 / 3
K *= 2 / h
In [8]:
K
Out[8]:
In [9]:
def calc_M(a, b):
M = numpy.zeros((3, 3))
M[0][0] = (7 * a + b) / 30
M[1][0] = M[0][1] = 2 * a / 15
M[2][0] = M[0][2] = (-a - b) / 30
M[1][1] = 8 * (a + b) / 15
M[1][2] = M[2][1] = 2 * b / 15
M[2][2] = (a + 7 * b) / 30
M *= h / 2
return M
In [10]:
def calc_I(x_j1, x_j):
f_j1 = f(x_j1)
f_j = f(x_j)
I = numpy.zeros(3)
I[0] = 1/3 * f_j1
I[1] = 2/3 * (f_j1 + f_j)
I[2] = 1/3 * f_j
I *= h/2
return I
In [11]:
calc_I(0,0.1).shape
Out[11]:
In [12]:
K
Out[12]:
In [13]:
global_size = 3
for _ in range(number_of_elements - 1):
global_size += 2
In [14]:
K_global = numpy.zeros((global_size, global_size))
In [15]:
K_global[0:3, 0:3] = K
shift = 2
for _ in range(number_of_elements - 1):
temp = K_global[shift, shift]
K_global[shift:shift + 3, shift:shift + 3] = K
K_global[shift, shift] += temp
shift += 2
#K_global[global_size - 1, global_size - 1] += K[0][0]
In [16]:
K_global
Out[16]:
In [17]:
M_global = numpy.zeros((global_size, global_size))
In [18]:
M_global[0:3, 0:3] = calc_M(x[0], x[1])
shift = 2
for i in range(1, number_of_elements):
temp = M_global[shift, shift]
M_global[shift:shift + 3, shift:shift + 3] = calc_M(x[i], x[i + 1])
M_global[shift, shift] += temp
shift += 2
In [19]:
I_global = numpy.zeros(global_size)
In [20]:
I_global[0:3] = calc_I(x[0], x[1])
shift = 2
for i in range(1, number_of_elements):
temp = I_global[shift]
I_global[shift:shift+3] = calc_I(x[i], x[i + 1])
I_global[shift] += temp
shift += 2
In [21]:
K_global_bc = K_global[1:,1:]
In [22]:
M_global_bc = M_global[1:, 1:]
In [23]:
I_global_bc = I_global[1:]
In [24]:
I_global_bc.shape
Out[24]:
In [25]:
matr = K_global_bc + M_global_bc
In [26]:
matr.shape
Out[26]:
In [27]:
solution = linalg.solve(matr, I_global_bc)
In [28]:
solution.shape
Out[28]:
In [29]:
coeffs = numpy.concatenate(([0], solution))
#coeffs = solution
In [30]:
coeffs.shape
Out[30]:
In [31]:
def mapping(calc_x, xj_1, xj):
return (xj_1 + xj - 2 * calc_x) / (xj_1 - xj)
In [32]:
def mapping2(z, xj_1, xj):
return (1-z)/2*xj_1 + (1+z)/2*xj
In [33]:
def N0(qsi):
return 1/2 * ((qsi ** 2) - qsi)
In [34]:
def N1(qsi):
return 1 - qsi ** 2
In [35]:
def N2(qsi):
return 1/2 * ((qsi ** 2) + qsi)
In [36]:
def calc_solution(calc_x):
interval_num = 0
xj_1 = 0
xj = 0
if (calc_x < 0) and (calc_x > 1):
return False
if calc_x == 0:
return 0
for index, item in enumerate(x):
if calc_x <= item:
interval_num = index
xj_1 = x[index - 1]
xj = x[index]
break
c = [0, 0, 0]
for i in range(3):
c[i] = coeffs[(interval_num - 1) * 2 + i]
cur_qsi = mapping(calc_x, xj_1, xj)
res = c[0] * N0(cur_qsi) + c[1] * N1(cur_qsi) + c[2] * N2(cur_qsi)
return res
In [45]:
aaa = numpy.linspace(x[2],x[3], 60, endpoint=True)
ress = []
resN2 = []
for item in aaa:
cur_qsi = mapping(item, x[2], x[3])
print(cur_qsi)
c = coeffs[2*2:2*2+3]
ress.append(c[0] * N0(cur_qsi) + c[1] * N1(cur_qsi) + c[2] * N2(cur_qsi))
resN2.append(N2(cur_qsi)*c[1])
plt.plot(aaa,ress)
plt.plot(resN2)
Out[45]:
In [40]:
eq_solution = []
new_x = numpy.linspace(0, 1, 10000, endpoint=True)
x_for_ex_solution = numpy.linspace(0, 1, 10)
for item in x:
#print(item)
temp = calc_solution(item)
eq_solution.append(temp)
#print(temp)
In [38]:
expected_solution = expected_solution_calc(x_for_ex_solution)
In [41]:
plt.plot(x, eq_solution, 'o')
plt.plot(x_for_ex_solution, expected_solution, 'x')
Out[41]: