<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Цветков Андрей, 405 группа
"Одномерный метод конечных элементов"
Постановка задачи
Step1:  Матрица $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}$
Step2:  Аналогично
Step3:  Вектор нагрузки также можно записать в виде скалярного произведения
Step4:  Глобальные матрицы конструируются комбинированием элементарных. Они перекрываются по элементам с целыми индексами, а именно последний элемент последней строки матрицы прибавляется к первому элементу первой строки матрицы. В случае с вектором нагрузки последний элемент предыдущего вектора прибаляется к первому элементу следующего.
В силу граничных условий первый столбец и первую строку всех глобальный матриц, а так же первый элемент вектора нагрузок необходимо убрать.
Step5:  В итоге, задача может быть переписана в виде
Step6:  Проведем анализ погрешности. Введем величину погрешности $\Delta(x) = u(x) - u_{ex}(x)$, где $u_{ex}(x)$ - точное решение уравнения.
Найдем две нормы погрешности
Step7:  Построим зависимость обеих норм погрешности от числа точек
<ASSISTANT_TASK:>
Python Code:
import numpy
import matplotlib.pyplot as plt
from matplotlib import rc
from scipy import linalg
from scipy.integrate import simps
%matplotlib inline
def f(x):
    return x * numpy.sin(3 * numpy.pi * x / 2) + (9/4) * (numpy.pi ** 2) * numpy.sin(3 * numpy.pi * x /2)
def expected_solution_calc(x):
    return numpy.sin(3 * numpy.pi * x / 2)
Explanation: Цветков Андрей, 405 группа
"Одномерный метод конечных элементов"
Постановка задачи:
Используя одномерный метод конечных элементов, решить уравнение $-u''(x) + xu(x) = f(x)$ на интервале [0, 1]. $u(0)=u'(1)=0$. Функция $f(x)$ определяется из условия, что $u(x) = sin(\frac{3 \pi x}{2})$ является точным решением.  
Описание метода:
При решении были использованны лагранжевы КЭ второго порядка. Лагранжев базис на каноническом элементе представляет из себя набор из трех функций:
\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}
End of explanation
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
Explanation: Матрица $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}$
End of explanation
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
Explanation: Аналогично:
\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}$
End of explanation
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
Explanation: Вектор нагрузки также можно записать в виде скалярного произведения:
$(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}
End of explanation
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
def mapping(calc_x, xj_1, xj):
    return (xj_1 + xj - 2 * calc_x) / (xj_1 - xj)
def mapping2(z, xj_1, xj):
    return (1-z)/2*xj_1 + (1+z)/2*x
def N0(qsi):
    return 1/2 * ((qsi ** 2) - qsi)
def N1(qsi):
    return 1 - qsi ** 2
def N2(qsi):
    return 1/2 * ((qsi ** 2) + qsi)
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
Explanation: Глобальные матрицы конструируются комбинированием элементарных. Они перекрываются по элементам с целыми индексами, а именно последний элемент последней строки матрицы прибавляется к первому элементу первой строки матрицы. В случае с вектором нагрузки последний элемент предыдущего вектора прибаляется к первому элементу следующего.
В силу граничных условий первый столбец и первую строку всех глобальный матриц, а так же первый элемент вектора нагрузок необходимо убрать.
End of explanation
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
x, eq_solution, h = numpy.array(solve(100))
expected_solution = expected_solution_calc(x)
plt.plot(x, expected_solution, 'x', label='Точное решение')
plt.plot(x, eq_solution, label='Полученное решение')
plt.title('N=100')
plt.legend()
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)
Explanation: В итоге, задача может быть переписана в виде:
\begin{equation}
    \mathbf{c(K+M)d^T = Id^T}
\end{equation}
Так как уравнение должно выполняться для любого $\mathbf{d}$, то достаточно решить линейную задачу:
\begin{equation}
    \mathbf{(K+M)с^T = I^T}
\end{equation}
End of explanation
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)
Explanation: Проведем анализ погрешности. Введем величину погрешности $\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}
End of explanation
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()
Explanation: Построим зависимость обеих норм погрешности от числа точек:
End of explanation
<END_TASK>