Численные методы решения уравнений эллиптического типа

от автора

Введение

Наиболее распространённым уравнением эллиптического типа является уравнение Пуассона.
К решению этого уравнения сводятся многие задачи математической физики, например задачи о стационарном распределении температуры в твердом теле, задачи диффузии, задачи о распределении электростатического поля в непроводящей среде при наличии электрических зарядов и многие другие.

Для решения эллиптических уравнений в случае нескольких измерений используют численные методы, позволяющие преобразовать дифференциальные уравнения или их системы в системы алгебраических уравнений. Точность решения опреде­ляется шагом координатной сетки, количеством итераций и разрядной сеткой компьютера [1]

Цель публикации получить решение уравнения Пуассона для граничных условий Дирихле и Неймана, исследовать сходимость релаксационного метода решения на примерах.

Уравнение Пуассона относится к уравнениям эллиптического типа и в одномерном случае имеет вид [1]:

(1)

где x – координата; u(x) – искомая функция; A(x), f(x) – некоторые непрерывные функции координаты.

Решим одномерное уравнение Пуассона для случая А = 1, которое при этом принимает вид:

(2)

Зададим на отрезке [xmin, xmax] равномерную координатную сетку с шагом ∆х:

(3)

Граничные условия первого рода (условия Дирихле) для рассматривае­мой задачи могут быть представлены в виде:

(4)

где х1, xn – координаты граничных точек области [xmin, xmax]; g1, g2 – некоторые
константы.

Граничные условия второго рода (условия Неймана) для рассматривае­мой задачи могут быть представлены в виде:

(5)

Проводя дискретизацию граничных условий Дирихле на равномерной координатной сетке (3) с использованием метода конечных разностей, по­лучим:

(6)

где u1, un – значения функции u(x) в точках x1, xn соответственно.

Проводя дискретизацию граничных условий Неймана на сетке (3), по­лучим:

(7)

Проводя дискретизацию уравнения (2) для внутренних точек сетки, по­лучим:

(8)

где ui, fi – значения функций u(x), f(x) в точке сетки с координатой xi.

Таким образом, в результате дискретизации получим систему линейных алгебраических уравнений размерностью n, содержащую n – 2 уравнения вида (8) для внутренних точек области и уравнения (6) и (7) для двух граничных точек [1].

Ниже приведен листинг на Python численного решения уравнения (2) с граничными условиями (4) – (5) на координатной сетке (3).

Листинг решения

from numpy import* from numpy.linalg import solve import matplotlib.pyplot as plt x0=0#Начальная координата области решения xn=5#Конечная координата области решения n=100#Число точек координатной сетки dx=(xn-x0)/(n-1)#Задание равномерной координатной сетки с шагом dx x=[i*dx+x0 for i in arange(0,n,1)]#Задание равномерной координатной сетки с шагом dx def f(i):#Функция правой части уравнения          return 2*sin(x[i]**2)+cos(x[i]**2) v1=1.0#Вид ГУ на левой границе (1 - Дирихле, 2 - Неймана) g1=0.0#Значение ГУ на левой границе v2=2.0#'Вид ГУ на правой границе (1 - Дирихле, 2 - Неймана) g2=-0.5#Значение ГУ на правой границе a=zeros([n,n])#Задание матрицы коэффициентов СЛАУ размерностью n x n b=zeros([1,n])# Задание матрицы-строки свободных членов СЛАУ размерностью 1 x n  #Определение коэффициентов и свободных членов СЛАУ, # соответствующих граничным условиям и проверка корректности #значений параметров v1, v2 b[0,n-1]=g1; if v1==1:          a[0,0]=1 elif v1==2:          a[0,0]=-1/dx          a[0,1]=1/dx; else:          print('Параметр v1 имеет неправильное значение') b[0,n-1]=g2; if v2==1:          a[n-1,n-1]=1          elif v2==2:          a[n-1,n-1]=1/dx          a[n-1,n-2]=-1/dx; else:          print('Параметр v2 имеет неправильное значение') #Определение коэффициентов и свободных членов СЛАУ, # соответствующих внутренним точкам области          for i in arange(1, n-1,1):          a[i,i]=-2/dx**2          a[i,i+1]=1/dx**2          a[i,i-1]=1/dx**2          b[0,i]=f(i) u=linalg.solve(a,b.T).T#Решение СЛАУ def viz(v1,v2):          if v1==v2==1:                   return "ГУ  Дирихле на левой и ГУ Дирихле на правой  границе "          elif v1==1 and v2==2:                   return "ГУ  Дирихле на левой и ГУ Неймана на правой  границе "          elif v2==1 and v2==1:                   return "ГУ  Неймана на левой и ГУ Дирихле  на правой  границе " plt.figure() plt.title("График функции правой части уравнения Пуассона") y=[f(i) for i in arange(0,n,1)] plt.plot(x,y) plt.grid(True) plt.xlabel('x') plt.ylabel('f(x)') plt.figure() plt.title("График искомой функции уравнения Пуассона") plt.xlabel('x') plt.ylabel('u(x)') plt.plot(x,u[0,:],label='%s'%viz(v1,v2)) plt.legend(loc='best') plt.grid(True) plt.show()

