Цветков Андрей, 405 группа
При решении были использованны лагранжевы КЭ второго порядка. Лагранжев базис на каноническом элементе представляет из себя набор из трех функций: \begin{equation*} N_0 = \frac{\xi(\xi - 1)}{2}, \ N_1 = 1 - \xi^2, \ N_2 = \frac{\xi(\xi + 1)}{2} \end{equation*}
Решение уравнения на каноническом элементе представляет из себя линейную комбинацию базисных функций.
Пересчет с канонического элемента на физический осуществляется при помощи следующего соотношения: \begin{equation*} x(\xi) = \frac{1 - \xi}{2} x_{j-1} + \frac{1 + \xi}{2} x_j \end{equation*}
Следующим этапом в решении будет переход от краевой задачи к вариационной. В такой терминологии задача формулируется так: необходимо найти такую функцию $u(x)$, что для любой $v(x)$ будет выполнено следующее равенство: \begin{equation*} \int \limits_0^1 u'(x)v'(x) + x u(x) v(x) dx = \int \limits_0^1 f(x)v(x), \ \forall v(x) \end{equation*}
Заменяя интеграл по всему интервалу на сумму интегралов по интервалам $[x_{j-1}, x_j]$ и подставляя вместо функций $u(x)$ и $v(x)$ элементные функции, пересчитанные с канонического элемента, получим следующее выражение:
\begin{equation*} \sum \limits_{j=1}^N A^S_j(U,V) + A^M_j(U,V) - (f, V)_j = 0 \end{equation*}Были использованы следующие обозначения:
Внутренняя энергия $A^S_j = \frac{h}{2} \int \limits_{-1}^1 \frac{dU(\xi)}{d\xi} \frac{dV(\xi)}{d\xi} d\xi$
Внешняя энергия $A^M_j = \int \limits_{-1}^1 x(\xi) U(\xi) V(\xi) d\xi$
Вектор наргузки $(f, V)_j = \int \limits_{-1}^1 f(x\left(\xi)\right) V(\xi) d\xi$
Функции $U(\xi)$ и $V(\xi)$ можно представить в виде скалярных произведений: \begin{equation*} U(\xi) = \begin{bmatrix} c_{j-1} & c_{j - \frac{1}{2}} & c_j \end{bmatrix} \begin{bmatrix} N_0 \\ N_1 \\ N_2 \end{bmatrix} \end{equation*} \begin{equation*} V(\xi) = \begin{bmatrix} d_{j-1} & d_{j - \frac{1}{2}} & d_j \end{bmatrix} \begin{bmatrix} N_0 \\ N_1 \\ N_2 \end{bmatrix} \end{equation*}
Тогда представление для внутренней энергии запишется в виде: \begin{equation*} A^S_j = \begin{bmatrix}c_{j-1} & c_{j - \frac{1}{2}} & c_j \end{bmatrix} K_j \begin{bmatrix} d_{j-1} \\ d_{j - \frac{1}{2}}\\ d_j \end{bmatrix} \end{equation*}
In [1]:
import numpy
import matplotlib.pyplot as plt
from matplotlib import rc
from scipy import linalg
from scipy.integrate import simps
%matplotlib inline
In [2]:
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 [3]:
def expected_solution_calc(x):
return numpy.sin(3 * numpy.pi * x / 2)
Матрица $K_j = \frac{2}{h} \int \limits_{-1}^{1} \frac{d}{d\xi} \begin{bmatrix} N_0 \\ N_1 \\ N_2 \end{bmatrix} \frac{d}{d\xi} \begin{bmatrix} N_0 & N_1 & N_2 \end{bmatrix} d\xi = \frac{2}{h} \begin{bmatrix} \frac{19}{6} & -\frac{16}{3} & \frac{13}{6} \\ -\frac{16}{3} & \frac{32}{3} & -\frac{16}{3} \\ \frac{13}{6} & -\frac{16}{3} & \frac{19}{6} \end{bmatrix}$
In [4]:
def calc_K(h):
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
return K
Аналогично: \begin{equation*} A^M_j = \begin{bmatrix}c_{j-1} & c_{j - \frac{1}{2}} & c_j \end{bmatrix} M_j \begin{bmatrix} d_{j-1} \\ d_{j - \frac{1}{2}}\\ d_j \end{bmatrix} \end{equation*}
$M_j = \frac{h}{2} \int \limits_{-1}^{1} x(\xi) \begin{bmatrix} N_0 \\ N_1 \\ N_2 \end{bmatrix} \begin{bmatrix}N_0 & N_1 & N_2 \end{bmatrix} d\xi = \frac{h}{2} \begin{bmatrix} \frac{7x_{j-1} + x_j}{30} & \frac{2x_{j-1}}{15} & -\frac{x_{j-1} + x_j}{30} \\ \frac{2x_{j-1}}{15} & \frac{8(x_{j-1} + x_j)}{15} & \frac{2x_j}{15} \\ -\frac{x_{j-1} + x_j}{30} & \frac{2x_j}{15} & \frac{7x_j + x_{j - 1}}{30} \end{bmatrix}$
In [5]:
def calc_M(xj_1, xj, h):
M = numpy.zeros((3, 3))
M[0][0] = (7 * xj_1 + xj) / 30
M[1][0] = M[0][1] = 2 * xj_1 / 15
M[2][0] = M[0][2] = (-xj_1 - xj) / 30
M[1][1] = 8 * (xj_1 + xj) / 15
M[1][2] = M[2][1] = 2 * xj / 15
M[2][2] = (xj_1 + 7 * xj) / 30
M *= h / 2
return M
Вектор нагрузки также можно записать в виде скалярного произведения:
$(f, V)_j = I_j \begin{bmatrix}d_{j-1} & d_{j - \frac{1}{2}} & d_j \end{bmatrix}^T$
где $I_j = \frac{h}{2} \int \limits_{-1}^1 f(x(\xi)) \begin{bmatrix}N_0 & N_1 & N_2 \end{bmatrix} d\xi$
Для вычисления интеграла воспользуемся линейной апроксимацией функции $f(x)$: \begin{equation*} f(x) \approx \begin{bmatrix}f_{j-1} & f_j \end{bmatrix} \begin{bmatrix} N_0 \\ N_2 \end{bmatrix} \end{equation*} где $f_j = f(x_j)$
Тогда: \begin{equation*} I_j = \frac{h}{2} \int \limits_{-1}^1 \begin{bmatrix}f_{j-1} & f_j \end{bmatrix} \begin{bmatrix} N_0 \\ N_2 \end{bmatrix} \begin{bmatrix}N_0 & N_1 & N_2 \end{bmatrix} d\xi = \frac{h}{2} \begin{bmatrix}\frac{1}{3}f_{j-1} \\ \frac{2}{3} \left( f_{j-1} + f_j \right) \\ \frac{1}{3} f_j \end{bmatrix}^T \end{equation*}
In [6]:
def calc_I(xj_1, xj, h):
fj_1 = f(xj_1)
fj = f(xj)
I = numpy.zeros(3)
I[0] = 1/3 * fj_1
I[1] = 2/3 * (fj_1 + fj)
I[2] = 1/3 * fj
I *= h/2
return I
Глобальные матрицы конструируются комбинированием элементарных. Они перекрываются по элементам с целыми индексами, а именно последний элемент последней строки матрицы прибавляется к первому элементу первой строки матрицы. В случае с вектором нагрузки последний элемент предыдущего вектора прибаляется к первому элементу следующего.
В силу граничных условий первый столбец и первую строку всех глобальный матриц, а так же первый элемент вектора нагрузок необходимо убрать.
In [7]:
def create_global_matrixes(number_of_points):
x = numpy.linspace(0, 1, number_of_points)
number_of_elements = number_of_points - 1
h = 1 / number_of_elements
global_size = 3
for _ in range(number_of_elements - 1):
global_size += 2
K_global = numpy.zeros((global_size, global_size))
M_global = numpy.zeros((global_size, global_size))
I_global = numpy.zeros(global_size)
K = calc_K(h)
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
M_global[0:3, 0:3] = calc_M(x[0], x[1], h)
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], h)
M_global[shift, shift] += temp
shift += 2
I_global[0:3] = calc_I(x[0], x[1], h)
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], h)
I_global[shift] += temp
shift += 2
return x, K_global, M_global, I_global
In [8]:
def mapping(calc_x, xj_1, xj):
return (xj_1 + xj - 2 * calc_x) / (xj_1 - xj)
In [9]:
def mapping2(z, xj_1, xj):
return (1-z)/2*xj_1 + (1+z)/2*x
In [10]:
def N0(qsi):
return 1/2 * ((qsi ** 2) - qsi)
In [11]:
def N1(qsi):
return 1 - qsi ** 2
In [12]:
def N2(qsi):
return 1/2 * ((qsi ** 2) + qsi)
In [13]:
def calc_solution(calc_x, x, coeffs):
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
В итоге, задача может быть переписана в виде: \begin{equation*} \mathbf{c(K+M)d^T = Id^T} \end{equation*}
Так как уравнение должно выполняться для любого $\mathbf{d}$, то достаточно решить линейную задачу: \begin{equation*} \mathbf{(K+M)с^T = I^T} \end{equation*}
In [14]:
def solve(number_of_points):
x, K_global, M_global, I_global = create_global_matrixes(number_of_points)
h = 1 / (number_of_points - 1)
K_global_bc = K_global[1:,1:]
M_global_bc = M_global[1:, 1:]
I_global_bc = I_global[1:]
matr = K_global_bc + M_global_bc
solution = linalg.solve(matr, I_global_bc)
coeffs = numpy.concatenate(([0], solution))
eq_solution = []
for item in x:
temp = calc_solution(item, x, coeffs)
eq_solution.append(temp)
return x, eq_solution, h
In [15]:
x, eq_solution, h = numpy.array(solve(100))
In [16]:
expected_solution = expected_solution_calc(x)
In [17]:
plt.plot(x, expected_solution, 'x', label='Точное решение')
plt.plot(x, eq_solution, label='Полученное решение')
plt.title('N=100')
plt.legend()
Out[17]:
In [18]:
def calc_derivatives(delta, h):
der_list = [0]
l = len(delta)
for i in range(1, l - 1):
der_list.append((delta[i + 1] - delta[i - 1]) / (2 * h))
der_list.append(0)
return numpy.array(der_list)
Проведем анализ погрешности. Введем величину погрешности $\Delta(x) = u(x) - u_{ex}(x)$, где $u_{ex}(x)$ - точное решение уравнения.
Найдем две нормы погрешности: \begin{equation*} ||\Delta||_0 = \left( \int |\Delta(x)|^2 dx \right)^{\frac{1}{2}} \end{equation*}
\begin{equation*} ||\Delta||_1 = \left( \int \left| \frac{d\Delta(x)}{dx} \right|^2 dx \right)^{\frac{1}{2}} \end{equation*}
In [19]:
e0_list = []
e1_list = []
N_list = []
for i in range(50, 501, 50):
x, eq_solution, h = numpy.array(solve(i))
expected_solution = expected_solution_calc(x)
delta = eq_solution - expected_solution
d_delta = calc_derivatives(delta, h)
e0 = numpy.sqrt(simps(delta ** 2, x))
e1 = numpy.sqrt(simps(d_delta ** 2, x))
e0_list.append(e0)
e1_list.append(e1)
N_list.append(i)
Построим зависимость обеих норм погрешности от числа точек:
In [20]:
rc('text', usetex=True)
plt.plot(N_list, e0_list, 'x', label='$||\Delta||_0$')
plt.plot(N_list, e1_list, 'x', label='$||\Delta||_1$')
plt.xlabel('N')
plt.legend()
Out[20]:
Как можно видеть из графика, оби нормы погрешности уменьшаются с ростом числа точек.