Пример 4. Уравнение Эйлера для двойного интеграла

Решить вариационную задачу для функционала, если на контуре области функция принимает заданные значения:

Решение. Будем решать задачу в среде Maple. Определим подынтегральную функцию. Для удобства обозначим производные функции по и по соответственно, через и :

> F:=proc(x,y,u,ux,uy) options operator, arrow;

ux^2+uy^2+2*y*u*cos(x)

end proc;

Уравнение Эйлера в этих обозначениях имеет вид:

Вычисляем последовательно производные в этом уравнении

> a1:=diff(F(x,y,u,ux,uy),u);

a2:=diff(F(x,y,u,ux,uy),ux);

a3:=diff(F(x,y,u,ux,uy),uy);

Составляем уравнение Эйлера для функционала

> EulerEq:=a1-(diff(subs(ux=diff(u(x,y),x),a2),x))-

(diff(subs(uy=diff(u(x,y),y),a3),y))=0;

Таким образом, получили следующее уравнение

> pde:=-(1/2)*op(2,lhs(EulerEq))-

(1/2)*op(3,lhs(EulerEq))=(1/2)*op(1,lhs(EulerEq));

Сформируем теперь граничные условия. Для этого определим граничную функцию

> f:=proc(x,y) options operator, arrow;

(1/10)*x+(1/50)*y^2

end proc;

Определяем граничные условия:

> bc:=u(0,y)=f(0,y),u(1,y)=f(1,y),

u(x,0)=f(x,0),u(x,2)=f(x,2);

Таким образом, задача математической физики поставлена: найти функцию , удовлетворяющую дифференциальному уравнению

в прямоугольнике , и принимающую заданные значения на границе этого прямоугольника

Сформулированная задача классифицируется как задача Дирихле для уравнения Пуассона в прямоугольнике [1]. Это — неоднородная задача, причем неоднородности присутствуют как в граничных условиях, так и в уравнении.

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

Задача 1: найти функцию , удовлетворяющую дифференциальному уравнению

и граничным условиям

Задача 2: найти функцию , удовлетворяющую дифференциальному уравнению

и граничным условиям

Определим эти задачи в Maple.

Задача 1:

> pde1:=diff(u1(x,y),x,x)+diff(u1(x,y),y,y)=

(1/2)*y*cos(x);

bc1:=u1(0,y)=0,u1(1,y)=0,

u1(x,0)=f(x,0),u1(x,2)=f(x,2);

Задача 2:

> pde2:=diff(u2(x,y),x,x)+diff(u2(x,y),y,y)=

(1/2)*y*cos(x);

bc2:=u2(0,y)=f(0,y),u2(1,y)=f(1,y),

u2(x,0)=0,u2(x,2)=0;

Очевидно, решение исходной задачи будет

> u:=proc(x,y) options operator, arrow;

u1(x, y)+u2(x, y)

end proc;

Проверим это. Подставим функцию в уравнение

> simplify(pde, {pde1, pde2});

Как видим, получили тождество. Проверим выполнение граничных условий

> bc;

Таким образом, условия на функцию тоже выполняются.

Рассмотрим задачу 1. Будем решать ее методом Гринберга [1]. Соответствующая задача Штурма-Лиувилля

очевидно, имеет решение[1]

> X:=proc(x,n) options operator, arrow;

sin(n*Pi*x)

end proc;

lambda:=proc(n) options operator, arrow;

n^2*Pi^2

end proc;

Действительно, подставим это решение в уравнение

> Diff(X(x,n),x,x)+lambda(n)*X(x,n) = 0;value(%);

Проверим выполнение граничных условий

>X(0,n)=0,X(1,n)=0;simplify(%,assume=integer);

Решение задачи 1, в соответствии с методом Гринберга, представим в виде ряда по собственным функциям задачи Штурма-Лиувилля

Здесь учтено, что квадрат нормы собственных функций равен :

Действительно,

>int(X(x, n)^2, x = 0.. 1);

simplify(%, assume = integer);

Уравнение для трансформанты получим, умножив исходное уравнение в частных производных на и проинтегрировав на отрезке :

>int(lhs(pde1)*X(x,n),x = 0.. 1) =

int(rhs(pde1)*X(x,n),x = 0.. 1);

>intPDE1:= simplify(%, assume = integer);

Преобразуем правую часть полученного уравнения

>LHS1:=IntegrationTools[Expand](lhs(intPDE1));

Первый интеграл берем по частям два раза:

>IntegrationTools[Parts](op(1, LHS1), sin(n*Pi*x));

I1:= simplify(%, assume = integer);