Получим:

Разработанная мною на Python программа удобна для анализа граничных условий.Приведенный алгоритм решения на Python использует функцию Numpy — u=linalg.solve(a,b.T).T для решения системы алгебраических уравнений, что повышает быстродействие при квадратной матрице {a}. Однако при росте числа измерений необходимо переходить к использованию трех диагональной матрицы решение для которой усложняется даже для очень простой задачи, вот нашёл на форуме такой пример:

Пример решения с трёх диагональной матрицей

from __future__ import print_function  from __future__ import division  import numpy as np  import time  ti = time.clock()  m = 1000  A = np.zeros((m, m))  B = np.zeros((m, 1))   A[0, 0] = 1  A[0, 1] = 2  B[0, 0] = 1  for i in range(1, m-1):      A[i, i-1] = 7      A[i, i] = 8       A[i, i+1] = 9      B[i, 0] = 2  A[m-1, m-2] = 3  A[m-1, m-1] = 4  B[m-1, 0] = 3  print('A \n', A)  print('B \n', B)  x = np.linalg.solve(A, B)  # solve A*x = B for x  print('x \n', x)  print('NUMPY time', time.clock()-ti, 'seconds') 

Программа численного решения на равномерной по каждому направлению сетки задачи Дирихле для уравнения конвекции-диффузии

[2]

(9)

Используем аппроксимации центральными разностями для конвективного слагаемого и итерационный метод релаксации.для зависимость скорости сходимости от параметра релаксации при численном решении задачи с /(х) = 1 и 6(х) = 0,10. В сеточной задаче:

(10)

Представим матрицу А в виде суммы диагональной, нижней треугольной и верхней треугольных матриц:

(10)

Метод релаксации соответствует использованию итерационного метода:

(11)

При \ говорят о верхней релаксации, при — о нижней релаксации.

Листинг програмы

from numpy import * """ Численное решение задачи Дирихле для уравнения конвекции-диффузии в прямоугольнике.Метод релаксации.""" def relaxation(b, f, I1, I2, n1, n2, omega, tol = 1.e-8):          h1 = I1 / n1          h2 = I2 / n2          d = 2. / h1**2 + 2. / h2**2          y = zeros([n1+1, n2+1])          ff = zeros([n1+1, n2+1])          bb = zeros([n1+1, n2+1])          for j in arange(1,n2,1):                   for i in arange(1,n1,1):                            ff [i,j] = f(i*h1, j*h2)                            bb[i,j] = b(i*h1, j*h2)                                     #максимальное число итераций - 10000          for k in arange(1, 10001,1):                   rn = 0.                   for j in arange(1,n2,1):                            for i in arange(1,n1,1):                                                                         rr = - (y[i-1,j] - 2.*y [i, j] + y[i+1,j]) / h1**2 \                                          - (y[i,j-1] - 2.*y [i,j] + y[i,j+1]) / h2**2 \                                          + bb[i,j]*(y [i+1,j] - y [i-1,j]) / (2.*h1) - ff [i,j]                                                                         rn = rn + rr**2                                     y[i,j] = y[i,j] - omega * rr / d                   rn = rn*h1*h2                   if rn < tol**2: return y, k                            print ('Метод релаксации не сходиться:')          print ('после 10000 итерации остаток=',sqrt(rn)) import matplotlib.pyplot as plt bcList = [0., 10.] sglist = ['-','--'] kk = 0 for bc in bcList:                   I1 = 1.          I2 = 1.          def f(x,y):                   return 1.          def b(x,y):                                    return bc          n1 = 25          n2 = 25          m = 20          om = linspace(1., 1.95, m)          it = zeros(([m]))          for k in arange(0,m,1):                   omega = om[k]                   y, iter = relaxation(b, f, I1, I2, n1, n2, omega, tol=1.e-6)                   it[k] = iter          s1= 'b =' + str(bc)          sg = sglist[kk]          kk = kk+1          plt.plot( om,it, sg, label = s1) plt.title("Число итераций метода релаксации\n для приближённого решения эллиптической задачи\n с использованием заданного параметра релаксации $\\omega$") plt.xlabel('$\\omega$') plt.ylabel('iterations') plt.legend(loc=0) plt.grid(True) plt.show(

)

Получим:

На графике показана зависимость числа итераций от параметра релаксации для уравнения Пуассона (b(х) = 0) и уравнения конвекции-диффузии (b(х) = 10). Для сеточного уравнения Пуассона оптимальное значении параметра релаксации находится аналитически, а итерационный метод сходиться при .

Выводы:

  1. Приведено решение эллиптической задачи на Python с гибкой системой установки граничных условий
  2. Показано что метод релаксации имеет оптимальный диапазон () параметра релаксации.

Ссылки:

  1. Рындин Е.А. Методы решения задач математической физики. – Таганрог:
    Изд-во ТРТУ, 2003. – 120 с.
  2. Вабищевич П.Н.Численные методы: Вычислительный практикум. — М.: Книжный дом
    «ЛИБРОКОМ», 2010. — 320 с.


ссылка на оригинал статьи https://habr.com/post/418981/


Комментарии

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *