7.6. Граничная задача для дифференциального уравнения в частных производных

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

Рассмотрим в качестве примера задачу [64] об определении напряженности статического электрического поля в бесконечном металлическом желобе прямоугольного поперечного сечения (рис. 7.7). Верхняя стенка желоба имеет потенциал остальные стенки заземлены и находятся под нулевым потенциалом. Необходимо найти функцию описывающую рвспределение потенциала в прямоугольной области.

Рис. 7.7. Поперечное сечение металлического желоба

Из электростатики известно, что, если в замкнутом объеме отсутствуют объемные электрические заряды, то распределение потенциала удовлетворяет уравнению Лапласа

В декартовой системе координат для двумерного случая последнее уравнение принимает вид

Необходимо найти решения уравнения (7.40), удовлетворяющие граничным условиям

После определения потенциала можно найти напряженность электрического поля используя известное соотношение

Поставленная задача является частным случаем задачи Дирихле, когда необходимо найти решение дифференциального уравнения в частных производных в замкнутой области при заданном распределении искомой функции на границах этой области (рис. 7.7). В случае прямоугольной области и несложных граничных условий (7.41) можно методом разделения переменных решить уравнение Лапласа в аналитическом виде [64]

Однако решение (7.42) в отличие от численных алгоритмов не обобщается на более сложные границы и граничные условия. В некоторых задачах аналитические решения настолько громоздки, что невозможно осуществить их анализ без проведения численных расчетов, в этом случае численные методы могут оказаться более предпочтительными вследствие их универсальности.

Для представления производных, входящих в уравнение и в граничные условия, в области изменения переменных вводят некоторую сетку. В двумерных областях наиболее часто используют прямоугольную сетку (рис. 7.8), разбивая область вдоль координаты х на частей, а вдоль координаты у - на частей. Приближенные значения производных в каждом узле выбранной сетки записывают, используя информацию о значениях искомой функции в соседних узлах. Конфигурацию расположения информационных узлов называют шаблоном данной разностной схемы. Один из простейших шаблонов на прямоугольной сетке - шаблон типа «крест» (рис. 7.9).

Составление разностных схем на заданном шаблоне осуществляется тремя основными способами: методом разностной аппроксимации, интегро-интерполяционным методом и методом неопределенных коэффициентов [1].

Разностную схему уравнения Лапласа запишем с помощью первого из названных методов, который заключается в замене производной в каждом узле разностным выражением, включающим только узлы шаблона.

Рис. 7.8. Область поиска решений задачи Дирихле

Рис. 7.9. Шаблон типа «крест»

По аналогии со вторым соотношением (7.5) представим вторые частные производные на шаблоне типа «крест» в виде

где и I - шаги сетки по координатам х

После подстановки конечно-разностных представлений производных в уравнение Лапласа (7.40) получим систему линейных алгебраических уравнений относительно искомых значений потенциала у в узлах. Вследствие большой размерности полученная система уравнений совместно с граничными условиями (7.41) решается итерационным методом Зейделя . Для осуществления итерационного процесса конечно-разностную схему уравнения Лапласа представляем в виде

где

Для квадратной сетки формула (7.43) упрощается, и ее знаменатель будет равен четырем, в этом случае итерационное соотношение называют формулой [65].

Итерационный процесс по формуле (7.43) завершается при выполнении в каждом узле сетки условия

где заданная относительная погрешность вычисления потенциала

Общая погрешность решения задами Дирихле обусловлена тремя составляющими: 1) погрешностью аппроксимации частных производных; 2) погрешностью представления граничных условий при несовпадении узлов сетки с границей области погрешностью решения системы разностных уравнений методом Зейделя. В большинстве задач основной вклад в общую погрешность вносит погрешность аппроксимации производных. При известном порядке последней погрешности можно провести уточнение результатов сгущением сетки по формулам Рунге.