>IntegrationTools[Parts](I1, cos(n*Pi*x));

I1:= simplify(%, assume = integer);

Учтем граничные условия

>I1:= simplify(I1, {bc1});

Итак, первый интеграл с учетом преобразовался в выражение

> I1:= subs(op(4, I1) = U1(y, n), I1);

Второй интеграл в выражении LHS1 есть, очевидно,

> I2:=subs(op(2,LHS1)=diff(U1(y,n),y,y),op(2,LHS1));

Таким образом, уравнение для трансформанты получено:

> ode1:= I2+I1 = rhs(intPDE1);

Сформируем граничные условия

>int(lhs(bc1[3])*X(x,n),x = 0.. 1) =

int(rhs(bc1[3])*X(x,n),x = 0.. 1);

>bcODE1:=simplify({subs(lhs(%) = U1(0,n), %)},

assume = integer);

>int(lhs(bc1[4])*X(x,n),x = 0.. 1) =

int(rhs(bc1[4])*X(x,n),x = 0.. 1);

>simplify({subs(lhs(%) = U1(2,n), %)},

assume = integer);

>bcODE1:= `union`(bcODE1, %);

Находим трансформанту

> res:= dsolve(ode1, U1(y, n));

Удобно записать общее решение ОДУ так

>V1:=proc(y) options operator, arrow;

C1*sinh(n*Pi*y)+C2*sinh(n*Pi*(2-y))+

(-1+cos(1)*(-1)^n)*y/(-2*n*Pi+2*n^3*Pi^3)

end proc;

Проверим

>subs(U1(y,n) = V1(y), ode1);

simplify(subs(U1(y,n) = V1(y), ode1));

Все в порядке!

Формируем систему уравнений по граничным условиям

>subs({U1(0,n) = V1(0), U1(2,n) = V1(2)}, bcODE1);

Решаем систему

>solve(%, {C1, C2}); assign(%):

Проверка:

>subs(U1(y,n) = V1(y), ode1): simplify(%);

Итак, трансформанта найдена:

>U1:= unapply(V1(y), y, n);

Решение задачи 1 дается следующим рядом

>u1:=proc(x,y) options operator, arrow;

2*(Sum(U1(y, n)*X(x, n), n = 1.. infinity))

end proc;

> u1(x,y);

Проверим найденное решение. Подставим решение в уравнение; отдельно рассмотрим левую и правую части уравнения. Левая часть:

> LHSpde:= combine(lhs(pde1)); RHSpde:= rhs(pde1);

Упростим левую часть уравнения:

> op(1, LHSpde);

> numer(op(1, LHSpde));

> expand(%);

> factor(%);

Таким образом, левая часть ДУЧП имеет вид

>LHSpde:=Sum(-2*Pi*n*y*sin(n*Pi*x)*(-1+cos(1)*

(-1)^n)/(-2+2*n^2*Pi^2), n=1..infinity);

Правая часть уравнения:

> RHSpde:= rhs(pde1);

Покажем, что левая часть уравнения совпадает с правой частью, т. е.

>1/2*y*cos(x)=Sum(-2*Pi*n*y*sin(n*Pi*x)*(-1+cos(1)*

(-1)^n)/(-2+2*n^2*Pi^2),n = 1.. infinity);

Для этого разложим функцию в ряд по собственным функциям на отрезке . Коэффициенты разложения:

>2*int(1/2*y*cos(x)*sin(n*Pi*x), x = 0.. 1);

simplify(%, assume = integer);

Таким образом, уравнение выполняется.

Проверим выполнение граничных условий.

Граничные условия по переменной :

> simplify(u1(0,y), assume = integer);

simplify(u1(1,y), assume = integer);

Граничные условия по переменной :

> u1(x,0);

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

> 2*int(f(x,0)*sin(n*Pi*x), x = 0.. 1);

simplify(%, assume = integer);

Далее,

> u1(x,2);

Это — разложение граничной функции в ряд по функциям на отрезке . Действительно, преобразуем общий член ряда

> q:=op(2, u1(x,2));

> q:=op(1,q);

> op(1,q);

> q:=%:combine(q);

Таким образом, мы имеем ряд

> 'u1(x,2)'=2*Sum(%*sin(n*Pi*x),n=1..infinity);

Коэффициенты разложения граничной функции в ряд по функциям на отрезке :

>2*int(f(x,2)*sin(n*Pi*x), x = 0.. 1);

simplify(%, assume = integer);

что и требовалось доказать.

Итак, решение задачи 1 найдено. Читателю предлагается построить решение задачи 2 самостоятельно, в качестве упражнения. Приведем формулу этого решения, полученную в Maple