10 Решение двумерной краевой задачи методом ритца

Пусть требуется найти решение уравнения

, (10.1)

называемого уравнением Пуассона, внутри двумерной области , принимающее заданные значения на границе  (так называемая задача Дирихле).

В вариационном исчислении показано, что эта краевая задача равносильна задаче об исследовании на минимум функционала

(10.2)

на множестве функций u, удовлетворяющих на границе  области  граничному условию 1-го рода .

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

, (10.3)

удовлетворяющих граничному условию, то задача существенно упростится. Теперь требуется всего лишь решить задачу на экстремум функции нескольких переменных. Действительно, считая _i неизвестными параметрами (числами), а N_i – известными функциями координат (базисными функциями), подставим аппроксимацию u*(x,y) в функционал (2):

.

Как видно, функционал F после такой подстановки превратился в функцию  переменных ₁,₂,… ,_m. Так как для существования экстремума функции n переменных необходимо выполнение условий

,

получим систему линейных алгебраических уравнений относительно параметров {_i}

,

; .

Эти параметры полностью определяют аппроксимацию искомого решения согласно (10.3). Чем больше берется слагаемых в (10.3), или, что то же самое, параметров аппроксимации, тем ближе полученное по методу Ритца решение приближается к точному решению. (Это верно при условии полноты системы базисных функций {N_i}).

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

и/или условие 3-го рода

(здесь q,  и  – заданные функции координат).

В этом случае учет этих условий обеспечивается за счет минимизации функционала, отличающегося от (2) появлением дополнительных граничных слагаемых – криволинейных интегралов 1-го рода по соответствующей части границы:

и

(10.4)

соответственно.

Пример. Найти стационарное распределение температуры в прямоугольной пластине 0xа, 0yb, все стороны которой поддерживаются при постоянной температуре u=20. Пластина нагревается от источников тепла, мощность которых описывается законом Q(x,y) = x² + y².

Решение. Задача сводится к решению уравнения

(k – коэффициент теплопроводности, положим k=1) с граничными условиями u(0,y) = u(5,y) = u(x,0) = = u(x,4) = 20.

Имеем краевые условия 1-го рода вдоль всей границы области . Чтобы обеспечить их выполнение, воспользуемся аппроксимацией вида

.

При этом если за базисные функции взять произведения , i, j=1, 2, 3,…, то на границах прямоугольной области они обращаются в нуль, и, тем самым, положив =20, граничные условия будут удовлетворены точно вне зависимости от значений {_k}. Итак,

, (10.5)

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

Поскольку в задаче есть только условия 1-го рода (Дирихле), то используем функционал

Выполняя затем стандартную процедуру Ритца, найдем оптимальные значения параметров {_ij}, а вместе с тем (см. (10.5)) и приближенное решение задачи.

Порядок решения в системе Maple

1. Зададим исходные данные задачи, а также верхние границы изменения индексов i, j – числа m и n, произведения которых определяет число параметров аппроксимации. Их увеличение приводит к повышению точности решения и одновременно усложняет расчеты.

> restart;

> a:=5: b:=4: f:=x^2+y^2; m:=2; n:=2;

2. Запишем аппроксимацию решения задачи

> z:=(x,y)->20+sum(sum(alpha[i,j]*x^i*y^j*

(x-a)*(y-b),i=1..m),j=1..n);

3. Зададим функционал в виде повторного интеграла

> F:=(1/2.)*int(int(diff(z(x,y),x)^2+

diff(z(x,y),y)^2-2*f*z(x,y),y=0..b),x=0..a);

4. Формируем уравнения системы

> eqns:=seq(seq(diff(F,alpha[i,j]),i=1..m),j=1..n); > p:=seq(seq(alpha[i,j],i=1..m),j=1..n);

5. Находим решение этой системы

> r:=solve({eqns},{p});

6. Определяем функцию (x, y), хранящую решение

> phi:=unapply(subs(r,z(x,y)),x,y);

7. Теперь решение задачи можно получить в любых точках прямоугольной области, например, вдоль линии x=2:

> [phi(2,0),phi(2,1),phi(2,2),phi(2,3),phi(2,4)];

8. Визуализация решения – вывод в виде цветовой карты и линий равного уровня

> plots[densityplot](phi(x,y),x=0..a,y=0..b,

grid=[49,49],colorstyle=HUE);

> plots[contourplot](phi(x,y),x=0..a,y=0..b,

contours=15);

Рис. 10.2.

Замечания. Если в задаче имеются краевые условия 2-го или 3-го рода вдоль какой-то стороны, то аппроксимация не должна удовлетворять каким бы то ни было условиям на этой стороне. При этом в функционал необходимо включить дополнительный интеграл по данной стороне (см. (10.4)).

Изучите внимательно нижеследующие примеры.

1) u(0,y) =1, u(a,y) =2; (x,0) =1, (x,b) = –1.

Аппроксимация учитывает только два первых условия как 1-го рода

В функционале к кратному интегралу по прямоугольнику добавляются интегральные слагаемые вдоль сторон x=0 и x=a:

Изменения в программе коснутся соответствующих строк:

> z:=(x,y)->1+x/a+sum(sum(alpha[i,j]*x^i*y^j*

(x-a),i=1..m),j=0..n);

> F:=(1/2.)*int(int(diff(z(x,y),x)^2+

diff(z(x,y),y)^2-2*f*z(x,y),y=0..b),x=0..a)-

int(z(x,0),x=0..a)+int(z(x,b),x=0..a);

> eqns:=seq(seq(diff(F,alpha[i,j]),i=1..m),0=1..n);

> p:=seq(seq(alpha[i,j],i=1..m),0=1..n);

2) u(0,y) =10, u(x,0) =10; (x,b) + u(x, b) =0, (a,y)+u(a, y) =0.

Аппроксимация ;

Функционал

.