пример работы. Задача (2), (3) называется задачей Дирихле для уравения Пуассона (Лапласса)
 Скачать 2.22 Mb.
|
|
Содержание 1 Постановка задачи 3 2 Построение разностной аппроксимации и тестовых примеров 3 2.1 Построение разностной аппроксимации . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Построение первого тестового примера . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Построение второго тестового примера . . . . . . . . . . . . . . . . . . . . . . . 5 3 Итерационные методы решения задачи Дирихле 6 3.1 Решение первого тестового примера . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 Решение второго тестового примера . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.3 Решение исходной задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4 Проекционные методы решения задачи Дирихле 13 5 Заключение 14 6 Список литературы 15 7 Приложение 16 2 1 Постановка задачи Прямоугольная металлическая пластина с вырезом используется как теплоотводящий эле- мент. В угловом вырезе пластины (границы ? 2 и ? 3 ) расположен источник тепла. Распреде- ление температуры T (x, y) по площади пластины описывается уравнением Лапласа: ? 2 T ?x 2 + ? 2 T ?y 2 = 0 (1) Рис. 1: Конфигурация пластины Требуется найти распределение T (x, y). 2 Построение разностной аппроксимации и тестовых при- меров 2.1 Построение разностной аппроксимации Имеем задачу вида: ?u(x, y) = ?f (x, y) (2) u| ? = ?(x, y) (3) Здесь ?u = u xx + u yy ? оператор Лапласа, уравнение (2) ? уравнение Пуассона отно- сительно неизвестной функции u(x, y), функция f(x, y) известна, (3) ? граничные условия первого рода. В нашем случае f(x, y) ? 0 и уравнение (2) называется уравением Лапласа. Задача (2) , (3) называется задачей Дирихле для уравения Пуассона (Лапласса). Для численного решения поставленной задачи Дирихле воспользуемся методом конеч- ных разностей. Для аппроксимации уравнения (2) возьмјм пятиточечный шаблон. 3 Рис. 2: Двумерная пространственная сетка и пятиточечный шаблон По каждой из пространственных переменных введјм операторы второй разностной про- изводной ? 1 [U ](x i , y j ) = U i?1,j ? 2U i,j + U i+1,j h 2 1 для i = 1, . . . , N 1 ? 1 и j = 0, . . . , N 2 и ? 2 [U ](x i , y j ) = U i,j?1 ? 2U i,j + U i,j+1 h 2 2 для i = 0, . . . , N 1 и j = 1, . . . , N 2 ? 1 В каждой внутренней точки пространственной сетки построим разностное уравнение ? 1 [U ] + ? 2 [U ] = ?F. (4) F i,j = f (x i , y i ) ? аппроксимация правой части уравнения (2) в случае непрерывности функции f. В развјрнутом виде уравнение (4) имеет вид U i?1,j ? 2U i,j + U i+1,j h 2 1 + U i,j?1 ? 2U i,j + U i,j+1 h 2 2 = ?F i,j Поскольку граничные условия первого рода, то они не нуждаются в аппроксимации. В случае прямоугольной области граничные условия можно записать в виде U i,0 = ?(x i ), U i,N 2 = ? 2 (x i ), i = 0, 1, . . . , N 1 , U 0,j = ? 1 (y j ), U N 1 ,j = ? 2 (y j ), j = 0, 1, . . . , N 2 Построенная дискретная задача называется пятиточечной разностной схемой для за- дачи Дирихле. 4 2.2 Построение первого тестового примера Необходимо строить тестовые примеры, для того чтобы можно было оценить корректность реализации методов решения задачи. Тестовые примеры будем строить по известному точ- ному решению. То есть зафиксируем некоторую функцию u(x, y) и границы прямоугольной области a, b. Получим явным образом функцию источника f(x, y) и граничные условия. Ес- ли точное решение и приближенное с приемлемой точностью близки, то можно считать, что метод и программа работают корректно. В качестве первого тестового примера возьмјм задачу, в которой решение будет пред- ставлять собой собственную функцию задачи Штурма-Лиувилля u(x, y) = 8 ? 2 · sin ?x 4 sin ?y 4 Найдјм правую часть уравнения (2): u ? x = 8 ? 2 · ? 4 cos ?x 4 sin ?y 4 u ?? xx = ? 8 ? 2 · ? 2 16 sin ?x 4 sin ?y 4 u ? y = 8 ? 2 · ? 4 sin ?x 4 cos ?y 4 u ?? yy = ? 8 ? 2 · ? 2 16 sin ?x 4 sin ?y 4 f (x, y) = ??u(x, y) = sin ?x 4 sin ?y 4 Пусть наша область будет представлять собой квадрат 8 Ч 8 с вырезом 4 Ч 4. Тогда граничные условия u(0, y) = 0, 0 ? y < 4 u(x, 4) = 0, 0 ? x < 4 u(4, y) = 0, 4 ? y < 8 u(x, 8) = 0, 4 ? x < 8 u(8, y) = 0, 0 < y ? 8 u(x, 0) = 0, 0 < x ? 8 И задача принимает вид: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 u ?x 2 + ? 2 u ?y 2 = ? sin ?x 4 sin ?y 4 , 0 < x < 8, 0 < y < 8 u(0, y) = 0, 0 ? y < 4 u(x, 4) = 0, 0 ? x < 4 u(4, y) = 0, 4 ? y < 8 u(x, 8) = 0, 4 ? x < 8 u(8, y) = 0, 0 < y ? 8 u(x, 0) = 0, 0 < x ? 8 (5) 2.3 Построение второго тестового примера В качестве второго тестового примера возьмјм задачу, в которой решение будет иметь вид u(x, y) = e ?x · sin (?y) Найдјм правую часть уравенения (2): 5 u ? x = ?e ?x · sin (?y) u ?? xx = e ?x · sin (?y) u ? y = ? · e ?x · cos (?y) u ?? yy = ?? 2 · e ?x · sin (?y) f (x, y) = ??u(x, y) = e ?x · sin (?y)(? 2 ? 1) Наша область будет представлять собой квадрат 2 Ч 2 с вырезом 1 Ч 1. Тогда граничные условия u(0, y) = sin (?y), 0 ? y < 1 u(x, 1) = 0, 0 ? x < 1 u(1, y) = e ?1 · sin (?y), 1 ? y < 2 u(x, 2) = 0, 1 ? x < 2 u(2, y) = e ?2 · sin (?y), 0 < y ? 2 u(x, 0) = 0, 0 < x ? 2 Задача принимает вид: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 u ?x 2 + ? 2 u ?y 2 = e ?x · sin (?y)(1 ? ? 2 ), 0 < x < 2, 0 < y < 2 u(0, y) = sin (?y), 0 ? y < 1 u(x, 1) = 0, 0 ? x < 1 u(1, y) = e ?1 · sin (?y), 1 ? y < 2 u(x, 2) = 0, 1 ? x < 2 u(2, y) = e ?2 · sin (?y), 0 < y ? 2 u(x, 0) = 0, 0 < x ? 2 (6) 3 Итерационные методы решения задачи Дирихле Как следует из развјрнутой записи уравнения (4), при больших N 1 и N 2 матрица получаю- щейся системы линейных алгебраических уравнений является сильно разреженной. Поэтому следует использовать итерационные методы. Для решения задачи будем использовать метод Зейделя. Так как метод является неявным, то целесообразно нумеровать точки сетки так, чтобы матрица системы сеточных уравнений имела наиболее простой вид. Одним из таких способов является нумерация точек по рядам, начиная с выбранной угловой приграничной точки. Если за начало принять точку (x 1 , y 1 ) и занумеровать, например так, как схематично изображено в виде матрицы на рис. 3 (по горизонтальным рядам), то двумерная задача может рассматриваться как одномерная. Рис. 3: Нумерация узлов сетки 6 Для построения метода Зейделя примем за основу нумерацию, приведјнную на рис. 3. Тогда при вычислении очередного приближения U i,j значения с индексами i + 1 и j ? 1 будут известными (уже вычислены или являются граничными условиями), а значения с индексами i?1 и j+1 берјм с предыдущей итерации. Поэтому имеем следующую запись метода Зейделя: U (k+1) i+1,j ? 2U (k+1) i,j + U (k) i?1,j h 2 1 + U (k+1) i,j?1 ? 2U (k+1) i,j + U (k) i,j+1 h 2 2 = ?F i,j Примем h 1 = h 2 = h . Тогда расчјтная формула для программирования будет иметь вид: U (k+1) i,j = 1 4 · (F i,j h 2 + U (k+1) i+1,j + U (k+1) i,j?1 + U (k) i?1,j + U (k) i,j+1 ) В качестве критерия окончания будет использована формула ?U (k+1) ? U (k) ? 1 ? ? Циклы будут организованы таким образом, что вырез пластины будет пропускаться. 3.1 Решение первого тестового примера Представлен график численного решения первого тестового примера (5), h = 0.5. С помощью метода Зейделя решение с точностью ? = 10 ?2 найдено за 34 итерации. Рис. 4: Температурная карта численного решения первого тестового примера Рис. 5: Трјхмерный график численного решения первого тестового примера 7 Представлен график точного решения по тем же узлам, что и приближјнное. Рис. 6: Трјхмерный график точного решения первого тестового примера Визуально решения совпадают. Построим график погрешности приближјнного решения: Рис. 7: Трјхмерный график погрешности численного решения первого тестового примера Оба решения с приемлемой точностью близки, поэтому можно говорить о корректности разностной схемы и реализации метода Зейделя. 8 3.2 Решение второго тестового примера Представлен график приближјнного решения второго тестового примера (6), h = 0.1. С помощью метода Зейделя решение с точностью ? = 10 ?2 найдено за 49 итераций. Рис. 8: Температурная карта приближјнного решения второго тестового примера Рис. 9: Трјхмерный график приближјнного решения второго тестового примера Представлен график точного решения. Рис. 10: Трјхмерный график точного решения второго тестового примера Построим график погрешности приближјнного решения: 9 Рис. 11: Трјхмерный график погрешности приближјнного решения второго тестового при- мера Можно сделать вывод, аналогичный выводу из пункта 3.1. 3.3 Решение исходной задачи Решать поставленную задачу методом Зейделя будем аналогичным образом, что и тестовые примеры. Находить погрешность будем по правилу Рунге. Представлен график приближјнного решения задачи из условия (1), h = 1. С помощью метода Зейделя решение с точностью ? = 10 ?1 найдено за 3697 итераций. Рис. 12: Температурная карта численного решения задачи 10 Рис. 13: Трјхмерный график приближјнного решения задачи Хорошо виден максимум решения на границе, на которой расположен источника тепла. Видно убывание температуры от этих границ к другим границам области. Расположим источник тепла на границах ? 5 , ? 6 и посмотрим на получившееся решение при таких условиях. Рис. 14: Численное решение задачи с изменјнными граничными условиями Снова виден максимум на границах, на которых расположен источник, и убывание тем- пературы по мере удаления от этих границ. 11 Представлен трјхмерный график погрешности приближјнного решения исходной задачи. Оценка погрешности произведена на основе правила Рунге ? h i,j = |u h i,j ?u h/2 i,j | 2 p ?1 , p = 2 Рис. 15: Трјхмерный график погрешности приближјнного решения задачи Время нахождения решения данной задачи с точностью ? = 10 ?1 ? около 5 минут. Это очень длительный промежуток времени и хочется найти решение гораздо быстрее. Для этого будем пользоваться более быстрыми проекционными методами численного решения сеточ- ного уравнения Лапласа. 12 4 Проекционные методы решения задачи Дирихле Мы решаем СЛАУ Ax = b, причјм, исходя из схемы, A ? R nЧn и A = A T > 0 . Поэтому для того чтобы ускорить нахождение решения исходной задачи можем применять метод сопряжјнных градиентов (CG). Метод удобен тем, что нет необходимости явно формировать матрицу системы A. Для того, чтобы реализовать метод, достаточно уметь вычислять произведение матрицы системы на произвольный массив, пользуясь формулами из левой части системы ? U i?1,j ? 2U i,j + U i+1,j h 2 1 + U i,j?1 ? 2U i,j + U i,j+1 h 2 2 = F i,j В случае однородных граничных условий всј довольно просто: в каждом внутреннем узле значение координаты вектора b вычисляется по правой части приведенной системы, значение координаты вектора Ax ? по левой части. В случае неоднородных граничных условий, ненулевые граничные условия переносят- ся в правые части уравнений, таким образом изменяется вычисление координат вектора b. Уравнения системы, соответствующие всем приграничным узлам слева и справа, имеют вид ? ?2U 1,j + U 2,j h 2 1 + U 1,j?1 ? 2U 1,j + U 1,j+1 h 2 2 = F 1,j + g(0, y j ) h 2 1 , ? U N 1 ?2,j ? 2U N 1 ?1,j h 2 1 + U N 1 ?1,j?1 ? 2U N 1 ?1,j + U N 1 ?1,j+1 h 2 2 = F N 1 ?1,j + g(X, y j ) h 2 1 , j = 2, . . . , N 2 ? 2. Уравнения системы, соответствующие всем приграничным узлам снизу и сверху, имеют вид ? U i?1,1 ? 2U i,1 + U i+1,1 h 2 1 + ?2U i,1 + U i,2 h 2 2 = F i,1 + g(x i , 0) h 2 2 , ? U i?1,N 2 ?1 ? 2U i,N 2 ?1 + U i+1,N 2 ?1 h 2 1 + U i,N 2 ?2 ? 2U i,N 2 ?1 h 2 2 = F i,N 2 ?1 + g(x i , Y ) h 2 2 , i = 2, . . . , N 1 ? 2. В узле сетки (x 1 , y 1 ) уравнение принимает вид ? ?2U 1,1 + U 2,1 h 2 1 + ?2U 1,1 + U 1,2 h 2 2 = F 1,1 + g(0, y 1 ) h 2 1 + g(x 1 , 0) h 2 2 В узлах сетки (x 1 , y N 2 ?1 ), (x N 1 ?1 , y 1 ), (x N 1 ?1 , y N 2 ?1 ) уравнения меняются аналогичным об- разом. После изменения правой части системы, еј можно решать как систему для однородной задачи. Приведјм алгоритм метода сопряжјнных градиентов: 13 Метод сопряженных градиентов (CG) 1: r 0 := b ? Ax 0 ; p 0 := r 0 2: for j = 0, 1, 2, . . . do 3: ? := (r j ,r j ) (Ap j ,p j ) 4: x j+1 := x j + ?p j 5: r j+1 := r j ? ?Ap j 6: ? := (r j+1 ,r j+1 ) (r j ,r j ) 7: p j+1 := r j+1 + ?p j 8: if достигнута заданная точность then 9: выйти из цикла 10: end if 11: end for В качестве критерия окончания будет использована формула ?U (k+1) ? U (k) ? 1 ? ? В качестве начального приближения U 0 во всех случаях, как и при решении итерацион- ными методами, возьмјм нулевой вектор. Решения методом сопряженных градиентов получены и они совпадают с теми решениями, которые получены методом Зейделя. Приведјм сводную таблицу результатов: Задача h ? Метод Зейделя Метод CG (5) 0.5 10 ?2 34 2 (6) 0.1 10 ?2 49 26 (1) 1 10 ?1 3697 213 Видно, что применение метода CG позволило ускорить нахождение решения для исходной задачи в 17 раз (по количеству итераций). 5 Заключение Мы рассмотрели некоторые методы для решения задачи Дирихле в области сложной формы. В целом, для таких задач очень хорошо подходит итерационный метод Зейделя. Однако, ко- гда требуется получить четкий портрет решения, мы уменьшаем шаг и сталкиваемся с тем, что метод Зейделя начинает сходится очень медленно. В таком случае для ускорения на- хождения решения применяют другие методы ? проекционные. Один из самых знаменитых проекционных методов ? метод сопряженных градиентов. Он отлично решает нашу задачу, намного быстрее, чем это делает итерационный метод. Одна итерация проекционного метода по вычислительной трудојмкости приблизительно сравнима с итерацией по итерационному методу, что дајт явное преимущество при выборе метода для решения сложной задачи, в которой возникает матрица системы большой размерности. 14 6 Список литературы 1. Амосов, А. А. Вычислительные методы: учебное пособие / А. А. Амосов, Ю. А. Дубин- ский, Н. В. Копченова. С-Пб.: Издательство ѕЛаньї, 2012. 672с. 2. Казјнкин, К.О. Численное решение задач математической физики. Стационарные урав- нения (уравнение Пуассона) : учебно-методическое пособие / К.О. Казјнкин, О.А. Амо- сова. М.: Издательство МЭИ, 2017. 36с. 3. Казјнкин, К.О. Численное решение стационарных уравнений математической физики: методические указания / К.О. Казјнкин, А.Е. Вестфальский, О.А. Амосова. М.: Издательство МЭИ, 2018. 48с. 4. Баландин, М. Ю. Методы решения СЛАУ болльшой размерности / М. Ю. Баландин, Э. П. Шурина. Новосибирск: Издательство НГТУ, 2000. 65с. 5. Самарский, А. А. Численные методы: учебное пособие / А. А. Самарский, А. В. Гулин. М.: Наука, 1989. 432с. 15 16 ?????????? 1. import numpy as np 2. from matplotlib import pyplot as plt 3. import seaborn as sns 4. 5. # ????? ??????? 6. # 1-?? ???????? ?????? 7. 8. A = 8 9. B = 4 10. C = 4 11. 12. def F(x , y): # ?????? ????? 13. return np.sin(np.pi * x * 0.25 ) * np.sin(np.pi * y * 0.25 ) 14. 15. def seidel(eps = 10 ** (- 2 ) , h = 0.5 ): 16. vert_1 = int (C / h - 1.0 ) 17. vert_2 = int (B / h) 18. start_1 = 1 19. start_2 = int ( 0.5 * (A / h) + 1 ) 20. T = np.zeros(( int ((B + C) / h + 1 ) , int (A // h + 1 ))) 21. T_ = np. copy (T) # T_ - ???????-????????? ?? ?????????? ???????? 22. it = 1 23. while True : 24. for i in reversed ( range (T_.shape[ 0 ] - 1 - vert_1 , T_.shape[ 0 ] - 1 )): 25. for j in range (start_1 , T_.shape[ 1 ] - 1 ): 26. T[i][j] = 0.25 * (h ** 2 * F(j * h , (B + C - i) * h) + T[i][j - 1 ] + T[i + 1 ][j] + T_[i][j + 1 ] + T_[i - 1 ][j]) 27. 28. for i in reversed ( range ( 1 , vert_2 + 1 )): 29. for j in range (start_2 , T_.shape[ 1 ] - 1 ): 30. T[i][j] = 0.25 * (h ** 2 * F(j * h , (B + C - i) * h) + T[i][j - 1 ] + T[i + 1 ][j] + T_[i][j + 1 ] + T_[i - 1 ][j]) 31. 32. if np.linalg.norm(np.subtract(T , T_) , ord = 1 ) < eps: 33. break 34. else : 35. it + = 1 36. T_ = np. copy (T) 37. 38. ( '?????????? ???????? ' , it) 39. 40. return T 41. 42. # 2-?? ???????? ?????? 43. # ????? ??????? T ???????????????? ?? ??????? ????????, ? ????????-????????????? ???????? 44. 45. A = 2 46. B = 1 47. C = 1 48. 49. def F(x , y): # ?????? ????? 50. return np.exp(-x) * np.sin(np.pi * y) * (np.pi ** 2 - 1 ) 51. 52. def g(y): # ????????? ??????? 53. return np.sin(np.pi * y) 54. 55. def initial(h = 0.1 ): 56. size_1 = int ((B + C) / h + 1 ) 57. size_2 = int (A / h + 1 ) 58. mtx = np.zeros((size_1 , size_2)) 59. for i in reversed ( range ( int (B / h + 1 ) , size_1 - 1 )): 60. mtx[i][ 0 ] = g(i * h) 61. for j in range ( 1 , int ((A / 2 ) / h + 1 )): 62. mtx[ int (B / h)][j] = 0 63. for i in reversed ( range ( 1 , int (B / h))): 64. mtx[i][ int ((A / 2 ) / h)] = g(i * h) / np.exp( 1 ) 17 65. for j in range ( int ((A / 2 ) / h + 1 ) , size_2 - 1 ): 66. mtx[ 0 ][j] = 0 67. for i in range ( 1 , size_1 - 1 ): 68. mtx[i][size_2 - 1 ] = g(i * h) / np.exp( 2 ) 69. for j in range ( 1 , size_2 - 1 ): 70. mtx[size_1 - 1 ][j] = 0 71. return mtx 72. 73. # ???????? ?????? 74. # ????????? ? ??????? initial, ??????????? F(x, y) 75. 76. A = 150 77. B = 75 78. C = 45 79. Gamma_1 = Gamma_4 = 28 # ????????? ??????? 80. Gamma_2 = Gamma_3 = 70 81. Gamma_5 = Gamma_6 = 25 82. 83. def initial(h = 1 ): 84. size_1 = int ((B + C) / h + 1 ) 85. size_2 = int (A / h + 1 ) 86. mtx = np.zeros((size_1 , size_2)) 87. mtx[size_1 - 2 : int (B / h) : - 1 , 0 ] = Gamma_1 88. mtx[ int (B / h) , : int ((A / 2 ) / h + 1 )] = Gamma_2 89. mtx[ int (B / h - 1 ) : : - 1 , int ((A / 2 ) / h)] = Gamma_3 90. mtx[ 0 , int ((A / 2 ) / h + 1 ) : size_2 - 1 ] = Gamma_4 91. mtx[ 1 : size_1 , size_2 - 1 ] = Gamma_5 92. mtx[size_1 - 1 , 1 : size_2 - 1 ] = Gamma_6 93. mtx[size_1 - 1 , 0 ] = mtx[ 0 , size_2 - 1 ] = (Gamma_1 + Gamma_6) / 2 94. return mtx 95. 96. # ????? CG 97. # 1-?? ???????? ?????? 98. 99. A = 8 100. B = 4 101. C = 4 102. h = 0.5 103. 104. vert_1 = int (C / h - 1 .) # ?????? 105. vert_2 = int (B / h) 106. start_1 = 1 107. start_2 = int ( 0.5 * (A / h) + 1 .) 108. 109. def F(x , y): 110. return np.sin(np.pi * x * 0.25 ) * np.sin(np.pi * y * 0.25 ) 111. 112. def A_(P): # ??????? ???????, ?????????? ?? ???????-?????? P (A*P) 113. P_ = np.zeros(( int ((B + C) / h + 1 ) , int (A / h + 1 ))) 114. for i in reversed ( range (P.shape[ 0 ] - 1 - vert_1 , P.shape[ 0 ] - 1 )): 115. for j in range (start_1 , P.shape[ 1 ] - 1 ): 116. P_[i][j] = (- 1 ) * (h ** (- 2 )) * (P[i - 1 ][j] + P[i + 1 ][j] + P[i][j + 1 ] + P[i][j - 1 ] - 4 * P[i][j]) 117. 118. for i in reversed ( range ( 1 , vert_2 + 1 )): 119. for j in range (start_2 , P.shape[ 1 ] - 1 ): 120. P_[i][j] = (- 1 ) * (h ** (- 2 )) * (P[i - 1 ][j] + P[i + 1 ][j] + P[i][j + 1 ] + P[i][j - 1 ] - 4 * P[i][j]) 121. 122. return P_ 123. 124. def b(): # ?????? b 125. mtx = np.zeros(( int ((B + C) / h + 1 ) , int (A / h + 1 ))) 126. for i in reversed ( range (mtx.shape[ 0 ] - 1 - vert_1 , mtx.shape[ 0 ] - 1 )): 127. for j in range (start_1 , mtx.shape[ 1 ] - 1 ): 128. mtx[i][j] = F(i * h , j * h) 129. 130. for i in reversed ( range ( 1 , vert_2 + 1 )): 131. for j in range (start_2 , mtx.shape[ 1 ] - 1 ): 132. mtx[i][j] = F(i * h , j * h) 18 133. 134. return mtx 135. 136. def CG(eps = 10 ** (- 2 )): 137. T = np.zeros(( int ((B + C) / h + 1 ) , int (A / h + 1 ))) 138. T_ = np. copy (T) 139. r0 = b() 140. p = np. copy (r0) 141. it = 1 142. while True : 143. tmp = A_(p) 144. alp = np. sum (r0 * r0) / np. sum (tmp * p) 145. T = T_ + alp * p 146. r = r0 - alp * tmp 147. bett = np. sum (r * r) / np. sum (r0 * r0) 148. p = r + bett * p 149. 150. if np.linalg.norm(np.subtract(T , T_) , ord = 1 ) < eps: 151. break 152. else : 153. it + = 1 154. T_ = np. copy (T) 155. r0 = np. copy (r) 156. 157. ( '?????????? ???????? ' , it) 158. 159. return T 160. 161. # 2-?? ???????? ?????? 162. # ?????????: ??????? T ???????????????? ?? ??????? ????????, ? ????????-????????????? ????????, 163. # ????????? ?????????? ????????? ????????? ???????, ?? ?????????? ??????? b() 164. 165. A = 2 166. B = 1 167. C = 1 168. h = 0.1 169. vert_1 = int (C / h - 1 .) 170. vert_2 = int (B / h) 171. start_1 = 1 172. start_2 = int ( 0.5 * (A / h) + 1 .) 173. 174. def g(y): 175. return np.sin(np.pi * y) 176. 177. def F(x , y): 178. return np.exp(-x) * np.sin(np.pi * y) * (np.pi ** 2 - 1 ) 179. 180. def initial(h = 0.1 ): 181. size_1 = int ((B + C) / h + 1 ) 182. size_2 = int (A / h + 1 ) 183. mtx = np.zeros((size_1 , size_2)) 184. for i in reversed ( range ( int (B / h + 1 ) , size_1 - 1 )): 185. mtx[i][ 0 ] = g(i * h) 186. for j in range ( 1 , int ((A / 2 ) / h + 1 )): 187. mtx[ int (B / h)][j] = 0 188. for i in reversed ( range ( 1 , int (B / h))): 189. mtx[i][ int ((A / 2 ) / h)] = g(i * h) / np.exp( 1 ) 190. for j in range ( int ((A / 2 ) / h + 1 ) , size_2 - 1 ): 191. mtx[ 0 ][j] = 0 192. for i in range ( 1 , size_1 - 1 ): 193. mtx[i][size_2 - 1 ] = g(i * h) / np.exp( 2 ) 194. for j in range ( 1 , size_2 - 1 ): 195. mtx[size_1 - 1 ][j] = 0 196. return mtx 197. 198. def b(): 199. mtx = np.zeros(( int ((B + C) / h + 1 ) , int (A / h + 1 ))) 200. for i in reversed ( range (mtx.shape[ 0 ] - 1 - vert_1 , mtx.shape[ 0 ] - 1 )): 201. for j in range (start_1 , mtx.shape[ 1 ] - 1 ): 202. if j == start_1: 19 203. mtx[i][j] = F(j * h , i * h) + g(i * h) / h ** 2 204. elif j == mtx.shape[ 1 ] - 2 : 205. mtx[i][j] = F(j * h , i * h) + np.exp(- 2 ) * g(i * h) / h ** 2 206. else : 207. mtx[i][j] = F(j * h , i * h) 208. 209. for i in reversed ( range ( 1 , vert_2 + 1 )): 210. for j in range (start_2 , mtx.shape[ 1 ] - 1 ): 211. if j == start_2: 212. mtx[i][j] = F(j * h , i * h) + np.exp(- 1 ) * g(i * h) / h ** 2 213. elif j == mtx.shape[ 1 ] - 2 : 214. mtx[i][j] = F(j * h , i * h) + np.exp(- 2 ) * g(i * h) / h ** 2 215. else : 216. mtx[i][j] = F(j * h , i * h) 217. 218. return mtx 219. 220. # ???????? ?????? 221. # ????????? ? ??????? initial, ??????????? F(x, y), ???????? ? ??????? b(), ????????? ??? ????????? ??????? ????????? 222. 223. A = 150 224. B = 75 225. C = 45 226. Gamma_1 = Gamma_4 = 28 # ????????? ??????? 227. Gamma_2 = Gamma_3 = 70 228. Gamma_5 = Gamma_6 = 25 229. h = 1 230. 231. vert_1 = int (C / h - 1 .) 232. vert_2 = int (B / h) 233. start_1 = 1 234. start_2 = int ( 0.5 * (A / h) + 1 .) 235. 236. def initial(): 237. size_1 = int ((B + C) / h + 1 ) 238. size_2 = int (A / h + 1 ) 239. mtx = np.zeros((size_1 , size_2)) 240. mtx[size_1 - 2 : int (B / h) : - 1 , 0 ] = Gamma_1 241. mtx[ int (B / h) , : int ((A / 2 ) / h + 1 )] = Gamma_2 242. mtx[ int (B / h - 1 ) : : - 1 , int ((A / 2 ) / h)] = Gamma_3 243. mtx[ 0 , int ((A / 2 ) / h + 1 ) : size_2 - 1 ] = Gamma_4 244. mtx[ 1 : size_1 , size_2 - 1 ] = Gamma_5 245. mtx[size_1 - 1 , 1 : size_2 - 1 ] = Gamma_6 246. mtx[size_1 - 1 , 0 ] = mtx[ 0 , size_2 - 1 ] = (Gamma_1 + Gamma_6) / 2 247. return mtx 248. 249. def b(): 250. mtx = np.zeros(( int ((B + C) / h + 1 ) , int (A / h + 1 ))) 251. for i in reversed ( range (mtx.shape[ 0 ] - 1 - vert_1 , mtx.shape[ 0 ] - 1 )): 252. for j in range (start_1 , mtx.shape[ 1 ] - 1 ): 253. if j == start_1 and i == mtx.shape[ 0 ] - 2 : 254. mtx[i][j] = (Gamma_1 + Gamma_6) / h ** 2 255. elif j == start_1 and i == mtx.shape[ 0 ] - vert_1 - 1 : 256. mtx[i][j] = (Gamma_1 + Gamma_2) / h ** 2 257. elif j == mtx.shape[ 1 ] - 2 and i == mtx.shape[ 0 ] - 2 : 258. mtx[i][j] = (Gamma_6 + Gamma_5) / h ** 2 259. elif i == mtx.shape[ 0 ] - 2 : 260. mtx[i][j] = Gamma_6 / h ** 2 261. elif i == mtx.shape[ 0 ] - vert_1 - 1 and j < start_2: 262. mtx[i][j] = Gamma_2 / h ** 2 263. elif j == start_1: 264. mtx[i][j] = Gamma_1 / h ** 2 265. elif j == mtx.shape[ 1 ] - 2 : 266. mtx[i][j] = Gamma_5 / h ** 2 267. else : 268. mtx[i][j] = 0 269. 270. for i in reversed ( range ( 1 , vert_2 + 1 )): 271. for j in range (start_2 , mtx.shape[ 1 ] - 1 ): 20 272. if j == start_2 and i == 1 : 273. mtx[i][j] = (Gamma_3 + Gamma_4) / h ** 2 274. elif j == mtx.shape[ 1 ] - 2 and i == 1 : 275. mtx[i][j] = (Gamma_4 + Gamma_5) / h ** 2 276. elif i == 1 : 277. mtx[i][j] = Gamma_4 / h ** 2 278. elif j == start_2: 279. mtx[i][j] = Gamma_3 / h ** 2 280. elif j == mtx.shape[ 1 ] - 2 : 281. mtx[i][j] = Gamma_5 / h ** 2 282. else : 283. mtx[i][j] = 0 284. 285. return mtx |