Введение


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

Для решения эллиптических уравнений в случае нескольких измерений используют численные методы, позволяющие преобразовать дифференциальные уравнения или их системы в системы алгебраических уравнений. Точность решения опреде­ляется шагом координатной сетки, количеством итераций и разрядной сеткой компьютера [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 с.

Комментарии (5)


  1. Andy_U
    03.08.2018 00:14

    Странный выбор. Вообще-то, один из самых популярных методов в современных коммерческих пакетах — метод Ньютона-Крылова. Очень быстро сходится.

    Вот пример использования scipy.optimize.newton_krylov (даже без preconditioner'а) для решения аж интегро-дифференциального уравнения:

    https://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html


    1. rafuck
      03.08.2018 00:50

      Ну, справедливости ради, бывает так, что релаксация может оказаться не хуже, а то и лучше прочих продвинутых методов (в нелинейных задачах правда). Я, во всяком случае, такое встречал.


      1. Andy_U
        03.08.2018 00:59

        Ну вот пусть автор публикации и сравнит. Хоть на своем примере, хоть на том, что из моей ссылки.


  1. rafuck
    03.08.2018 00:26
    -1

    Нехорошо как-то граничные условия Неймана аппроксимированы в одномерной задаче, в двумерной, я так полагаю, кругом Дирихле? Второе, что такое «квадратичная матрице {a}»? Ну и центральные разности для конвекции в сочетании с релаксацией… Не очень понял, кому предназначена эта публикация? Школьникам старших классов? Студентам первого курса?


  1. ooby
    03.08.2018 08:00

    У Петра Николаевича Вабищевича защищался как раз по этой теме много лет назад. Эх, ностальгия