Метод канторовича для решения эллиптической краевой задачи

@Dzxtckfd · Регистрация: 17.11.2021

Студворк — интернет-сервис помощи студентам

Всем привет, нужна помощь в реализации решения эллиптической краевой задачи методом Канторовича.
Все формулы и вычисления взяты из книги методы конечных элементов для уравнений с частными производными Автор: Митчел,Уэйт( т е в формулах предоставленных мной нет ошибки они взяты из книги)
Задача заключается в следующем:
Дано эллиптическое уравнения с граничными условиями, если рассказывать коротко, то задача приводиться к нахождению решения в виде
U(x,y)=∑_(i=1)^n α_i (x) φ_i (y)
n задается самому
φ_i (y)= матрица базисных функций размером n*n ( вычисленная в коде приведенном ниже)
шаг сетки равен pi/n
α_i (x) находиться из системы диф уравнений следующего вида с граничными условиями в точках:

∑_(j=1)^n (ai*cij- ai''*dij) =bi α_i (-π/2)=α_i (π/2)=0 i=1,2,3……n

коэффициенты cij dij и вектор значений bi вычисляются по приведенным в коде(взято так же из книги )

в моем случае cij и dij матрицы размерностью n*n, bi- вектор размерностью (1,n)
матрица cij dij являются три диагональными.
У матрицы cij на диагонали элементы имеют вид ( 1/h, 2/h,2/h......1/h), над диагональные и под диагональные элементы имеют вид -1/h
У матрицы dij на диагонали элементы имеют вид (h/3,2h/3,2h/3.....h/3) над диагональные и под диагональные элементы имеют вид -h/6 ( взято из теории кусочно-линейных функций из книги под авторством Сагдеева Копысова , введение в методы конечных элементов с14 )

Моя проблема заключается в нахождении a_i(x), он должен получиться в виде матрицы размерностью (1,n) где каждая строка матрица это вектор a_i(xi) для каждого x в области [-pi/2, pi/2] с шагом h=pi/N_x( в данном случае возможно братm количество разбиений не N а например N_x другому числу, но для начала хотелось бы решить с N=N_x).
То есть в конечном итоге при N=5 ( к примеру) должна получиться матрица (1,5) где каждая строка это вектор из 5 значений.
С учетом граничных условий вид данной матрицы будет таков что первый и последний вектор это 0, а в каждом другом векторе первое и последнее значение тоже 0( получается своего рода матрица где граница из 0 по периметру. Так же она должна получиться симметричной относительно центра матрицы.
Я пробовал решать двумя способами, первый это используя матричную прогонку находя решение в виде ai(x) = Ai+1* ai+1 - Bi+1
A1=0 ( нулевая матрица)
где Ai= (dij/h^2+ cij- (dij/h^2*Ai-1))^-1 * dij/h^2
B1=0
Bi=(dij/h^2+ cij- (dij/h^2*Ai-1))^-1 * (bi +(dij/h^2)* Bi-1)
Взято из книги Самарского " Теория разностных схем", параграф метод матричной прогонки

Так же пробовал решить через систему диф уравнений вида :
a_i '(x) = q
a_i ''(x) = dij ^-1 * cij * q - dij * bi
пробовал решать через scipy.integrate.solve_bvp
Во всех моих попытках решения я не получил желаемого результата в виде симметричной относительно центра матрицы, мои графики точного и приближенного решения не сходятся по одной из сторон

Точное решение данного уравнения это 3d график функции на участке по x и y [-pi/2, pi/2] ( приведен на фото )

Мой код вычислений коэффициентов cij dij и вектора bi ниже :

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
import numpy as np
from scipy.integrate import quad
 
def phi_i(y, i, h):
    if i == 0:
        if -np.pi/2 <= y <= -np.pi/2 + h:
            return (-np.pi/2 + h - y) / h
        else:
            return 0
    elif i == N:
        if np.pi/2 - h <= y <= np.pi/2:
            return (y - (np.pi/2 - h)) / h
        else:
            return 0
    else:
        if -np.pi/2 + (i-1)*h <= y <= -np.pi/2 + i*h:
            return (y - (-np.pi/2 + (i-1)*h)) / h
        elif -np.pi/2 + i*h <= y <= -np.pi/2 + (i+1)*h:
            return ((-np.pi/2 + (i+1)*h) - y) / h
        else:
            return 0
 
def d_phi_dy(y, i, h):
    if i == 0:
        if -np.pi/2 <= y <= -np.pi/2 + h:
            return -1 / h
        else:
            return 0
    elif i == N:
        if np.pi/2 - h <= y <= np.pi/2:
            return 1 / h
        else:
            return 0
    else:
        if -np.pi/2 + (i-1)*h <= y <= -np.pi/2 + i*h:
            return 1 / h
        elif -np.pi/2 + i*h <= y <= -np.pi/2 + (i+1)*h:
            return -1 / h
        else:
            return 0
 
def calculate_c_ij(i, j, h):
    integrand = lambda y: d_phi_dy(y, i, h) * d_phi_dy(y, j, h)
    result, _ = quad(integrand, -np.pi/2, np.pi/2, limit=1000)
    if j == N - 1 and i == N - 1 :
        result /= 2
    return result
 
def calculate_d_ij(i, j, h):
    integrand = lambda y: phi_i(y, i, h) * phi_i(y, j, h)
    result, _ = quad(integrand, -np.pi/2, np.pi/2, limit=1000)
    if j == N - 1 and i == N - 1 :
        result /= 2
    return result
 
def calculate_b_i(i, h):
    integrand = lambda y: 2 * phi_i(y, i, h)
    result, _ = quad(integrand, -np.pi/2, np.pi/2, limit=1000)
    if i == N - 1:
        result /= 2
    return -result
 
 
N = 5
h = np.pi / N
y_values = np.linspace(-np.pi/2, np.pi/2, N + 1)
c_matrix = np.zeros((N, N))
d_matrix = np.zeros((N, N))
b_vector = np.zeros(N)
 
 
for i in range(N):
    for j in range(N):
        c_matrix[i, j] = calculate_c_ij(i, j, h)
        d_matrix[i, j] = calculate_d_ij(i, j, h)
    b_vector[i] = calculate_b_i(i, h)
 
 
print("Matrix C:")
print(c_matrix)
print("Matrix Alpha:")
print(d_matrix)
print("Vector b:")
print(b_vector)

Если кто то поможет решить данную задачу или подскажет какие шаги решения можно предпринять для поиска необходимой матрицы ai(x) буду очень признателен
Так же предоставляю функцию вычисления ai(x) методом матрчиной прогонки который работает но не так как нужно :

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
N_x=5
a_list=[]
b_list=[]
alpha_list=[]
alpha_1=[0] * N
a = np.zeros((N, N))
b = np.zeros((1,N ))
 
a_list.append(a)
b_list.append(b)
 
lambda_ = 0.01
c_matrix += lambda_ * np.eye(N)
 
a_prev = np.zeros((N, N))
b_prev = np.zeros(N)
 
for i in range(1, N_x):
    h_i = np.pi / N_x
    a_i = np.linalg.inv(((d_matrix * 2 )/ (h ** 2)) + c_matrix - ((d_matrix / (h ** 2)).dot(a_prev))).dot((d_matrix / (h ** 2)))
    b_i = np.linalg.inv(((d_matrix * 2 )/ (h ** 2)) + c_matrix - ((d_matrix / (h ** 2)).dot(a_prev))).dot((b_vector + ((d_matrix / (h ** 2)) @ b_prev)))
    a_list.append(a_i)
    b_list.append(b_i)
    a_prev = a_i 
    b_prev=b_i
 
def compute_alpha(a_list, b_list):
 
    alpha_list = []
    alpha_next = np.zeros(N_x)
 
    for i in range(N_x ,0, -1):
        if i == N :
            alpha = np.zeros(N_x)
        else:
            alpha = alpha_next @  a_list[i-1] + b_list[i-1]
            alpha_next = alpha.T
 
        alpha = np.zeros(N_x) if i == 1 else alpha 
 
        alpha_list.insert(0, alpha)
        
        
 
    for i, alpha in enumerate(alpha_list):
        alpha_list[i][0] = 0
        alpha_list[i][-1] = 0
 
    return alpha_list
 
alpha_values = compute_alpha(a_list, b_list)
print('alpha_i(x)',np.array(alpha_values))

Новые блоги и статьи Все статьи Все блоги /
Подстановка значения реквизита справочника в табличную часть документа Maks 10.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "ПланированиеПерсонала", разработанного в конфигурации КА2. Задача: при выборе сотрудника (справочник Сотрудники) в ТЧ документа. . .	Очистка реквизитов документа при копировании Maks 09.04.2026 Алгоритм из решения ниже применим как для типовых, так и для нетиповых документов на самых различных конфигурациях. Задача: при копировании документа очищать определенные реквизиты и табличную. . .	модель ЗдравоСохранения 8. Подготовка к разному выполнению заданий anaschu 08.04.2026 https:/ / github. com/ shumilovas/ med2. git main ветка * содержимое блока дэлэй из старой модели теперь внутри зайца новой модели 8ATzM_2aurI	Блокировка документа от изменений, если он открыт у другого пользователя Maks 08.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа, разработанного в конфигурации КА2. Задача: запретить редактирование документа, если он открыт у другого пользователя. / / . . .
Система безопасности+живучести для сервера-слоя интернета (сети). Двойная привязка. Hrethgir 08.04.2026 Далее были размышления о системе безопасности. Сообщения с наклонным текстом - мои. А как нам будет можно проверить, что ссылка наша, а не подделана хулиганами, которая выбросит на другую ветку и. . .	Модель ЗдрввоСохранения 7: больше работников, больше ресурсов. anaschu 08.04.2026 работников и заданий может быть сколько угодно, но настроено всё так, что используется пока что только 20% kYBz3eJf3jQ	Дальние перспективы сервера - слоя сети с космологическим дизайном интефейса карты и логики. Hrethgir 07.04.2026 Дальнейшее ближайшее планирование вывело к размышлениям над дальними перспективами. И вот тут может быть даже будут нужны оценки специалистов, так как в дальних перспективах всё может очень сильно. . .	Горе от ума kumehtar 07.04.2026 Эта мне ментальная установка, что вот прямо сейчас, мол, мне для полного счастья не хватает (нужное вписать), и когда я этого достигну - тогда и полный кайф. Одна из самых сильных ловушек на пути. . . .