Для уменьшения погрешности аппроксимации производных используют шаблоны с болчшим количеством узлов по сравнению с шаблоном типа «крест». Так, используя шаблон с девятью узлами (рис. 7.10), вторые производные аппроксимируют с погрешностью четвертого порядка и на квадратной сетке получают итерационную формулу [65]

которую называют формулой

Рис. 7.10. Шаблон для формулы

Рис. 7.11. Шаблоны для формул

С более крупным шаблоном (рис. 7.11) получают так называемую формулу [65]

Если потенциал скачком изменяется в углах за счет граничных условий, то рекомендуется исключить из шаблона (рис. 7.10) угловые точки. 13-16 и использовать формулу

В [65] со ссылкой на [66] отмечается, что при одинаковой погрешности граничной задачи для уравнения Лапласа по формуле получается в пять раз быстрее, чем по формуле и в два раза быстрее, чем по формуле

Рис. 7.12. Блок-схема программы решения задачи Дирихле методом конечных разностей

Структура программ 7.5 представлена на рис. 7.12. В основном блоке Каждой из программ 7.5 в диалоговом режиме осуществляется ввод следующих переменных: условное число, при нулевом значении которого Осуществляется расчет потенциала по формуле (7.42); размеры системы; погрешность решения системы разностных уравнений; потенциал верхней стенки; число разбиений области вдоль координат х и у. Затем вычисляются шаги сетки по координатам х и у, а также значение переменной необходимое в подпрограммах. После Обращения к подпрограммам аналитического или разностного методов юуществляется вывод результатов на дисплей.

В блоке 1 реализована подпрограмма расчета распределения потенциала в заданной прямоугольной области по формуле (7.42). Если слагаемые суммы (7.42) вычислять непосредственно с использованием гиперболических функций, то погрешность расчета по аналитической формуле будет превышать погрешность разностного алгоритма даже со сравнительно крупной сеткой. Поэтому в формуле (7.42) переходим к вычислению отношения гиперболических синусов в виде

где что значительно уменьшает погрешность округления. Выражения, не зависящие от индекса суммирования к, вычисляются вне цикла для накопления суммы. Суммирование прекращается, когда абсолютная величина отношения очередного слагаемого к накопленной сумме оказывается меньше заданной погрешности. Вследствие симметрии системы распределение потенциала вычисляется только в точках левой половины прямоугольной области (рис. 7.6). Переменная К введена для получения изменяющейся величины 2 к входящей в аргументы функций. Результаты вычисления потенциала в узлах размещаются в двумерном массиве

В блоке 2 задаются граничные условия и нулевые начальные приближения к потенциалу в каждом узле выбранной сетки. Затем с помощью

вложенных циклов организован итерационный метод Зейделя по формуле (7.43). Программная реализация окончания процесса уточнения решений по условиям (7.44) осуществлена с помощью способа, использованного в программах 2.2.

Чтобы не нарушать общности алгоритма решения уравнения Лапласа, в приведенных подпрограммах блока 2 не учитываются свойства симметрии системы, что позволит перейти к более сложным границам и неоднородным граничным условиям. Если же при подобных усложнениях задачи сохраняются геометрическая и «электрическая» симметрии, то можно почти вдвое сократить объем вычислений, а также размеры массива для значений потенциала в узлах сетки. Для реализации таких возможностей в программе достаточно изменить следующие строки:

и в тело цикла по переменной включить условие симметрии

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

Таблица 7.2 (см. скан)

В программе блоки 1 и 2 оформлены в виде подпрограмм, размещенных в строках 100-190 и 200-390 соответственно.

В программах расчет потенциала в узлах по формуле (7.42) осуществляется в подпрограммах с именем входными параметрами которых являются геометрические размеры системы, размеры сетки, допустимая погрешность, а выходным параметром - двумерный массив для значений потенциала в узлах сетки. Разностный метод реализован в подпрограммах с именем

В табл. 7.2 приведена выборка результатов расчета потенциала при следующих значениях параметров: Вследствие симметрии системы в таблице приведены значения потенциала в узлах левой половины прямоугольной области.

(см. скан)

(см. скан)

(см. скан)