Решается следующее уравнение на $x \in [0,1]$: $$ -u''(x) + x u(x) = f(x)$$
Граничные условия: $$ u'(0)=u(1)=0$$
Точное решение: $$ u(x) = cos(\frac{5 \pi x}{2})$$
Правая часть: $$ f(x) = (x+\frac{25 \pi^2}{4})cos(\frac{5 \pi x}{2}) $$
$ N^1_{-1}(\xi) = \frac{\xi-1}{2} \text{; }$ $N^2(\xi) = \frac{3(\xi^2 - 1)}{2 \sqrt{6}} \text{; }$ $N^1_{1}(\xi) = \frac{\xi+1}{2}$
Все интегрирования проводились в Wolfram Alpha, а матричные операции - с помощью MatrixCalc.org. Проводим дискретизацию равномерно, чтобы упросить вычисления, имея общий для всех элементов $h$.
$$ K_j = \frac{1}{h} \begin{bmatrix} 1 & 0 & -1\\ 0 & 2& 0\\ -1 & 0 & 1 \end{bmatrix} $$$$ M_j(x_{j-1}, x_j) = \frac{h}{2} \begin{bmatrix} \frac{1}{8} (3 x_{j-1} + x_j) & -\frac{3 x_{j-1} + 2 x_j}{5\sqrt{6}} & \frac{x_{j-1} + x_j}{6}\\ \cdots & \frac{x_{j-1}+x_j}{5} & -\frac{2 x_{j-1}+3 x_j}{5\sqrt{6}}\\ \cdots & \cdots & \frac{1}{6}(x_{j-1}+3 x_j) \end{bmatrix} $$Матрица $M_j$ симметрична, т.е. $M_j^{ab}=M_j^{ba}$поэтому используется символ $\cdots$ для улучшения читаемости.
$$ I_j(f,x_{j-1}, x_j) = \frac{h}{6} \begin{bmatrix} 2 f_{j-1} + f_j \\ -\sqrt{\frac{3}{2}}(f_{j-1} + f_j) & \\ f_{j-1} + 2 f_j \end{bmatrix} $$При расчете $I_j$ использовалась линейная интерполяция вместо расчёта интегралов(Wolfram alpha отказался бесплатно считать такие интегралы).
Комбинируются из $M_j$, $K_j$ и вектора $I_j$, так, что каждый следующий блок перекрывает предыдущий с суммированием в общем элементе. Для учета граничных условий выбрасывается последний столбец и последняя строка для матриц и последний элемент для вектора.
Решаем линейную систему $(K+M)c=I$. В полученном векторе коэфф-в каждые 3 соотв-т своему элементу. Используем разложение по базису $[N^1_{-1}(\xi), N^2(\xi), N^1_1(\xi)]$, транслируя $\forall$ x в соотв-й ему $\xi$ на каноническом элементе через формулу: $$ \xi(x,x_{j-1},x_j) = \frac{x_{j-1}+x_j-2x}{x_{j-1}-x_j} $$
Введем величину погрешности $\Delta(x) = U_{\text{fem}}(x) - U(x)$, где $U_{\text{fem}}(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 [1]:
import numpy as np
from matplotlib import pylab as plt
import math
from matplotlib import rc
from scipy.integrate import simps
%matplotlib inline
a = 0
b = 1
# m-количество точек
m = 20
n_el = m -1
x = np.linspace(a,b,m, endpoint = True)
h = x[1]-x[0]
# Расчет размера глобальных матриц(до учета гр. условий). Расчёт "в лоб"
global_m_size = 3
for i in range(n_el-1):
global_m_size += 2
# Точное решение
def u(x):
return np.cos(5*math.pi*x/2)
# Функция правой части
def f(x):
return (x+25*math.pi*math.pi/4)*np.cos(5*math.pi*x/2)
# Базисные функции
def N_1_m1(xi):
return (xi-1)/2
def N_2(xi):
return 3*(np.square(xi)-1)/(2*np.sqrt(6))
def N_1_1(xi):
return (xi+1)/2
K_j = np.zeros((3,3))
K_j[0][0] = 1
K_j[1][1] = 2
K_j[2][2] = K_j[0][0]
K_j[0][1] = 0
K_j[1][0] = K_j[0][1]
K_j[2][0] = -1
K_j[0][2] = K_j[2][0]
K_j[2][1] = 0
K_j[1][2] = K_j[2][1]
K_j = K_j * 1 /h
def M_j(a, b):
M = np.zeros((3,3))
M[0][0] = (3*a+b)/6
M[1][1] = (a+b)/5
M[2][2] = (3*b+a)/6
M[0][1] = -(3*a+2*b)/(5*np.sqrt(6))
M[1][0] = M[0][1]
M[0][2] = (a+b)/6
M[2][0] = M[0][2]
M[1][2] = -(2*a+3*b)/(5*np.sqrt(6))
M[2][1] = M[1][2]
M = M*h/2
return M
def I_j(a,b):
I_j = np.zeros(3)
I_j[0] = 2*f(a)+f(b)
I_j[1] = -np.sqrt(3/2)*(f(a)+f(b))
I_j[2] = f(a) + 2*f(b)
I_j = I_j * h/6
return I_j
In [2]:
# Решим СЛАУ:
def get_fem_u_can(xi,c):
return c[0]*N_1_m1(xi) + c[1]*N_2(xi) + c[2]*N_1_1(xi)
# преобразования координат на канонический элемент и обратно
def xi_from_x(calc_x, xj_1, xj):
return (xj_1 + xj - 2 * calc_x) / (xj_1 - xj)
def x_from_xi(z, xj_1, xj):
return (1-z)/2*xj_1 + (1+z)/2*xj
def calc_solution(calc_x, x, coeffs):
interval_num = 1
xj_1 = x[0]
xj = x[1]
if calc_x == 0:
return 1
for index, item in enumerate(x):
if calc_x <= item:
interval_num = index
#print(index)
xj_1 = x[index - 1]
xj = x[index]
break
c = [0, 0, 0]
c = coeffs[(interval_num-1)*2:(interval_num-1)*2+3]
xi = xi_from_x(calc_x, xj_1, xj)
res = get_fem_u_can(xi,c)
return res
In [3]:
def solve(m):
n_el = m -1
x = np.linspace(0,1,m, endpoint = True)
h = x[1]-x[0]
print(m)
# Расчет размера глобальных матриц(до учета гр. условий). Расчёт "в лоб"
global_m_size = 3
for i in range(n_el-1):
global_m_size += 2
# Посчитаем глобальные K,M,J
K = np.zeros((global_m_size, global_m_size))
M = np.zeros((global_m_size, global_m_size))
I = np.zeros(global_m_size)
# K
K[0:3, 0:3] = K_j
s = 2
for i in range(n_el-1):
temp = K[s][s]
K[s:s+3,s:s+3] = K_j
K[s][s] += temp
s += 2
# M
M[0:3, 0:3] = M_j(x[0], x[1])
s = 2
for i in range(1,n_el):
temp = M[s][s]
M[s:s+3,s:s+3] = M_j(x[i], x[i+1])
M[s][s] += temp
s += 2
# I
I[0:3] = I_j(x[0], x[1])
s = 2
for i in range(1, n_el):
temp = I[s]
I[s:s+3] = I_j(x[i], x[i + 1])
I[s] += temp
s += 2
# учтем граничные условия
I = I[0:-1]
M = M[0:-1,0:-1]
K = K[0:-1,0:-1]
coeffs = np.concatenate((np.linalg.solve(M+K, I),[0] ))
fem_sol = [calc_solution(item, x, coeffs) for item in x]
return x, fem_sol, h
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 np.array(der_list)
In [4]:
x, sol, h = solve(50)
# строим график
fig = plt.figure(figsize=(15,10))
plt.plot(x, u(x), 'r', label = 'u')
plt.plot(x, sol, 'bo-', label = 'fem')
#plt.title('Результат при N='+str(N), fontsize=20)
plt.legend(fontsize=20)
#plt.xlim(0, Z)
plt.grid(True)
In [5]:
e0_list = []
e1_list = []
N_list = []
for i in range(50, 501, 50):
x, fem_sol, h = solve(i)
x = np.array(x)
fem_sol = np.array(fem_sol)
u_arr = u(x)
delta = fem_sol - u_arr
d_delta = calc_derivatives(delta, h)
e0 = np.sqrt(simps(delta ** 2, x))
e1 = np.sqrt(simps(d_delta ** 2, x))
e0_list.append(e0)
e1_list.append(e1)
N_list.append(i)
In [6]:
fig = plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Errors of FEM solution', fontsize=26)
plt.xlabel('Number of points', fontsize=20)
plt.ylabel('Value of error', fontsize=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.legend()
Out[6]:
In [7]:
def e0(N):
#первая норма
#считаем интеграл от квадрата ошибки методом трапеций
# N - число узлов
n = 2*N - 2 # размерность задачи
# x для построения точного решения
x = np.linspace(0, Z, N)
# nodes для построения нашего решения
h = Z/(N - 1)
nodes = np.zeros(N)
nodes[0] = 0.0
for i in range(1, N):
nodes[i] = nodes[i-1] + h
# получаем коэффициенты в узлах
coeffs = solve(N)
sol = []
for i in range(n):
if (i%2 == 0):
sol.append(coeffs[i])
sol.append(0)
# получаем квадрат ошибки
e = []
for i in range(N):
e.append(u(x[i]) - sol[i])
e2 = [i**2 for i in e]
return (np.trapz(e2, x)**(1/2))
In [8]:
def e1(N):
#вторая норма
#считаем интеграл от квадрата производной ошибки методом трапеций
n = 2*N - 2 # размерность задачи
# x для построения точного решения
x = np.linspace(0, 1, N)
# nodes для построения нашего решения
h = 1/(N - 1)
nodes = np.zeros(N)
nodes[0] = 0.0
for i in range(1, N):
nodes[i] = nodes[i-1] + h
# получаем коэффициенты в узлах
_,sol,_ = solve(N)
e = []
for i in range(N):
e.append(u(x[i]) - sol[i])
# сначала считаем производную
dx = x[1]-x[0] #const
dedx = np.gradient(e, dx)
dedx2 = [i**2 for i in dedx]
#потом считаем интеграл от квадрата
return (np.trapz(dedx2, x)**(1/2))
In [9]:
x = np.array([100, 200, 300, 400, 500])
#y = np.array([e0(100),e0(200),e0(300),e0(400),e0(500)])
y = np.array([e1(100),e1(200),e1(300),e1(400),e1(500)])
z = np.polyfit(np.log10(x), np.log10(y), 1)
# строим график
plt.plot(np.log10(x), np.log10(y), 'go')
plt.plot(np.log10(x), z[0]*np.log10(x) + z[1], '--')
plt.grid(True)
print(z)
In [ ]:
In [ ]: