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

    Введение


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

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

    Подробнее
    Реклама

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

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

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

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

        –1
        Ну, справедливости ради, бывает так, что релаксация может оказаться не хуже, а то и лучше прочих продвинутых методов (в нелинейных задачах правда). Я, во всяком случае, такое встречал.
          –1
          Ну вот пусть автор публикации и сравнит. Хоть на своем примере, хоть на том, что из моей ссылки.
        –1
        Нехорошо как-то граничные условия Неймана аппроксимированы в одномерной задаче, в двумерной, я так полагаю, кругом Дирихле? Второе, что такое «квадратичная матрице {a}»? Ну и центральные разности для конвекции в сочетании с релаксацией… Не очень понял, кому предназначена эта публикация? Школьникам старших классов? Студентам первого курса?
          –1
          У Петра Николаевича Вабищевича защищался как раз по этой теме много лет назад. Эх, ностальгия

          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

          Самое читаемое