8.3. Решение дифференциальных уравнений эллиптического типа

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

Из электростатики известно, что распределение потенциала удовлетворяет уравнению Пуассона.

, (8.24)

которое в декартовой системе координат будет иметь вид

. (8.25)

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

(8.26)

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

.

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

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

, (8.28)

(8.29)

где и– шаг сетки по координатами, соответственно, а. Откуда

. (8.30)

Это один из простейших шаблонов, часто называемый шаблоном типа «крест», будет иметь вид:

(i,j+1)

|

(i-1,j) — (i,j) — (i+1,j)

|

(i,j-1)

После преобразований (8.30) получаем систему уравнений относительно :

, (8.31)

где . Для квадратной сетки, которая наиболее удобна для вычислений, формула упрощается, и ее знаменатель будет равен четырем. Система (8.31) является неявной и обычно решается итерационными методами, например методом Зейделя. Итерационный процесс завершается при выполнении в каждом узле сетки условия

, (8.32)

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

В качестве примера реализации программы рассмотрим случай соответствующий рисунку 8.2 при , а длина стенок равна 1.

Программа 8.2

h=0.01; %Задание шага сетки.

x=0:h:1;

y=0:h:1;

u=zeros(101,101); %Задание начальных значений потенциала.

u(:,101)=10; %Задание граничного условия.

for k=1:20, %k – число итераций.

for i=2:100,

for j=2:100,

u(i,j)=(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))/4;

end;

end;

end;

meshgrid(x,y); %Построение поверхности решения.

mesh(u);