Решается следующее уравнение на $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$.
Комбинируются из $M_j$, $K_j$ и вектора $I_j$, так, что каждый следующий блок перекрывает предыдущий с суммированием в общем элементе. Для учета граничных условий выбрасывается последний столбец и последняя строка для матриц и последний элемент для вектора.
In [1]:
import numpy as np
from matplotlib import pylab as plt
import math
%matplotlib inline
a = 0
b = 1
# m-количество точек
m = 50
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
# Посчитаем глобальные 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]
In [2]:
# Решим СЛАУ:
vect = np.linalg.solve(M+K, I)
coeffs = np.concatenate((vect,[0]))
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):
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
def calc_solution2(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
#print(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 = xi_from_x(calc_x, xj_1, xj)
res = c[0]*N_1_m1(cur_qsi) + c[1]*N_2(cur_qsi) + c[2]*N_1_1(cur_qsi)
return res
In [3]:
coeffs[4:4+3]
coeffs
Out[3]:
In [4]:
K
Out[4]:
In [5]:
aaa = np.linspace(x[2],x[3], 60, endpoint=True)
ress = []
resN2 = []
for item in aaa:
xi = xi_from_x(item, x[2], x[3])
#print(xi)
c = coeffs[2*2:2*2+3]
ress.append(get_fem_u_can(xi,c))
resN2.append(N_2(xi)*c[1])
plt.plot(aaa,ress)
plt.plot(aaa,resN2)
Out[5]:
In [6]:
FEM_SOL = []
x_show = np.linspace(0,1, 1000, endpoint=True)
#x_show = x
for item in x_show:
FEM_SOL.append(calc_solution(item))
plt.plot(x_show,FEM_SOL)
Out[6]:
In [7]:
fig = plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('U(x) and FEM solution', fontsize=26)
plt.xlabel('x', fontsize=20)
plt.ylabel('U(x)', fontsize=20)
plt.plot(x, u(x), 'rx--', label = 'U(x)')
#plt.plot(x, f(x), 'bx--', label = 'f(x)')
plt.plot(x_show, FEM_SOL, 'gx', label = 'FEM solution')
#plt.plot(x_show, FEM_SOL - u(x_show), 'bx--', label = 'FEM solution')
plt.legend()
Out[7]:
In [ ]: