Файл: Учебное пособие СанктПетербург 2016.pdf

Скачать файл
Заказать решение

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 30.11.2023

Просмотров: 139

Скачиваний: 2

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

СОДЕРЖАНИЕ

= Jp(√λx). Подставив решение в краевое условие и учитывая, что dJp(√λx)dx=√λJ0p(√λx),получимR√λJ0p(√λT ) + SJp(√λT ) = 0.(1.26)Обозначим√λT =γи перепишем уравнение (1.26) в видеRγJ0p(γ) + ST Jp(γ) = 0.(1.27)Как известно из 1.5, это уравнение имеет множество простых корнейγk, k = 1, 2, ... . В случае, если R = 0 или S = 0, уравнение (1.27) превраща- ется в Jp(γ) = 0 или J0p(γ) = 0 с аналогичными свойствами. Следовательно,λ будет собственным числом задачи, если√λT =γk. Таким образом, по- лучено множество собственных чиселλk=γkT2,гдеγk– корень уравнения (1.27), k = 1, 2, ..., и множество собственных функций yk(x) = JpγkTxСистема собственных функций nJpγkTxo+∞k=1образует полную ортого- нальную систему в пространстве L2[0, T ; x].Для того чтобы разложить функцию в ряд по системе ортогональных функций Бесселя, потребуется несколько интегральных тождеств.26 Утверждение 1.5. Для любыхα,β∈ R,α6=β, выполняетсяTZ0JpαTxJpβTxx dx = T2βJp(α)J0p(β) −αJ0p(α)Jp(β)α2−β2(1.28)Доказательство. Заметим, что функция y(x) = JpαTxявляется решением уравнения с λ =α2T2, т. е. справедливо равенство1xxJpαTx00+α2T2−p2x2JpαTx= 0.Аналогично, если λ =β2T2, то1xxJpβTx00+β2T2−p2x2JpβTx= 0.Умножим первое равенство на JpβTxx, а второе – на JpαTxx и вы- чтем из первого второе:xJpαTx00JpβTx−xJpβTx00JpαTx++α2−β2T2JpαTxJpβTxx = 0.Проинтегрировав полученное выражение, получимα2−β2T2TZ0JpαTxJpβTxx dx ==TZ0xJpβTx00JpαTxdx −TZ0xJpαTx00JpβTxdx.Интегралы в правой части выражения возьмем по частям:α2−β2T2TZ0JpαTxJpβTxx dx =27 = xJpβTx0JpαTxT0−TZ0xJpβTx0JpαTx0dx −− xJpαTx0JpβTx T0−TZ0xJpαxT0JpβTx0dx ==βJ0p(β)Jp(α) −αJ0p(α)Jp(β).Поделив наα2−β2T2, получим формулу (1.28).Утверждение 1.6. Для любогоα6= 0 справедлива формулаTZ0J2pαTxx dx =T2 2J0p(α)2+1 −p2α2J2p(α)(1.29)Доказательство. В выражении (1.28) сделаем предельный переход приβ→αи применим правило Лопиталя:TZ0J2pαTxx dx = T2limβ→αβJp(α)J0p(β) −αJ0p(α)Jp(β)α2−β2== T2limβ→αJp(α)J0p(β) +βJp(α)J00p(β) −αJ0p(α)J0p(β)−2βИз уравнения Бесселя (1.19) следует, что J00p(β) = −1βJ0p(β)−1 −p2β2Jp(β),поэтомуTZ0J2pαTxdx == T2limβ→αJp(α)J0p(β) − Jp(α)J0p(β) −β1 −p2β2Jp(α)Jp(β) −αJ0p(α)J0p(β)−2β==T2 2J0p(α)2+1 −p2α2J2p(α)28 1.7. Оператор Лежандра. Многочлены ЛежандраОператором Лежандра называется дифференциальный оператор видаL(y) = − (1 − x2)y00За область определения оператора возьмем функции, дважды дифферен- цируемые на (−1; 1). Так как на концах интервала функция p(x) = 1 − x2обращается в нуль, то в качестве однородных краевых условий потребуем ограниченность y(x) при x → ±1.Поставим для оператора Лежандра задачу на собственные значения(− (1 − x2)y00= λy,y(x) − ограничена при x → ±1.(1.30)Уравнение задачи можно переписать в виде(x2− 1)y00+ 2xy0− λy = 0.(1.31)Рассмотрим теперь выражения видаPn(x) =1 2n n!(x2− 1)n(n)Достаточно очевидно, что Pn(x) является многочленом степени n,P0(x) ≡ 1, P1(x) =1 2(x2− 1)0= x, P2(x) =1 4 · 2(x2− 1)200=1 2(3x2− 1).Функции Pn(x) называются многочленами Лежандра.Утверждение 1.7. Функция Pn(x) является собственной функци- ей оператора Лежандра, соответствующей собственному числу λn== n(n + 1).Доказательство. При доказательстве используется формула диффе- ренциального бинома(uv)(n)= uv(n)+ C1n u0v(n−1)+ C2n u00v(n−2)+ ... + u(n)v.(1.32)Здесь Ck n=n(n − 1)...(n − k + 1)k!– биномиальные коэффициенты.Запишем очевидное равенство(x2− 1)n0= n(x2− 1)n−1 2x.Умножим левую и правую части этого равенства на x2− 1 и n + 1 раз продифференцируем:h(x2− 1) (x2− 1)n0i(n+1)=n(x2− 1)n2x(n+1)29 Последнее равенство распишем по дифференциальному биному, учитывая,что (m + 1)-я производная от многочлена степени m равна нулю:(x2−1)(x2− 1)n(n+2)+(n+1)2x(x2− 1)n(n+1)+(n + 1)n2!2(x2− 1)n(n)== 2nx(x2− 1)n(n+1)+ 2n(n + 1)(x2− 1)n(n),или после приведения подобных(x2− 1)(x2− 1)n(n+2)+ 2x(x2− 1)n(n+1)− n(n + 1)(x2− 1)n(n)= 0.Умножив последнее соотношение на1 2n n!, получим(x2− 1)P00n(x) + 2xP0n(x) − n(n + 1)Pn(x) = 0.Отсюда видно, что Pn(x) удовлетворяет уравнению (1.31) с λ == n(n + 1).Можно доказать [2], что при других значениях λ задача (1.30) имеет только тривиальное решение y ≡ 0.Таким образом, {λn= n(n + 1)}+∞n=0образуют систему собственных чи- сел оператора Лежандра, а многочлены {Pn(x)}+∞n=0– систему собственных функций. Система {Pn(x)}+∞n=0полна в пространстве L2[−1; 1].Перед тем как вычислять норму Pn(x), докажем вспомогательное утвер- ждение.Утверждение 1.8. Многочлен Q(x) =(x2− 1)n(m), где m < n, об- ращается в нуль на концах промежутка [−1; 1].Доказательство. Распишем Q(x), используя формулу дифференци- ального бинома (1.32):Q(x) = [(x − 1)n(x + 1)n](m)== ((x − 1)n)(m)(x+1)n+C1n((x − 1)n)(m−1)((x+1)n)0+...+(x−1)n((x + 1)n)(m)Так как m < n, все слагаемые этого выражения содержат множители (x−1)и (x + 1) и, следовательно, Q(−1) = Q(1) = 0.Найдем теперь норму собственной функции Pn(x):kPn k2=1Z−1P2n(x)dx =1(2n n!)2 1Z−1(x2− 1)n(n)(x2− 1)n(n)dx.Применим формулу интегрирования по частям:kPn k2=1(2n n!)(x2− 1)n(n)(x2− 1)n(n−1)1−1−30 −1(2n n!)2 1Z−1(x2− 1)n(n+1)(x2− 1)n(n−1)dx.Согласно утверждению 1.8 первое слагаемое равно нулю, а интеграл опять преобразуем по частям и т. д., пока не получим kPn k2=1(2n n!)2(−1)n1Z−1(x2− 1)n(2n)(x2− 1)n dx ==(−1)n(2n)!(2n n!)2 1Z−1(x2− 1)n dx.Последний интеграл снова возьмем по частям:1Z−1(x2− 1)n dx =1Z−1(x − 1)n(x + 1)n dx == (x − 1)n(x + 1)n+1n + 1 1−1−n n + 1 1Z−1(x − 1)n−1(x + 1)n+1dx ==(−1)n n!(n + 1)(n + 2)...2n1Z−1(x + 1)2n dx =(−1)n n!n!(2n)!2 2n+1 2n + 1В итоге получим kPn k2=(−1)n(2n)!(2n n!)2(−1)n(n!)2(2n)!2 2n+1(2n + 1)=2 2n + 1(1.33)Таким образом, любую функцию f ∈ L2[−1; 1] можно разложить в рядФурье по многочленам Лежандра f (x) =∞Xn=0c nPn(x),где коэффициент Фурье вычисляется по формуле cn=(f, Pn)kPn k2=2n + 1 21Z−1f (x)Pn(x)dx.Отметим еще несколько свойств многочленов Лежандра, которые яв- ляются общими для любых систем ортогональных многочленов.31 Утверждение 1.9. Пусть Q(x) – произвольный многочлен степени n, тогда его можно представить в виде Qn(x) =nXk=0αkPk(x), при этомαk=(Qn, Pk)kPk k2Доказательство. Этот факт очевидным образом следует из того, что система P0(x), ..., Pn(x) ортогональна и, следовательно, линейно независи- ма, а потому ее можно принять за базис линейного пространства много- членов степени не выше n. Поскольку многочлены ортогональны, то коэф- фициенты разложенияαk вычисляются по формуламαk=(Qn, Pk)kPk k2Утверждение 1.10. Пусть Pn(x) – многочлен Лежандра, Qm(x) –произвольный многочлен степени m < n. Тогда Pn и Qm ортогональны.Доказательство. Запишем Qm(x) =mXk=0αkPk(x). Скалярное произ- ведение(Qm, Pn) =mXk=0αkPk, Pm!=mXk=0αk(Pk, Pn) = 0.Утверждение 1.11. Все корни многочлена Лежандра простые и ле- жат в интервале (−1; 1).Доказательство. Обозначим через x1, x2, ..., x m∈ (−1; 1) точки, где многочлен Pn(x) меняет знак, и предположим, что m < n. Тогда много- член Qm(x) = (x − x1)(x − x2)...(x − x m) меняет знак в тех же точках и произведение Qm(x)Pn(x) на промежутке [−1; 1] знак не меняет. А тогда(Qm, Pn) =1Z−1Qm(x)Pn(x)dx 6= 0,что противоречит утверждению 1.8 .Для многочленов Лежандра выполняются рекуррентные соотноше- ния. Приведем без доказательства две рекуррентные формулы:(n + 1)Pn+1(x) − (2n + 1)xPn(x) + Pn−1(x) = 0,(1.34)(2n + 1)Pn(x) = P0n+1(x) − P0n−1(x).(1.35)Вывод этих формул можно найти, например, в [7] или [8].32 1.8. Присоединенные функции ЛежандраПрисоединенной функцией Лежандра степени n порядка k называется функция видаPk n(x) = (1 − x2)k/2P(k)n(x),k = 0, 1, ...,n = k, k + 1, ...,где Pn(x) – многочлен Лежандра. Заметим, что для четных k функцияPk n(x) будет многочленом степени n, а для нечетных – многочленом степени n − 1, умноженным на√1 − x2. Например,P0 0(x) = 1,P0 1(x) = P1(x) = x,P1 1(x) = (1 − x2)1/2P0 1(x) =√1 − x2,P0 2(x) = P2(x) =3 2x2−1 2,P1 2(x) = (1 − x2)1/2P0 2(x) =√1 − x2 3x,P2 2(x) = (1 − x2)P00 2(x) = (1 − x2)3,P2 3(x) = (1 − x2)P00 3(x) = (1 − x2)15x и т. д.Перед тем как исследовать свойства присоединенных функций Ле- жандра, получим уравнение для производных многочленов Лежандра.Утверждение 1.12. Производные многочленов Лежандра P(k)n(x) удо- влетворяют уравнению−(1 − x2)z00+ 2(k + 1)xz0+ k(k + 1)z = n(n + 1)z.(1.36)Доказательство. Многочлен Лежандра удовлетворяет уравнению(1.31)(x2− 1)P00n(x) + 2xP0n(x) = n(n + 1)Pn(x).Продифференцируем данное соотношение k раз:(x2− 1)P00n(x)(k)+ (2xP0n(x))(k)= n(n + 1)P(k)n(x)и воспользуемся формулой дифференциального бинома (1.32):(x2−1)P(k+2)n(x) + k2xP(k+1)n(x) +k(k − 1)2 2P(k)n(x) + 2xP(k+1)n(x) + k2P(k)n== n(n + 1)P(k)nПосле приведения подобных получаем(x2− 1)(P(k)n)00+ 2(k + 1)x(P(k)n)0+ k(k + 1)P(k)n= n(n + 1)P(k)nУтверждение 1.13. Присоединенная функция Лежандра Pk n(x) удо- влетворяет уравнению−(1 − x2)y00+ 2xy0+k2 1 − x2y = n(n + 1)y.(1.37)33 Доказательство. Сделаем в уравнении замену y = (1 − x2)k/2z. Най- дем из этого равенства y0и y00и, подставив в (1.37), после преобразований получим−(1 − x2)z00+ 2(k + 1)xz0+ k(k + 1)z = n(n + 1)z,т. е. уравнение (1.36).Уравнение (1.37) можно записать в симметричной форме− (1 − x2)y00+k2 1 − x2y = n(n + 1)y.(1.38)Если к уравнению (1.38) добавить однородные краевые условия: y(x) огра- ничена при x → ±1, то получим задачу на собственные значения для опе- ратора Lk(y) = − (1 − x2)y00+k2 1 − x2y. Таким образом, для оператораLk(y) найден набор собственных чисел {λn= n(n + 1)}+∞n=k и набор соб- ственных функций {Pk n(x)}+∞n=k. Присоединенные функции Лежандра Pk nиPk m(n 6= m) ортогональны, так как они являются собственными функция- ми симметричного оператора.Утверждение 1.14. Для присоединенных функций Лежандра спра- ведлива формула(Pk m, Pk n) = (n − k + 1)(n + k)(Pk−1m, Pk−1n).(1.39)Доказательство. Найдем скалярное произведение(Pk m, Pk n) =1Z−1P(k)m(x)P(k)n(x)(1 − x2)k dx.Данный интеграл возьмем по частям, считая, что v0(x) = P(k)m(x), а u(x) == P(k)n(x)(1 − x2)k:(Pk n, Pk m) = P(k−1)m(x)P(k)n(x)(1 − x2)k1−1−−1Z−1hP(k+1)n(x)(1 − x2)k− 2kxP(k)n(x)(1 − x2)k−1iP(k−1)m(x)dx.Внеинтегральное слагаемое очевидно равно нулю, поэтому(Pk n, Pk m) =1Z−1h−(1 − x2)P(k+1)n+ 2kxP(k)n iP(k−1)m(1 − x2)k−1dx.(1.40)34 Согласно формуле (1.36) P(k−1)n будет удовлетворять уравнению−(1 − x2)P(k−1)n(x)00+ 2kxP(k−1)n(x)0+ (k − 1)kP(k−1)n(x) == n(n + 1)P(k−1)n(x).Отсюда−(1 − x2)P(k+1)n(x) + 2kxP(k)n(x) = (n2+ n − k2+ k)P(k−1)n(x).Подставив полученное соотношение в (1.40), получаем(Pk n, Pk m) = (n2+ n − k2+ k)1Z−1P(k−1)n(x)P(k−1)m(x)(1 − x2)k−1dx == (n − k + 1)(n + k)(Pk−1n, Pk−1m).Заметим, что в формуле (1.39) коэффициент не зависит от m. Это совершенно естественно, если вспомнить, что при m 6= n скалярное произ- ведение (Pk m, Pk n) = 0. При n = m имеем kPk nk2= (Pk n, Pk n) = (n − k + 1)(n + k)(Pk−1n, Pk−1n) == (n − k + 1)(n + k)(n − k + 2)(n + k − 1)(Pk−2n, Pk−2n) = ... == (n + k)(n + k − 1)...(n + 1)(n − k + 1)(n − k + 2)...n(P0n, P0n).Учитывая, что P0n(x) = Pn(x), и формулу (1.33), получим kPk nk2=(n + k)!(n − k)!2 2n + 1Для фиксированного значения k набор {Pk n}+∞n=k образует полную орто- гональную систему в пространстве L2[−1; 1]. Любую функцию f ∈ L2[−1; 1]можно разложить в ряд Фурье по присоединенным функциям Лежандра f (x) =∞Xn=k ckPk n(x),где c k– коэффициенты Фурье:c k=(f, Pk n)kPk nk2=(n − k)!(n + k)!(2n + 1)2 1Z−1f (x)Pk n(x)dx.35 1.9. Решение краевой задачи методом конечных разностейВ предыдущих параграфах было получено решение краевой задачи в виде аналитического выражения или в виде ряда. Теперь рассмотрим численный метод решения краевых задач, который называется методом конечных разностей, или методом сеток.Пусть отрезок [a; b] разбит на n частей. Введем шаг сетки h =b − a n,точки x0= a, x1= x0+ h, x2= x0+ 2h, ..., x n= x0+ nh = b назовем узлами сетки, а само множество узловωh= {x k} (k = 0, ..., n) – сеткой. Если на отрезке задана непрерывная функция y(x), то функцию {y k= y(x k)}n k=0естественно называть сеточным аналогом функции y(x). Для аппроксима- ции производных функции в узлах сетки будем рассматривать разностные отношения (линейные комбинации значений сеточной функции в несколь- ких узлах сетки). Из геометрических соображений понятно, что производ- ную функции в узле x kможно аппроксимировать следующим образом:y0(x k) ≈y k+1− y kh или y0(x k) ≈y k− y k−1hПервое из этих выражений принято называть разностным отношени- ем “вперед”, и его можно использовать в узлах x0, x1, ..., x n−1, второе –разностное отношение “назад” подходит для узлов x1, x2, ..., x n. Найдем связь между производными дифференцируемой функции y(x) и ее сеточ- ным аналогом.Утверждение 1.15. Пусть y ∈ C2[a,b], тогда справедливы оценки y0(x k) −y(x k+1) − y(x k)h≤M2h2,y0(x k) −y(x k) − y(x k−1)h≤M2h2,где M2= max x∈[a,b]|y00(x)|.Доказательство. В узле x kзапишем формулу Тейлора первого по- рядка с остатком в форме Лагранжа y(x) = y(x k) + y0(x k)(x − x k) +1 2y00(˜x)(x − x k)2,где ˜x ∈ [x k, x k+1]. Положим x = x k+1, тогда x k+1− x k= h и y(x k+1) = y(x k) + y0(x k)h +1 2y00(˜x)h2,36 откуда y(x k+1) − y(x k)h− y0(x k) =1 2y00(˜x)h.(1.41)Заменив y00(˜x) на максимум второй производной, получим оценку y0(x k) −y(x k+1) − y(x k)h≤M2h2Для разностного отношения “назад” доказательство аналогично.Отметим принципиальный момент. Для того чтобы аппроксимировать первую производную функции разностным отношением, нужно, чтобы функ- ция y(x) была дважды дифференцируемой. Для аппроксимации первой производной y0(x k) можно использовать и другие выражения, например симметричную разность (y k+1− y k−1)/2h.Утверждение 1.16. Пусть y ∈ C3[a,b], тогда справедлива оценка y0(x k) −y(x k+1) − y(x k−1)2h≤M3h2 6,где M3= max x∈[a,b]|y000(x)|, k = 1, 2, ..., n − 1.Доказательство. Запишем формулу Тейлора 2-го порядка с остатком в форме Лагранжа y(x) = y(x k) + y0(x k)(x − x k) +1 2y00(x k)(x − x k)2+1 6y000(˜x)(x − x k)3Для точек x k+1и x k−1получим:y(x k+1) = y(x k) + y0(x k)h +1 2y00(x k)h2+1 6y000(˜x)h3,y(x k−1) = y(x k) − y0(x k)h +1 2y00(x k)h2−1 6y000(≈x)h3,где ˜x ∈ [x k, x k+1],≈x∈ [x k−1, x k]. Вычтем из первой формулы вторую:y(x k+1) − y(x k−1) = 2y0(x k)h +1 6y000(˜x) + y000(≈x)h3,откуда y(x k+1) − y(x k−1)2h− y0(x k) =1 12y000(˜x) + y000(≈x)h2(1.42)Далее, заменив y000(˜x) и y000(≈x)) на M3, получим искомую оценку.Для аппроксимации второй производной обычно используется выра- жение y00(x k) ≈y k−1− 2y k− y k+1h2 37 Утверждение 1.17. Пусть y ∈ C4[a,b], тогда справедлива оценка y00(x k) −y(x k−1) − 2y(x k) + y(x k+1)h2≤M4h2 12,где M4= max x∈[a,b]|y(4)(x)|.Доказательство. Запишем формулу Тейлора 3-го порядка с остатком в форме Лагранжа y(x) = y(x k) + y0(x k)(x − x k) +1 2y00(x k)(x − x k)2+1 6y000(x k)(x − x k)3++1 24y(4)(˜x)(x − x k)4Подставив x k+1и x k−1вместо x, соответственно, получим:y(x k+1) = y(x k) + y0(x k)h +1 2y00(x k)h2+1 6y000(x k)h3+1 24y(4)(˜x)h4,y(x k−1) = y(x k) − y0(x k)h +1 2y00(x k)h2−1 6y000(x k)h3+1 24y(4)(≈x)h4,где ˜x ∈ [x k, x k+1],≈x∈ [x k−1, x k]. Сложим полученные формулы:y(x k+1) + y(x k−1) = 2y(x k) + y00(x k)h2+1 24y(4)(˜x) + y(4)(≈x)h4Отсюда y(x k−1) − 2y(x k) + y(x k+1)h2− y00(x k) =1 24y(4)(˜x) + y(4)(≈x)h2(1.43)Заменив в полученном выражении значения y(4)(˜x) и y(4)(≈x) на их макси- мум M4, получим нужную оценку.Из полученных равенств (1.41)–(1.43) следует, что соответствующее разностное отношение аппроксимирует производную в узле x kс погреш- ностью, которая ведет себя как произведение ограниченной функции на hm(m зависит от вида разностного отношения). Для описания таких по- грешностей удобно использовать символ O(h m). Опираясь на определение,известное из курса математического анализа, поясним смысл этого симво- ла. Символом O(h m) в точке x kобозначается любая функция, которая в окрестности точки x kведет себя как произведение ограниченной функции на h m, где h = x k+1− x kили h = x k− x k−1. Использование этого сим- вола позволит однотипно описывать погрешности аппроксимации во всех точках сеткиωh. Кроме того, при решении сеточных задач важно видеть,38 как зависит погрешность аппроксимации от шага сетки. Из полученных оценок ясно, что эта погрешность зависит только от h (множитель при h mс уменьшением шага существенно не влияет на изменение погрешности).Таким образом, в каждой внутренней точке сеткиωh справедливы равенства:y0k(x k) =y k+1− y kh+ O(h),(1.44)y0k(x k) =y k− y k−1h+ O(h),(1.45)y0k(x k) =y k+1− y k−1 2h+ O(h2),y00(x k) =y k−1− 2y k+ y k+1h2+ O(h2).(1.46)Основные идеи метода конечных разностей рассмотрим на примере краевой задачи−y00+ q(x)y = f (x),y(a) = 0,y(b) = 0,(1.47)где q(x) ≥ q0> 0. Введем на [a, b] сетку {x0, x1, ..., x n} с шагом h. Во внут- ренних узлах сетки x kзапишем уравнение, заменяя производную y00(x k)соответствующим разностным отношением (1.46), в граничных узлах сет- ки x0и x n– краевые условия:−y k−1− 2y k+ y k+1h2+ O(h2) + q(x k)y k= f (x k), k = 1, ..., n − 1,(1.48)y0= 0,y n= 0.Получилась система уравнений относительно y k= y(x k), k = 0, ..., n. Пред- полагая, что погрешности O(h2) малы, отбросим их. Тогда придем к систе- ме уравнений относительно сеточной функции ˜y−˜y k−1− 2˜y k+ ˜y kh2+ q(x k)˜y k= f (x k),k = 1, ..., n − 1,(1.49)˜y0= 0,˜y n= 0.Эта система называется разностной схемой для краевой задачи (1.47).Числа˜y k– это приближенные значения функции y в узлах сетки.Учитывая граничные условия, систему можно записать только для внутренних узлов сетки. Преобразуем k-е уравнение к виду−˜y k−1+ 2 + h2q(x k) ˜y k− ˜y k+1= h2f (x k).39 Последовательно подставим k = 1, 2, ..., n − 1. С учетом того, что ˜y0= ˜y n== 0, получим(2 + q(x1)h2)˜y1− ˜y2= f (x1)h2,−˜y1+ (2 + q(x2)h2)˜y2− ˜y3= f (x2)h2,−˜y2+ (2 + q(x3)h2)˜y3− ˜y4= f (x3)h2,−˜y n−2+ (2 + q(x n−1)h2)˜y n−1= f (x n−1)h2(1.50)Обозначим˜y =˜y1˜y2˜y n−1;f =f (x1)f (x2)f (x n−1);A =2 + q(x1)h2−1 00 0−1 2 + q(x2)h2−1 ...0 00 00... −1 2 + q(x n−1)h2и запишем систему (1.50) в матричном видеA˜y = f h2Матрица A – трехдиагональная с диагональным преобладанием (модуль диагонального элемента больше, чем сумма модулей остальных элементов строки или столбца). Справедливо следующее утверждение.Утверждение 1.18. Если A = [a ij], i = 1, ..., n, j = 1, ..., n – матри- ца со строгим диагональным преобладанием(|a ij|>Xj6=i|a ij|,i = 1, ..., n), тогда эта матрица невырожденна.Доказательство. Предположим, что матрица A вырожденна, тогда существует вектор x 6= 0, для которого Ax = 0. Пусть x p– максимальный по модулю элемент вектора x. Очевидно, что |x p| > 0. Из уравнения Ax = 0следует, что|a pp xp| =Xj6=p apj xj≤ |x p|Xj6=p|a pj|.Учитывая, что |a pp xp| = |a pp| |x p|, получим|a pp| ≤Xj6=p|a pj|.Это неравенство противоречит тому, что A – матрица со строгим диаго- нальным преобладанием. Значит, A – невырожденная матрица.40 Из этого утверждения следует, что система уравнений (1.50) однознач- но разрешима.Покажем теперь, что с уменьшением шага сетки h значения реше- ния системы (1.50) будут стремиться к значениям решения краевой задачи(1.47) в узлах сетки. Для этого введем понятие сеточной нормы.Сеточной нормой функции v, определенной в узлах сеткиωh= {x k},k = 0, ..., n, назовем наибольшее по модулю значение этой функции на сетке kvk h= max k=0,...,n|v k|.Решение ˜y разностной краевой задачи при измельчении сетки сходится к решению y краевой задачи (1.47), если ky − ˜yk h→ 0при h → 0.Если, кроме того, выполняется неравенство ky − ˜yk h≤ Ch k,где C не зависит от h, тогда разностная схема имеет k-й порядок точности.Сходимость разностной схемы связана с двумя понятиями: устойчиво- стью и точностью аппроксимации.Разностная схема называется устойчивой по правой части, если норма ее решения не превосходит нормы функции f , заданной в задаче, умножен- ной на число, не зависящее от шага сетки:k˜yk h≤ Ckf k h(1.51)Утверждение 1.19. Разностная схема (1.49) устойчива по правой части, т. е. справедлива оценка (1.51), где C не зависит от шага h.Доказательство. Пусть наибольшее среди чисел |˜y k|, k = 0, ..., n,есть число |˜y p|. Если p = 0 или p = n, то неравенство (1.51) выполняется(˜y0= 0, ˜y n= 0). Предположим, что 0 < p < n. Для всех k = 0, ..., n|˜y k| ≤ |˜y p|. Рассмотрим разностную схему (1.49) для узла сетки x p:(2 + h2q(x p))˜y p= ˜y p−1+ ˜y p+ h2f (x p).Так как 2 + h2q(x p) > 0, то(2 + h2q(x p)) |˜y p| = |˜y p−1+ ˜y p+ h2f (x p)| ≤ |˜y p−1| + |˜y p| + h2|f (x p)|.Учитывая, что |˜y k| ≤ |˜y p| для всех k = 0, ..., n, получим(2 + h2q(x p))|˜y p| ≤ 2|˜y p| + h2|f (x p)|,или|˜y p| ≤|f (x p)|q(x p)41 Из этого равенства следует оценка k˜yk h≤1q0kf k h,означающая устойчивость разностной схемы (1.49).Подставим в разностные уравнения (1.49) вместо значений сеточной функции ˜y kзначения точного решения дифференциальной задачи y(x k).Для сохранения равенства в правую часть (1.49) необходимо ввести допол- нительное слагаемое. Если это слагаемое с уменьшением h изменяется какO(h m), то говорят, что разностная схема аппроксимирует дифференциаль- ную задачу с погрешностью O(h m).Утверждение 1.20. Разностная схема (1.49) аппроксимирует кра- евую задачу (1.47) с погрешностью O(h2).Доказательство. Это утверждение следует из равенств (1.48). Во внутренних узлах сетки функция y удовлетворяет разностной схеме с по- грешностью O(h2). Граничные условия задаются точно.Из устойчивости по правой части и аппроксимации с погрешностьюO(h2) следует сходимость разностной схемы.Утверждение 1.21. Разностная схема (1.49) сходится к решению задачи (1.47) в сеточной норме при h → 0.Доказательство. Обозначим w k= y k− ˜y k. Из уравнений (1.48), (1.49)относительно функций y и ˜y следует, что сеточная функция w удовлетво- ряет разностной схеме−w k−1− 2w k+ w k+1h2+ q(x k)w k= O(h2),k = 1, ..., n − 1,(1.52)w0= 0,w n= 0.Таким образом, разность между точным и приближенными решениями краевой задачи удовлетворяет таким же сеточным уравнениям, что и функ- ция ˜y, только в правой части уравнений стоит погрешность аппроксимации второй производной.Согласно утверждению 1.19 об устойчивости разностная схема (1.52)устойчива и справедлива оценка kwk h≤ C kO(h2)k h,где C не зависит от h.Из этого следует, что kwk h→ 0 при h → 0. Так как справедлива оцен- ка |O(h2)| ≤ M h2, то построенная разностная схема (1.49) имеет второй порядок точности.Если краевые условия в задаче (1.47) более сложные, например:y0(a) − S1y(a) = t1,(1.53)42 y0(b) + S2y(b) = t2,(1.54)то их тоже надо аппроксимировать разностными уравнениями.Так, для условия (1.53) используют разностное отношение (1.44)y1− y0h+ O(h) − S1y0= t1,а для (1.54) – соотношение (1.45)y n− y n−1h+ O(h) + S2y n= t2Тогда, если отбросить погрешности, получается(1 + S1h)˜y0− ˜y1= −t1h,−˜y n−1+ (1 + S2h)˜y n= t2h.Эти уравнения нужно добавить к системе (1.50). Количество уравнений увеличится, однако матрица системы по-прежнему будет трехдиагональ- ной с диагональным преобладанием. Полученная разностная схема будет аппроксимировать краевую задачу с погрешностью O(h).2. УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ2.1. Вывод одномерного уравнения теплопроводностиРассмотрим однородный стержень длины l, боковая поверхность кото- рого теплоизолирована. Стержень предполагается настолько тонким, что в каждый момент времени температуру в любом его поперечном сечении можно считать постоянной. Если этот стержень неравномерно нагрет, то- гда в нем будет происходить процесс передачи тепла от более нагретых участков к менее нагретым. Для вывода уравнения теплопроводности бу- дем использовать следующую модель: ось 0x направим вдоль оси стержня,совместив начало координат с его левым концом. Тогда процесс распро- странения тепла в стержне может быть описан функцией u(x, t), представ- ляющей температуру в сечении стержня x в момент времени t. Частная производная∂u(x, t)∂x выражает скорость изменения температуры в направ- лении оси 0x. Если температура в направлении оси 0x растет, тогда∂u∂x> 0,если температура уменьшается, то∂u∂x< 0.Вывод дифференциального уравнения теплопроводности базируется на следующих экспериментальных положениях:43 1. Количество тепла ∆Q, которое необходимо сообщить однородному телу, чтобы повысить его температуру на ∆u, равно:∆Q = cρV ∆u,(2.1)где c – удельная теплоемкость;ρ– плотность; V – объем тела.2. Количество тепла q(x), протекающее через поперечное сечение x стержня за промежуток времени [t, t + ∆t], пропорционально скорости из- менения температуры в направлении, перпендикулярном к сечению, т. е.∂u∂x, площади сечения S и времени ∆t:q(x) = −k∂u∂xS ∆t,(2.2)где k – коэффициент теплопроводности.Величину −k∂u∂x называют удельным тепловым потоком. Знак минус объясняется тем, что величина теплового потока считается положитель- ной, когда тепло передается в сторону возрастания x. Если∂u∂x> 0, то в направлении оси 0x температура увеличивается, а так как тепло перехо- дит от более нагретых участков к менее нагретым, то тепловой поток будет направлен в сторону уменьшения x.Рассмотрим малый элемент стержня [x; x + ∆x]. Выведем уравнение,которому будет удовлетворять функция u(x, t). Для этого составим уравне- ние теплового баланса для выделенного элемента стержня за промежуток времени ∆t.Количество тепла q(x), входящее через поперечное сечение с абсцис- сой x за промежуток времени ∆t, определяется соотношением (2.2). Если воспользоваться формулой Тейлора f (x + ∆x) = f (x) + f0(x)∆x + o(∆x),величина теплового потока, выходящего через сечение x + ∆x, будет сле- дующей:q(x + ∆x) = −k∂u∂xS ∆t −∂∂xk∂u∂xS ∆t∆x + o(∆x).Считаем приращение ∆x настолько малым, что величиной o(∆x) можно пренебречь.Если найти разность значений входящего и выходящего тепловых по- токов на участке [x; x + ∆x] за время ∆t, получится количество тепла ∆Q,сообщенное выбранному участку за указанное время:∆Q = −k∂u∂xS ∆t +k∂u∂xS ∆t +∂∂xk∂u∂xS∆x∆t=44 =∂∂xk∂u∂xS∆x∆t.При этом за тот же промежуток времени температура в каждой точке выделенного участка изменилась на величину∆u = u(x, t + ∆t) − u(x, t) =∂u∂t∆t + o(∆t).Считаем,что ∆t мало. Пренебрегая величиной o(∆t), используя соотноше- ние (2.1) и учитывая, что V = S∆x, получим∆Q = cρS∆x∂u∂t∆t.Приравнивая полученные выражения для ∆Q, составим уравнение тепло- вого баланса:ρcS∆x∆t∂u∂t=∂∂xk∂u∂xS∆x∆t.(2.3)После сокращения на общий множитель S∆x∆t получим уравнение cρ∂u∂t=∂∂xk∂u∂x,которое называется линейным уравнением теплопроводности без тепловых источников. Если k, c,ρ– постоянные величины, то полученное уравнение можно записать в виде∂u∂t= a2∂2u∂x2(2.4)Число a2=k cρназывается коэффициентом температуропроводности.Предположим теперь, что в некоторых частях стержня находятся ис- точники тепла (выделение или поглощение тепла может происходить в ре- зультате прохождения электрического тока, вследствие химической реак- ции и т. п.). Пусть g(x, t) – функция, описывающая плотность тепловых ис- точников, т. е. такая функция, что на малом участке стержня [x; x + ∆x] за малый промежуток времени ∆t выделяется или поглощается тепло, равное g(x, t)S∆x∆t. Если тепло выделяется, тогда g(x, t) > 0, если поглощается,то g(x, t) < 0.Если в уравнении (2.3) учесть тепло, выделяемое (поглощаемое) ис- точниками на участке стержня [x; x + ∆x] за время ∆t, тогда уравнение теплового баланса после сокращения на S∆x∆t примет видρc∂u∂t=∂∂xk∂u∂x+ g(x, t).45 Пусть f (x, t) =1ρc g(x, t), тогда в случае постоянных k, c,ρполученное уравнение можно записать в виде∂u∂t= a2∂2u∂x2+ f (x, t).(2.5)Это линейное уравнение теплопроводности с учетом источников тепла.Уравнение теплопроводности (2.4) является линейным однородным, а уравнение (2.5) – линейным неоднородным дифференциальным уравнени- ем.Очевидно, что однородное уравнение (2.4) имеет тривиальное решение u = 0. Кроме того, любая функция, не зависящая от t и линейная относи- тельно x, удовлетворяет этому уравнению. Из этого следует, что уравнение теплопроводности имеет бесконечное множество решений.Для того чтобы получить единственное решение, к уравнению тепло- проводности добавляют обычно начальное и краевые условия.2.2. Постановка начально-краевых задачНачальное условие для уравнения теплопроводности состоит в задании значений функции u(x, t) в начальный момент времени t = 0:u(x, 0) =ϕ(x),гдеϕ(x) – функция, описывающая начальную температуру стержня. Кра- евые (граничные) условия задаются на торцевых концах стержня в точках x = 0 и x = l в соответствии с теплообменом стержня с окружающей средой.Рассмотрим различные случаи краевых условий.1. На концах стержня поддерживается заданная температура (в точке x = 0 – температура u0(t), а в точке x = l – температура u l(t)):u(0, t) = u0(t),u(l, t) = u l(t).Это условия первого рода, или условия Дирихле.2. Концы стержня теплоизолированы. Это означает, что поток тепла−k∂u∂x через единицу поверхности соответствующего конца стержня будет равен нулю. Тогда краевые условия на концах стержня можно задать сле- дующим образом:∂u(0, t)∂x= 0;∂u(l, t)∂x= 0.Полученные условия – это однородные условия второго рода, или условияНеймана.46 В общем случае могут быть известны значения теплового потока на границе: q0(t) (на конце x = 0) и q l(t) (на конце x = l). Тогда краевые условия второго рода записывают в виде k∂u(0, t)∂x= q0(t),−k∂u(l, t)∂x= q l(t).(Если тепло уходит в окружающую среду, то q0(t), q l(t) > 0.)3. На концах стержня происходит теплообмен с окружающей средой по закону Ньютона. В этом случае поток тепла пропорционален разности температур тела и окружающей среды и краевые условия имеют вид k∂u(0, t)∂x= h(u(0, t) − u0(t));−k∂u(l, t)∂x= h(u(l, t) − u l(t)).Здесь h – коэффициент теплоотдачи (h > 0); u0(t), u l(t) – температура среды на концах x = 0 и x = l соответственно. Такие краевые условия называются условиями третьего рода.Все рассмотренные краевые условия являются линейными. Они могут быть описаны одним соотношением на соответствующей границе:R1∂u(0, t)∂x− S1u(0, t) = g1(t),R2∂u(l, t)∂x+ S2(l, t) = g2(t).При R1= 0 (R2= 0) получается краевое условие первого рода, при S1= 0(S2= 0) – краевое условие второго рода, а при R1S1 6= 0 (R2S2 6= 0) –условие третьего рода.4. Если тело находится в вакууме, тогда изменение температуры на границе может происходить вследствие теплоизлучения по закону Стефана–Больцмана. Граничное условие в этом случае, например на конце x = 0,имеет вид k∂u(0, t)∂x=σu4(0, t) − u4 0(t) ,гдеσ– постоянная Стефана–Больцмана; u0– температура окружающей среды. Такое краевое условие, в отличие от рассмотренных ранее, является нелинейным.При решении прикладных задач на концах x = 0 и x = l могут быть поставлены краевые условия разного рода.5. Для бесконечного стержня (−∞ < x < +∞) предполагается, что температура в его бесконечно удаленных точках ограничена. Подобные за- дачи возникают при изучении процесса теплопроводности в очень длинном стержне. Температурный режим на концах такого стержня слабо влияет на температуру в его центральной части. Существенное значение имеет только начальное распределение температуры в стержне.47 В общем трехмерном случае, если изучается нагрев тела Ω с кусочно- гладкой границей Γ, уравнение теплопроводности с учетом источников или стоков тепла для функции u(M, t) (M ∈ Ω и t > 0) записывается в видеρc∂u∂t= div(k grad u) + g,где k – коэффициент теплопроводности; c – удельная теплоемкость;ρ–плотность тела; g – функция, описывающая объемную плотность источни- ков тепла. Если тело однородное и изотропное, т. е. можно считать, чтоρ,c, k – постоянные величины, то уравнение теплопроводности обычно пре- образуют к виду∂u∂t= a2∆u + f,где a2=kρc; f =gρc; ∆ = div(grad u) – оператор Лапласа. Вид оператора зависит от выбранной системы координат:∆u =∂2u∂x2+∂2u∂y2+∂2u∂z2(в декартовой),∆u =1ρ∂∂ρρ∂u∂ρ+1ρ2∂2u∂ϕ2+∂2u∂z2(в цилиндрической),∆u =1ρ2∂∂ρρ2∂u∂ρ+1ρ2sinθ∂∂θsinθ∂u∂θ+1ρ2sin2θ∂2u∂ϕ2(в сфери- ческой).Вывод трехмерного уравнения подробно описан в [2].К этому уравнению добавляется начальное условие при t = 0u(M, 0) =ϕ(M ),гдеϕ(M ) – функция, описывающая начальную температуру тела.На границе области Γ ставятся краевые условия в соответствии с усло- виями теплообмена тела с окружающей средой. Это может быть одно из трех условий вида:uΓ=µ(M, t),(2.6)∂u∂nΓ=ν(M, t),(2.7)∂u∂n+ huΓ=χ(M, t).(2.8)Здесьµ,ν,χ, h – заданные и непрерывные на границе Γ функции; M– точка границы области Γ;∂u∂n– производная по направлению внешней нормали к границе Γ.48 Условие (2.6) – это условие первого рода, которое ставится в том слу- чае, когда известна температура тела на границе. Условие (2.7) – это усло- вие второго рода, ставится тогда, когда на границе задан тепловой поток.Условие (2.8) – условие третьего рода, ставится, если теплообмен тела со средой происходит по закону Ньютона.Начально-краевая задача для уравнения теплопроводности с краевым условием первого рода на всей границе называется задачей Дирихле, а с краевым условием второго рода – задачей Неймана.На разных частях границы Γ могут быть заданы краевые условия раз- ного рода. В этом случае задача называется смешанной.2.3. Решение краевой задачи для уравнения теплопроводности методом ФурьеМетод Фурье является одним из наиболее распространенных методов решения линейных уравнений в частных производных.Рассмотрим начально-краевую задачу для одномерного уравнения теп- лопроводности. Пусть требуется найти функцию u(x, t), удовлетворяющую уравнению∂u∂t= a2∂2u∂x2+ f (x, t)(0 < x < l,t > 0),начальному условию u(x, 0) =ϕ(x),краевым условиямR1∂u(0, t)∂x− S1u(0, t) = g1(t),R2∂u(l, t)∂x+ S2u(l, t) = g2(t)(|R1| + |S1| 6= 0, |R2| + |S2| 6= 0). Решением поставленной задачи назовем функцию u(x, t), обладающую следующими свойствами:а) u(x, t) определена и непрерывна в области Ω = [0, l] × [0, T ];б) u(x, t) удовлетворяет уравнению теплопроводности при 0 < x < l,t > 0;в) u(x, t) удовлетворяет начальному условию и краевым условиям.Метод Фурье непосредственно применяется к задачам, в которых ис- комая функция удовлетворяет однородным краевым условиям (g1(t) ≡ 0 и g2(t) ≡ 0). Если краевые условия неоднородные, тогда сначала задачу сле- дует свести к задаче с однородными условиями. Для этого функцию u(x, t)представляют в виде u(x, t) = v(x, t) + w(x, t),49 где w(x, t) – дифференцируемая функция, удовлетворяющая тем же крае- вым условиям, что и функция u(x, t):R1∂w(0, t)∂x− S1w(0, t) = g1(t),R2∂w(l, t)∂x+ S2w(l, t) = g2(t).В этом случае функция v(x, t) будет удовлетворять однородным условиям(см. 1.1). Выбор функции w(x, t) зависит от типа граничных условий. Для многих краевых задач функция w(x, t) может быть задана как линейная по переменной x w(x, t) =α(t)x +β(t).Сумму функций v(x, t) + w(x, t) подставляют в дифференциальное уравнение и начальные условия, и всю краевую задачу записывают для функции v(x, t), которая удовлетворяет однородным краевым условиям.Задача сводится к нахождению функции v(x, t), удовлетворяющей уравне- нию∂v∂t= a2∂2v∂x2+ f v(x, t)(0 < x < l,t > 0),(2.9)начальному условию v(x, 0) =ϕv(x),(2.10)краевым условиямR1∂v(0, t)∂x− S1v(0, t) = 0,R2∂v(l, t)∂x+ S2v(l, t) = 0,(2.11)где f v(x, t) = f (x, t) −∂w(x, t)∂t+ a2∂2w∂x2,ϕv(x) =ϕ(x) − w(x, 0).Далее для задачи с однородными краевыми условиями метод Фурье применяется по следующей схеме:1. Для линейного дифференциального оператора второго порядкаLx(u) = −∂2u∂x2решают соответствующую задачу Штурма–Лиувилля:−y00(x) = λy(x),0 < x < l(y(x) 6≡ 0),R1y0(0) − S1y(0) = 0,R2y0(l) + S2y(l) = 0.В результате находят собственные числа {λk}+∞k=1и систему собственных функций {y k(x)}+∞k=1оператора краевой задачи.2. Функцию v(x, t) представляют в виде ряда Фурье по собственным функциям {y k(x)}+∞k=1v(x, t) =+∞Xk=1c k(t)y k(x).Однородные краевые условия (2.11) для v(x, t) при этом автоматически выполняются.Этот ряд подставляют в уравнение (2.9) и начальное условие50 (2.10), предварительно разложив функции f v(x, t) иϕv(x) в ряды Фурье по той же системе функций {y k(x)}+∞k=1. Далее, для коэффициентов c k(t)получают и решают задачи Коши.3. Искомую функцию u(x, t) записывают в виде u(x, t) = w(x, t) ++∞Xk=1c k(t)y k(x).Покажем, как применяется этот метод, на примерах.Пример 2.1. Найти температуру тонкого однородного стержня дли- ны l с теплоизолированной поверхностью, если его начальная температура u(x, 0) =T0x l, конец x = 0 теплоизолирован, а конец x = l поддерживается при постоянной температуре T0Функция u(x, t), описывающая распределение температуры в стержне,удовлетворяет уравнению∂u∂t= a2∂2u∂x2(0 < x < l,t > 0),начальному условию u(x, 0) =T0x lи краевым условиям∂u(0, t)∂x= 0,u(l, t) = T0Метод Фурье в случае данных краевых условий непосредственно непри- меним, так как условие на конце стержня x = l неоднородное. Поэто- му сведем сначала поставленную задачу к задаче с однородными краевы- ми условиями. Представим функцию u(x, t) в виде суммы двух функций u(x, t) = v(x, t) + w(x, t).Пусть w(x, t) =αx+β. Подберемαиβтак, чтобы для функции w(x, t)выполнялись краевые условия∂w(0, t)∂x=α= 0,w(l, t) =αl +β= T0Тогдаα= 0,β= T0, w(x, t) = T0и, следовательно, u(x, t) = v(x, t) + T0В этом случае справедливы равенства∂u∂t=∂v∂t,∂2u∂x2=∂2v∂x2,u(x, 0) = v(x, 0) + T0 51 Тогда для функции v(x, t) получится следующая задача с однородны- ми краевыми условиями:∂v∂t= a2∂2v∂x2(0 < x < l,t > 0),v(x, 0) =T0x l− T0,∂v(0, t)∂x= 0,v(l, t) = 0.Эту задачу решим методом Фурье.1. Функцию v(x, t) будем искать в виде ряда Фурье по собственным функциям линейного дифференциального оператораLx(u) = −∂2u∂x2второго порядка краевой задачи. Эти функции найдем, ре- шая соответствующую задачу Штурма–Лиувилля−y00(x) = λy(x),0 < x < l(y(x) 6≡ 0),y0(0) = 0,y(l) = 0.Собственные числа оператора −y00удовлетворяют условию λ ≥ 0.При λ = 0 общее решение уравнения y(x) = C1+ C2x. Подставив y(x)в краевые условия, получим C1= 0, C2= 0, тогда y(x) ≡ 0 и λ = 0 не является собственным числом оператора.При λ > 0 обозначим λ =µ2. Общим решением уравнения будет функ- ция y(x) = C1cos(µx) + C2sin(µx).Функция y(x) удовлетворяет однородным краевым условиям:(y0(0) = C2µ= 0,y(l) = C1cos(µl) + C2sin(µl) = 0.Из этих равенств следует, что приµ6= 0C2= 0иC1cos(µl) = 0.Так как C1 6= 0 (функция y(x) не может быть нулевой), то cos(µl) = 0.Полученное равенство справедливо для чиселµk l =π2+πk, т. е.µk=(2k + 1)π2l,k = 0, 1, 2, ... .Таким образом, для данной задачи Штурма–Лиувилля собственными числами и собственными функциями являются:λk=µ2k= (2k + 1)π2l2,y k(x) = cos(µk x),k = 0, 1, 2, ... .52 2. Функцию v(x, t) будем искать в виде ряда по найденным собствен- ным функциям дифференциального оператора:v(x, t) =+∞Xk=0c k(t)y k(x).Подставим ряд в уравнение теплопроводности и в начальное условие:+∞Xk=0c0k(t)y k(x) = a2+∞Xk=0c k(t)y00k(x),+∞Xk=0c k(0)y n(x) =+∞Xk=0ϕk yk(x).Здесь+∞Xk=0ϕk yk(x) – ряд Фурье для функцииϕv(x) =T0x l− T0, коэффици- енты Фурье которого вычисляются по правилуϕk=(ϕv, y k)ky kk2Используя равенство y00k(x) = −µ2k yk(x) и свойство единственности раз- ложения функции в ряд Фурье, получим задачи Коши для коэффициентов ck(t):(c0k(t) = −a2µ2k ck(t),c k(0) =ϕk(2.12)Найдем сначала коэффициентыϕk:ky kk2=lZ0cos2(µk x)dx =l2,(ϕv, y k) =lZ0 T0x l− T0cos(µk x)dx = −T0lµ2k,тогдаϕk= −2T0l2µ2k, k = 0, 1, 2, ... .Общее решение уравнения (2.12) c k(t) = Ak e−a2µ2k t. Используя равен- ство c k(0) =ϕk, получим Ak=ϕk, а значит,c k(t) = −2T0l2µ2k e−a2µ2k tи v(x, t) = −2T0l2+∞Xk=0e−a2µ2k tµ2k cos(µk x).53 3. Окончательно решение исходной задачи получим как сумму функ- ций w(x, t) и v(x, t):u(x, t) = T0−8T0π2+∞Xk=0e−a2µ2k t(1 + 2k)2cos (1 + 2k)πx2lПример 2.2. В цилиндрическом проводнике радиуса R вследствие прохождения постоянного тока в соответствии с законом Джоуля–Ленца выделяется тепловая энергия с объемной плотностью Q0. Теплоотдача с поверхности проводника происходит по закону Ньютона. Найти распреде- ление температуры по сечению проводника, если его начальная темпера- тура и температура внешней среды равны T0В примере изучается нагрев трехмерного цилиндрического тела. Вы- берем произвольное поперечное сечение проводника. Введем цилиндриче- скую систему координат, поместив ее начало в центре выбранного сечения и направив ось 0z вдоль оси проводника. Тогда процесс нагрева проводника описывается функцией u(ρ,ϕ, z, t), удовлетворяющей уравнению теплопро- водности∂u∂t= a2∆u + Q,где Q =Q0ρc,ρ– плотность, c – удельная теплоемкость проводника.Запишем уравнение в цилиндрической системе координат∂u∂t= a2 1ρ∂∂ρρ∂u∂ρ+1ρ2∂2u∂ϕ2+∂2u∂z2+ Q,0 <ρ< R,0 ≤ϕ< 2π,−∞ < z < +∞,t > 0.Распределение температуры в каждом сечении проводника можно счи- тать одинаковым (функция u от z не зависит), т. е. в уравнении∂2u∂z2= 0.Можно заметить, что распределение температуры в выбранном сечении зависит только от времени и от расстояния до центра проводника, т. е.∂2u∂ϕ2= 0. Тогда функцию u можно рассматривать как функцию двух пе- ременных u(ρ, t).Функция u(ρ, t) удовлетворяет уравнению теплопроводности:∂u∂t= a2 1ρ∂∂ρρ∂u∂ρ+ Q,0 <ρ< R,t > 0,начальному условию: u(ρ, 0) = T0и краевым условиям:u(ρ, t) ограничена приρ→ 0 + 0,∂u(R, t)∂ρ+ h(u(R, t) − T0) = 0 54 (h =h0k, где h0– коэффициент теплоотдачи, k – коэффициент теплопро- водности материала проводника).Условие ограниченности функции u приρ→ 0 + 0 является естествен- ным однородным краевым условием. Условие приρ= R неоднородное.Заменой u(ρ, t) = v(ρ, t) + w(ρ, t) сведем данную краевую задачу к краевой задаче относительно функции v(ρ, t) с однородными краевыми условиями.Положим w(ρ, t) = C = const. Первое краевое условие будет выполняться.Подставив w во второе краевое условие, убеждаемся в том, что оно будет выполняться при C = T0. Таким образом, w(ρ, t) = T0. Для функции v(ρ, t)получим следующую начально-краевую задачу:∂v∂t= a2 1ρ∂∂ρρ∂v∂ρ+ Q(0 <ρ< R,t > 0),v(ρ, 0) = 0,v(ρ, t) ограничена приρ→ 0 + 0,∂v(R, t)∂ρ+ hv(R, t) = 0.Краевые условия для функции v однородные, и метод Фурье непосред- ственно применим к этой задаче.1. Сначала найдем собственные функции и собственные числа опера- тора B0(v) = −1ρ∂∂ρρ∂v∂ρ– оператора Бесселя нулевого порядка. Для этого решим соответствующую задачу Штурма–Лиувилля:−1ρ(ρy0(ρ))0= λy(ρ),0 <ρ< R(y(ρ) 6≡ 0),y(ρ) ограничена приρ→ 0 + 0, y0(R) + hy(R) = 0.Дифференциальное уравнение можно записать в виде y00+1ρy0+ λy = 0.Это уравнение Бесселя нулевого порядка. Из общей теории следует, что собственные числа рассматриваемого оператора неотрицательны.При λ = 0 уравнение Бесселя становится уравнением Эйлера:ρ2y00+ρy0= 0.Выполнив подстановкуρ= e t, найдем y(ρ) = C1+ C2lnρ. Учитывая кра- евые условия,получим y(ρ) ≡ 0. Значит, λ = 0 не является собственным числом оператора.При λ > 0 решением уравнения Бесселя является функция y(ρ) = C1J0(√λρ) + C2N0(√λρ),55 где J0и N0– функции Бесселя и Неймана соответственно.Функция N0(√λρ) не ограничена приρ→ 0 + 0 (см. 1.5). Ограни- ченным в нуле решением уравнения будет функция y(ρ) = C1J0(√λρ).Подставим ее во второе краевое условие(J0(√λρ))0+ hJ0(√λρ)ρ=R= 0.Так как J0 0(ρ) = −J1(ρ), справедливо равенство−√λJ1(√λR) + hJ0(√λR) = 0.Обозначим√λR =γи получим уравнение относительноγ−γJ1(γ) + hRJ0(γ) = 0.(2.13)Покажем, что это уравнение имеет корни. Для этого запишем его в видеJ1(γ)J0(γ)=hRγКорни полученного уравнения являются абсциссами точек пересечения гра- фиков функций y =J1(γ)J0(γ)и y =hRγприγ> 0 (рис. 2.1).γγγγγ2 34 561   2   3   4   5   6   7   8   9   10

xxxxxlxtttttTtinjm0 00 12 21Рис. 2.2Множество точек плоскостиωhτ= {(x i, t j)}, i = 0, ..., n, j = 0, ..., m,называется сеткой, покрывающей область Ω. Точки плоскости (x i, t j) –это узлы сетки. При фиксированном j совокупность узлов сетки (x i, t j),i = 0, ..., n, называется j-м временным слоем. Будем считать, что функ- ции f (x, t),ϕ(x), g1(t), g2(t) из условия задачи определены только в узлах сетки.Непрерывной в области Ω функции u(x, t) поставим в соответствие сеточную функцию u(x i, t j) = u jiЗаданное дифференциальное уравнение, начальное и краевые условия запишем в узлах сетки. При этом каждую производную заменим разност- ным отношением, связывающим значения сеточной функции в несколь- ких узлах сетки. Производную по времени∂u∂t в узле сетки (x i, t j) можно заменить разностным выражением различными способами. Простейшими являются замены разностным отношением “вперед”u j+1i− u jiτ+ O(τ)!или “назад”u j−1i− u jiτ+ O(τ)!. В зависимости от способа аппроксимации67 производных разностными отношениями получаются различные разност- ные схемы для уравнения теплопроводности.Явная разностная схема.Для узлов сетки(x i, t j),i = 1, ..., n − 1, j = 0, ..., m − 1, запишем уравнение теплопроводности,заменив∂u∂t разностным отношением “вперед”u j+1i− u jiτ+ O(τ), а∂2u∂x2–разностным отношением uj i−1− 2u ji+ u ji+1h2+ O(h2):u j+1i− u jiτ+ O(τ) = a2u ji−1− 2u ji+ u ji+1h2+ O(h2) + f (x i, t j),i = 1, ..., n − 1,j = 0, ..., m − 1.(2.23)Для нулевого временного слоя (x i, t0), i = 0, ..., n, запишем начальное усло- вие u0i=ϕ(x i),i = 0, ..., n.В граничных узлах сетки (x0, t j+1) и (x n, t j+1), j = 0, ..., m − 1, запишем краевые условия, заменив производные соответствующими разностными отношениями:R1u j+1 1− u j+1 0h+ O(h) − S1u j+1 0= g1(t j+1),R2u j+1n− u j+1n−1h+ O(h) + S2u j+1m= g2(t j+1).Считая, что погрешности аппроксимации O(τ), O(h2), O(h) в урав- нениях и граничных условиях малы, отбросим их. В результате получим систему уравнений относительно неизвестных ˜u ji(приближенных значений функции u(x, t) в узлах сетки), которая называется явной разностной схе- мой для уравнения теплопроводности:˜u j+1i− ˜u jiτ= a2˜u ji−1− 2˜u ji+ ˜u ji+1h2+f (x i, t j), i = 1, ..., n−1, j = 0, ..., m−1,˜u0i=ϕ(x i), i = 0, ..., n,R1˜u j+1 1− ˜u j+1 0h− S1˜u j+1 0= g1(t j+1),R2˜u j+1n− ˜u j+1n−1h+ S2˜u j+1n= g2(t j+1),j = 0, ..., m−1.На рис. 2.3 указаны узлы сетки, которые используются при аппрокси- мации производных, входящих в уравнение. Такой рисунок принято назы- вать шаблоном схемы. Данная разностная схема является четырехточечной двухслойной. Схема называется явной, поскольку значения искомой функ- ции ˜u j+1i на (j + 1)-м временном слое можно последовательно в явном виде68 определить, если известны значения ˜u ji на j-м временном слое. Для этого полученные уравнения преобразуем к виду, удобному для вычислений:˜u j+1i= ˜u ji+τa2˜u ji−1− 2˜u ji+ ˜u ji+1h2+ f (x i, t j)!,i = 1, 2, ..., n − 1,˜u j+1 0=R1R1+ S1h˜u j+1 1−hR1+ S1h g1(t j+1),˜u j+1n=R2R2+ S2h˜u j+1n−1+hR2+ S2h g2(t j+1),˜u0i=ϕ(x i),i = 0, ..., n.(2.24)(xi −1, tj( xi, tj(xi,tj)))+1( xi,tj+1)Рис. 2.3Сначала находятся значения ˜u0i== u0i=ϕ(x i), i = 0, ..., n, на нулевом слое. Затем последовательно определяют- ся значения ˜u j+1i на всех временных сло- ях.Исследуем построенную разностную схему и найдем условия, при выполнении которых решение сеточной задачи схо- дится к решению исходной задачи для уравнения теплопроводности. Исследование сходимости проведем для за- дачи (2.20),(2.21), когда искомая функция u(x, t) удовлетворяет краевым условиям первого рода u(0, t) = g1(t),u(l, t) = g2(t).Разностная схема в этом случае примет вид˜u j+1 0= g1(t j+1),˜u j+1i= ˜u ji+τa2˜u ji−1− 2˜u ji+ ˜u ji+1h2+ f (x i, t j)!,i = 1, 2, ..., n − 1,˜u j+1n= g2(t j+1),˜u0i=ϕ(x i),i = 0, ..., n.(2.25)Сеточной нормой функции называется максимальное значение модуля69 функции на той сетке, на которой она задана:kuk h,τ= max i=0,...,n;j=0,...,m|u(x i, t j)|,kϕk h= max i=0,...,n|ϕ(x i)|.Разностная схема называется сходящейся к решению дифференциаль- ной задачи в сеточной норме, если ku − ˜uk h,τ→ 0 при h → 0 иτ→ 0.Сходимость разностной схемы связана с двумя важными понятиями:устойчивостью разностной схемы и точностью аппроксимации (см. 1.9).Разностная схема называется устойчивой по начальным данным и пра- вой части, если норма ее решения не превосходит суммы норм заданных в задаче функций, умноженных на число, не зависящее от шагов сетки:k˜uk h,τ≤ C(kϕk h+ kf k h,τ+ kg1kτ+ kg2kτ).(2.26)Это условие показывает, что при малом изменении исходных данных или при уменьшении шагов сетки решение разностной задачи будет меняться мало. Выполнение этого условия означает непрерывную зависимость реше- ния, полученного с помощью разностной схемы, от входных данных. При решении разностных уравнений, например (2.25), естественно использовать вычислительные средства. А значит, расчеты будут выполняться не точ- но, а с погрешностью, связанной с конечным представлением чисел (длина разрядной сетки ограничена).Если малые ошибки округления чисел приводят к большим искаже- ниям решения, то это явление называется вычислительной неустойчиво- стью разностной схемы. Разностная схема, обладающая таким свойством,не пригодна для вычислений. При исследовании разностной схемы обычно используется следующий подход. Считается, что ошибки, возникающие при выполнении арифметических операций, можно рассматривать как возму- щения, которые вносятся в начальные данные и в правую часть уравнения.От разностной схемы требуют, чтобы для нее выполнялось условие устой- чивости по начальным данным и правой части. Выполнение этого условия обеспечивает ее вычислительную устойчивость, при этом разностную схему называют устойчивой, а условие (2.26) – условием устойчивости разностной схемы.Докажем устойчивость явной разностной схемы при условии, что функ- ция u(x, t) удовлетворяет однородным краевым условиям первого рода:u(0, t) = 0,u(l, t) = 0(g1(t) = 0, g2(t) = 0).Утверждение 2.3. Еслиτ≤h2 2a2, то явная разностная схема (2.25)устойчива по начальным данным и правой части, т. е.k˜uk h,τ≤ C(kϕk h+ kf k h,τ),70 где C не зависит от h иτДоказательство. Обозначим k˜u jk h, kf jk hнаибольшие по модулю зна- чения сеточных функций ˜u и f на j-м временном слое:k˜u jk h= max i=0,...,n|˜u ji|,kf jk h= max i=0,...,n|f (x i, t j)|.Из равенств (2.24) следует, что|˜u j+1i| =τa2h2˜u ji−1+1 −2τa2h2˜u ji+τa2h2˜u ji+1+τf (x i, t j)Тогда|˜u j+1i| ≤τa2h2k˜u jk h+1 −2τa2h2k˜u jk h+τa2h2k˜u jk h+τkf jk hТак как ˜u j+1 0= ˜u j+1n= 0, то неравенство справедливо для всех i = 0, ..., n.Пусть 1−2τa2h2≥ 0 (т. е.τ≤h2 2a2), тогда, заменив1 −2τa2h2на 1−2τa2h2,получим|˜u j+1i| ≤ k˜u jk h+τkf jk hпри всех i = 0, ..., n. Значит,k˜u j+1k h≤ k˜u jk h+τkf jk hИспользуя оценку k˜u jk h≤ k˜u j−1k h+τkf j−1k hи т. д., придем к цепочке неравенств k˜u j+1k h≤ k˜u j−1k h+τkf j−1k h+ kf jk h ≤≤ k˜u j−1k h+τkf j−2k h+ kf j−1k h+ kf jk h ≤ ... ≤≤ k˜u0k h+τkf0k h+ kf1k h+ ... + kf jk h .Поскольку kf jk h≤ kf k h,τдля всех j = 0, 1, ..., m − 1 и k˜u0k h= kϕk h,получим k˜u j+1k h≤ kϕk h+τ(j + 1)kf k h,τ,гдеτ(j + 1) ≤ T (j = 0, ..., m − 1). Тогда для всех j = 0, ..., m − 1k˜u j+1k h≤ kϕk h+ T kf k h,τ,а значит,k˜uk h,τ≤ kϕk h+ T kf k h,τЯвная разностная схема будет устойчивой только при выполнении условияτ≤h2 2a2. В противном случае явная разностная схема не пригодна71 для вычислений. В случае a2= 1 условие устойчивости записывается в видеτ≤h2 2Если в систему разностных уравнений вместо значений сеточной функ- ции ˜u ji подставить значения точного решения дифференциальной задачи u(x i, t j), то для сохранения равенства в правую часть разностных урав- нений необходимо ввести дополнительное слагаемое. Если с уменьшением шагов сетки это слагаемое меняется как O(h k) + O(τm), то говорят, что разностная схема аппроксимирует задачу с погрешностью O(h k) + O(τm).Утверждение 2.4. Явная разностная схема (2.25) аппроксимирует задачу для уравнения теплопроводности (2.20), (2.21) с краевыми услови- ями u(0, t) = g1(t), u(l, t) = g2(t) с погрешностью O(h2) + O(τ).Доказательство. Запишем разностную схему для указанных крае- вых условий:˜u j+1i− ˜u jiτ= a2˜u ji−1− 2˜u ji+ ˜u ji+1h2+ f (x i, t j),i = 1, ..., n − 1,j = 1, ..., m − 1,˜u0i=ϕ(x i),i = 0, ..., n,˜u j+1 0= g1(t j+1),˜u j+1n= g2(t j+1),j = 0, ..., m − 1.Функции u(x, t) и ˜u ji в граничных узлах сетки (при t = 0, x = 0, x = l)удовлетворяют одинаковым условиям. Поэтому погрешность аппроксима- ции в данной задаче связана только с заменой производных разностными отношениями в уравнении теплопроводности. Из равенств (2.23) следует,что погрешность аппроксимации для рассматриваемой разностной схемы есть O(h2) + O(τ).Заметим, что разностная схема (2.24) аппроксимирует начально-крае- вую задачу для уравнения теплопроводности (2.20)–(2.22) с краевыми усло- виями третьего рода с погрешностью O(h) + O(τ). Это связано с тем, что появляется погрешность O(h) при аппроксимации производных в гранич- ных условиях.Докажем теперь сходимость явной разностной схемы (2.25).Утверждение 2.5. Приτ≤h2 2a2явная разностная схема (2.25) схо- дится в сеточной норме к решению задачи (2.20), (2.21) с краевыми усло- виями u(0, t) = g1(t), u(l, t) = g2(t).Доказательство. Доказательство проведем, используя точность ап- проксимации и устойчивость разностной схемы. Обозначим w = u − ˜u раз- ность между точным и приближенным решениями задачи. Очевидно, что72 во внутренних узлах сетки функция w будет удовлетворять уравнениям wj+1i− w jiτ= a2w ji−1− 2w ji+ w ji+1h2+ O(h2) + O(τ),i = 1, ..., n − 1,j = 0, ..., m − 1;на нулевом временном слое – начальному условию w0i= 0,i = 0, ..., n;в граничных узлах (x0, t j+1), (x n, t j+1) – краевым условиям wj+1 0= 0,w j+1n= 0,j = 0, ..., m − 1.В данном случае f (x, t) = O(h2)+O(τ) иϕ(x) = 0. Применив утверждение2.3 к разностной схеме для функции w, при выполнении условияτ≤h2 2a2получим kwk h,τ≤ T kO(h2) + O(τ)k h,τ≤ T kO(h2)k h,τ+ kO(τ)k h,τ .Поскольку справедливы оценки kO(h2)k h,τ≤ M1h2,kO(τ)k h,τ≤ M2τ(см. 1.9), то kwk h,τ≤ M (h2+τ),где M = T max(M1, M2), а значит, kwk h,τ→ 0 при h → 0 иτ→ 0.Сходимость доказана только при условии, что шаги сетки по времениτдостаточно малы (τ≤h2 2a2). Если условие устойчивости нарушается, то явная схема не будет сходиться к решению поставленной задачи. От слоя к слою значения ˜u j+1i начинают неограниченно расти. Схема “разваливается”.Неявная разностная схема. Построим для задачи (2.20)–(2.22) неяв- ную разностную схему. Для внутренних узлов сетки запишем уравнение теплопроводности, заменив∂u(x i, t j)∂t разностным отношением “назад”u ji− u j−1iτ+ O(τ), а∂2u(x i, t j)h2τ∂x2, как и для явной схемы, разностным от- ношением uj i−1− 2u ji+ u ji+1h2+ O(h2):u ji− u j−1iτ+ O(τ) = a2u ji−1− 2u ji+ u ji+1h2+ O(h2) + f (x i, t j),i = 1, ..., n − 1,j = 1, ..., m.Для слоя (x i, t0) запишем начальное условие:u0i=ϕ(x i),i = 0, ..., n,73 а для граничных узлов – краевые условия:R1u j1− u j0h+ O(h) − S1u j0= g1(t j),R2u jn− u jn−1h+ O(h) + S2u jn= g2(t j),j = 1, ..., m.Отбросив в уравнениях неизвестные погрешности O(h), O(τ), O(h2)получим систему относительно приближенных значений ˜u ji искомой функ- ции в узлах сетки:˜u ji− ˜u j−1iτ= a2˜u ji−1− 2˜u ji+ ˜u ji+1h2+ f (x i, t j),i = 1, ..., n − 1,j = 1, ..., m;начальное условие:˜u0i= u0i=ϕ(x i),i = 0, ..., n;краевые условия:R1˜u j1− ˜u j0h− S1˜u j0= g1(t j),R2˜u jn− ˜u jn−1h+ S2˜u jn= g2(t j),j = 1, ..., m.(xi −1, tj( xi, tj(xi +1,tj)(xi,tj −1)))Рис. 2.4Построенная разностная схема, как и явная, является четырехточечной двух- слойной разностной схемой. Ее шаблон изображен на рис. 2.4.Схема называется неявной, посколь- ку значения функции ˜u ji в данном случае не могут быть вычислены последователь- но через значения функции на (j − 1)-м временном слое. Для их нахождения на каждом временном слое следует решать систему уравнений−(R1+ S1h)˜u j0+ R1˜u j1= hg1(t j),˜u ji−1− (2 +h2τa2)˜u ji+ ˜u ji+1= −h2τa2˜u j−1i− h2f (x i, t j),i = 1, ..., n − 1,−R2˜u jn−1+ (R2+ S2h)˜u jn= hg2(t j).Значение ˜u0i для нулевого временного слоя определяется из начального условия ˜u0i=ϕ(x i) (i = 0, ..., n).74 Систему уравнений для j-го временного слоя удобно записать в мат- ричном видеB0C0 00 0A1B1C1 00 0A2B2C2 00 0... An−1Bn−1Cn−1 00 0AnBn˜u j0˜u j1˜u j2˜u jn−1˜u jn=F0F1F2Fn−1FnЗдесьB0= R1+ S1h, C0= R1, F0= hg1(t j),Ai= 1,Bi= −2 +h2τa2, Ci= 1,Fi= −h2τa2˜u j−1i− h2f (x i, t j),i = 1, ..., n − 1An= −R2, Bn= R2+ S2h,Fn= hg(t j).Кратко эту систему можно записать в видеA˜u j= Fj,где A – матрица коэффициентов; ˜u j– вектор значений функции ˜u; Fj–вектор значений правой части уравнений для j-го временного слоя. ВекторFj меняется от слоя к слою. Матрица коэффициентов системы A трехдиа- гональная. Для строк этой матрицы выполняется условие строгого диаго- нального преобладания |Bi| > |Ai| + |Ci| (i = 0, ..., n). Это условие обес- печивает однозначную разрешимость системы (см. 1.9). Кроме того, для решения системы уравнений с подобной матрицей коэффициентов можно применить метод Гаусса без перестановки строк и столбцов матрицы. Такой метод решения системы уравнений с трехдиагональной матрицей коэффи- циентов называется методом прогонки. Сначала выполняется прямой ход метода прогонки: вычисляются прогоночные коэффициенты˜C0=C0B0,˜F0=F0B0,˜Ci=CiBi− Ai˜Ci−1,i = 1, ..., n − 1,˜Fi=Fi− Ai˜Fi−1Bi− Ai˜Ci−1,i = 1, ..., n.При этом матрица коэффициентов системы преобразуется к верхней тре- угольной матрице, на главной диагонали которой стоят единицы; ˜Ci– нену- левые элементы матрицы, стоящие выше главной диагонали; ˜Fi– числа,стоящие в правой части системы уравнений.75 Затем выполняется обратный ход (решается система уравнений с верх- ней треугольной матрицей коэффициентов):˜u jn= ˜Fn,˜u ji= ˜Fi− ˜Ci˜u ji+1,i = n − 1, ..., 0.Метод прогонки можно применять только в том случае, когда в фор- мулах для вычисления прогоночных коэффициентов ˜Ci, ˜Fi знаменатели не обращаются в нуль. Это требование будет выполняться, если для матри- цы коэффициентов справедливо условие строгого диагонального преобла- дания. Таким образом, при применении неявной разностной схемы для на- хождения приближенного решения поставленной задачи значения ˜u ji будут определяться последовательно по временным слоям. Для каждого времен- ного слоя необходимо будет решать систему уравнений с трехдиагональной матрицей коэффициентов.Преимущество неявной схемы заключается в том, что в отличие от явной разностной схемы построенная неявная схема является безусловно устойчивой, т. е. она устойчива при любыхτи h. Устойчивость и сходи- мость неявной разностной схемы рассмотрим для задачи (2.20),(2.21) при условии,что функция u(x, t) удовлетворяет краевым условиям первого ро- да:u(0, t) = g1(t),u(l, t) = g2(t).Неявная разностная схема в этом случае имеет вид˜u j0= g1(t j),˜u ji−1− (2 +h2τa2)˜u ji+ ˜u ji+1= −h2τa2˜u j−1i− h2f (x i, t j),i = 1, ..., n − 1,˜u jn= g2(t j),˜u0i=ϕ(x i),i = 0, ..., n.(2.27)Пусть функция u(x, t) удовлетворяет однородным краевым условиям g1(t) = 0 и g2(t) = 0. Справедливо следующее утверждение.Утверждение 2.6. Неявная разностная схема (2.27) устойчива по начальным данным и правой части, т. е.k˜uk h,τ≤ C(kϕk h+ kf k h,τ),где C не зависит от h иτ76 Доказательство. Для всех внутренних узлов j-го временного слоя выполняется равенство˜u ji−1−2 +h2 2τ˜u ji+ ˜u ji+1= −h2τ˜u j−1i− h2f (x i, t j)или2 +h2τ˜u ji=˜u ji−1+ ˜u ji+1+h2τ˜u j−1i+ h2f (x i, t j)Тогда справедлива оценка2 +h2τ|˜u ji| ≤ |˜u ji−1| + |˜u ji+1| +h2τ|˜u j−1i| + h2|f (x i, t j)|.Пусть k˜u jk hи kf jk h– наибольшие по модулю значения сеточных функций˜u и f на j-м временном слое, тогда2 +h2τ|˜u ji| ≤ 2k˜u ji kh+h2τk˜u j−1i kh+ h2kf jk hФункция ˜u удовлетворяет однородным краевым условиям ˜u j0= ˜u jn= 0,поэтому последнее неравенство выполняется для всех i = 0, ..., n, и тогда2 +h2τk˜u jk h≤ 2k˜u ji kh+h2τk˜u j−1i kh+ h2kf jk hПреобразуем это неравенство:k˜u jk h≤ k˜u j−1k h+τkf jk hПоследовательно применяя его к (j−1)-му, (j−2)-му и т. д. слоям, получим:k˜u jk h≤ k˜u j−2k h+τkf j−1k h+ kf jk h ≤ ... ≤≤ k˜u0k h+τkf1k h+ kf2k h+ ... + kf jk h .Учитывая равенство k˜u0k h= kϕk hи оценки kf jk h≤ kf k h,τ,τj ≤ T для j = 1, ..., m, получим неравенство k˜u jk h≤ kϕk h+ T kf k h,τ,которое будет выполняться для всех j = 1, ..., m, а значит,k˜uk h,τ≤ kϕk h+ T kf k h,τОчевидно, что неявная разностная схема аппроксимирует задачу для уравнения теплопроводности с той же погрешностью, что и соответствую- щая явная схема. Если функция u(x, t) удовлетворяет краевым условиям77 первого рода u(0, t) = g1(t) и u(l, t) = g2(t), то погрешность аппроксимации разностной схемы (2.27) будет O(h2) + O(τ).Из устойчивости и аппроксимации следует сходимость неявной раз- ностной схемы. Доказательство сходимости абсолютно аналогично доказа- тельству сходимости явной разностной схемы.Схема Кранка–Николсона. Неявная разностная схема хотя и не имеет ограничения на шаги по времениτ, но поскольку ее погрешность аппроксимации равна O(τ) + O(h2), то разумно выбрать шагиτдостаточ- но малыми (порядка O(h2)). Более высокую точность аппроксимации по времени имеет схема Кранка–Николсона. Это шеститочечная двухслойная неявная разностная схема (рис. 2.5).1   2   3   4   5   6   7   8   9   10

ji−1+1( xij)( xi,tj)(xi +1, tj)( xi −1, tj −1)( xi, t)(x,tj −1−1, t)Рис. 2.5Построим эту схему для уравнения теплопроводности (2.20) при условии, что функция u(x, t) удовлетворяет начально- му условию (2.21) и краевым условиям первого рода u(0, t) = g1(t), u(l, t) = g2(t).Между узлами сетки (x i, t j) и (x i, t j+1)выбирается точка (x i, t j+τ2) (рис. 2.5).Для этой точки записывается уравнение теплопроводности. При этом∂u∂t заменя- ется центральным разностным отношением uj i− u j−1i2(τ/2)+Oτ22(см. 1.9),а∂2u∂x2в этой точке находится как полусумма uj i−1− 2u ji+ u ji+1h2+ O(h2) и uj−1i−1− 2u j−1i+ u j−1i+1h2+ O(h2). После отбрасывания погрешностей аппрокси- мации производных получается неявная разностная схема Кранка–Никол- сона:˜u ji− ˜u j−1iτ==˜u ji−1− 2˜u ji+ ˜u ji+1+ ˜u j−1i−1− 2˜u j−1i+ ˜u j−1i+1 2h2+ f (x i, t j−1+τ2),i = 1, ..., n − 1,j = 1, ..., m.Начальное условие:˜u0i=ϕ(x i),i = 0, ..., n.Граничные условия:˜u j0= g1(t j),˜u jn= g2(t j).Погрешность аппроксимации такой разностной схемы равна O(h2) + O(τ2).78 Иногда значение f (x i, t j−1+τ2) заменяют полусуммой f (x i, t j) + f (x i, t j−1)2Тогда для внутренних узлов сетки (x i, t j) схема Кранка–Николсона полу- чается как полусумма неявной схемы, записанной для узлов (x i, t j), и явной схемы, построенной для узлов (x i, t j−1). Погрешность аппроксимации раз- ностной схемы в этом случае также O(h2) + O(τ2).3. УРАВНЕНИЕ КОЛЕБАНИЙ СТРУНЫ.ВОЛНОВОЕ УРАВНЕНИЕ3.1. Уравнение колебаний струныРассмотрим туго натянутую струну, т. е. тонкую упругую (подчиняю- щуюся закону Гука) гибкую нить. Ось 0x направим вдоль струны. Предпо- ложим, что струну вывели из положения равновесия, например оттянули или ударили по ней. Будем считать, что возникающие при этом поперечные колебания струны происходят в одной плоскости 0xu и все точки струны движутся перпендикулярно оси 0x. Обозначим через u(x, t) отклонение то- чек струны от прямолинейной формы. Функция u(x, t) описывает попереч- ные колебания струны при t ≥ 0. График функции u(x, t) (рис. 3.1) зависит от момента времени t и представляет собой форму колеблющейся струны в данный момент времени t.xuBAαT( x , t)(∆,xt)+xα∆xT+F0xx1 2Рис. 3.1При фиксированном x функция u(x, t) описывает закон движения, про- изводная∂u(x, t)∂t– скорость движения, а вторая производная∂2u(x, t)∂t2–ускорение выбранной точки струны.Для вывода уравнения поперечных колебаний струны сделаем следу- ющие предположения:1. Будем считать, что струна совершает только малые колебания, т. е.ее форма незначительно отличается от прямой u = 0. Касательная в каж- дой точке x, проведенная к графику u(x, t), почти параллельна оси абсцисс.Пустьα(x, t) – угол, который образует касательная с осью 0x. Условие ма-79 лости колебаний означает, что величинойα2можно пренебречь. Разложим функции sinα, tgα, cosαпо формуле Тейлора первого порядка:sinα=α+ o(α),tgα=α+ o(α),cosα= 1 + o(α).При сделанном предположенииα≈ sinα≈ tgα=∂u∂x и cosα≈ 1.2. На любой выбранный участок струны действуют упругие силы на- тяжения. Выделим на струне произвольный малый участок [gAB](рис. 3.1).Его проекция на ось 0x есть [x, x + ∆x]. Будем считать, что упругие силы,приложенные к концам выбранного участка, направлены по касательной к графику u(x, t), их модули | T (x, t)| = | T (x + ∆x, t)| = T и не зависят от x и t.При изучении малых поперечных колебаний струны такое предполо- жение допустимо. Обозначимα1иα2углы, которые образуют касательные,проведенные в точках A и B соответственно к графику струны в момент времени t, с осью 0x. По условию все точки струны движутся параллельно оси 0u. Значит, силы, действующие на участок [gAB], таковы, что сумма их проекций на ось 0x должна быть равна нулю:−| T (x, t)| cosα1+ | T (x + ∆x, t)| cosα2= 0.Поскольку cosα1≈ cosα2≈ 1, то можно считать, что | T (x, t)| = T = const.3. В положении равновесия масса участка струны [x, x + ∆x] равнаρ∆x, гдеρ– линейная плотность струны. Считаем, что в процессе колеба- ний масса этого участка не меняется и масса [gAB] тоже равнаρ∆x.Такое предположение также допустимо, так как длина участка струны[gAB] в момент времени t|gAB| =x+∆xZx s1 + ∂u∂x2dx.При этом∂u∂x≈αи величинойα2можно пренебречь. Значит, |gAB| ≈ ∆x.4. Предположим,что на струну в плоскости колебаний действуют непре- рывно распределенные внешние силы, перпендикулярные оси 0x. Плот- ность распределения этих сил, рассчитанную на единицу длины, обозначим g(x, t). Если сила F , приложенная к участку длины ∆x, направлена вверх,то F (x, t) = g(x, t)∆xj, если – вниз, то F (x, t) = −g(x, t)∆xj.Рассмотрим выбранный малый участок [gAB]. Удалим части струны,расположенные справа и слева от него. Воздействие отброшенных частей струны заменим соответствующими силами натяжения. Тогда выделенный80 участок можно рассматривать как материальную точку, которая находится под воздействием трех сил T (x, t), T (x + ∆x, t) и F (x, t) = g(x, t)∆xj.Проецируя на ось 0u действующие на участок [gAB] силы и применяя закон Ньютона, получимρ∆x∂2u∂t2= T sinα2− T sinα1+ g(x, t)∆x.Согласно сделанным предположениям углыα1иα2малы и sinα1≈ tgα1=∂u(x, t)∂x,sinα2≈ tgα2=∂u(x + ∆x, t)∂xТогда, если использовать формулу Тейлора∂u(x + ∆x, t)∂x=∂u(x, t)∂x+∂2u(x, t)∂x2∆x + o(∆x)и отбросить бесконечно малую o(∆x), придем к равенствуρ∆x∂2u∂t2= T∂2u(x, t)∂x2∆x + g(x, t)∆x.Сокращая на ∆x, получим уравнение малых вынужденных поперечных колебаний струныρ∂2u∂t2= T∂2u∂x2+ g(x, t).Обычно это уравнение записывают в виде∂2u∂t2= a2∂2u∂x2+ f (x, t),(3.1)где a2=Tρ; f (x, t) =g(x, t)ρ. Полученное уравнение называется также одномерным волновым уравнением, или уравнением Даламбера.Если f (x, t) = 0, уравнение (3.1) называется однородным. Оно описы- вает свободные колебания струны без воздействия внешних сил.В случае f (x, t) 6= 0 уравнение называется неоднородным и описыва- ет вынужденные колебания струны. В частности, вынужденные колебания могут происходить под действием силы тяжести. Если же натяжение стру- ны T велико и на нее действует только сила тяжести, действием последней обычно пренебрегают и считают, что струна совершает свободные колеба- ния.81 Отметим, что многие физические задачи приводят к полученному вол- новому уравнению. Точно так же выглядит уравнение продольных колеба- ний тонкого упругого стержня. Функция u(x, t) при изучении таких коле- баний описывает продольные смещения точек стержня, имеющих в поло- жении равновесия абсциссу x. В уравнении (3.1) в этом случае a2=Eρ,f (x, t) =g(x, t)ρ, где E – модуль Юнга; g(x, t) – функция, описывающая плотность сил, действующих вдоль оси стержня. Вывод уравнения подроб- но описан в [9]. Кроме того, уравнение (3.1) совпадает с уравнением кру- тильных колебаний вала. Подобные уравнения появляются при изучении электрических колебаний и во многих других случаях. Меняется только физический смысл функций и коэффициентов, входящих в уравнение.В случае двух или трех пространственных переменных волновое урав- нение имеет вид∂2u∂t2= a2∆u + f,(3.2)где ∆ – оператор Лапласа; u и f – функции пространственных переменных и времени.Двумерное волновое уравнение получается, например, при изучении поперечных колебаний мембраны. Если мембрана прямоугольная, то урав- нение записывается в декартовой системе координат:∂2u∂t2= a2 ∂2u∂x2+∂2u∂y2+ f (x, y, t).При описании колебаний круглой мембраны переходят к полярной си- стеме координат и уравнение принимает вид∂2u∂t2= a2 1ρ∂∂ρρ∂u∂ρ+1ρ2∂2u∂ϕ2+ f (ρ,ϕ, t).Трехмерное волновое уравнение возникает, например, в задачах, свя- занных с изучением колебаний газа, находящегося в ограниченном объеме,в задачах распространения акустических волн и во многих других. В декар- товой системе координат это уравнение записывается следующим образом:∂2u∂t2= a2 ∂2u∂x2+∂2u∂y2+∂2u∂z2+ f (x, y, z, t).3.2. Постановка начальных и краевых условийВолновое уравнение имеет бесчисленное множество решений. Для од- нозначного описания колебательного процесса на функцию u следует на-82 ложить дополнительные условия, вытекающие из физического смысла за- дачи. Дополнительные условия – это начальное и краевые условия для функции u.Покажем как ставятся эти условия при изучении поперечных колеба- ний струны.Начальные условия обычно задают в момент времени t = 0. Условия описывают начальное положение точек струны и ее начальную скорость:u(x, 0) =ϕ(x),∂u(x, 0)∂t=ψ(x),гдеϕ(x) иψ(x) – заданные функции.Граничные условия показывают, что происходит на концах струны в течение всего времени колебаний. Предположим, что струна имеет конеч- ную длину l (0 ≤ x ≤ l).1. Если концевые точки струны движутся по определенному закону,то краевые условия имеют вид u(0, t) =ω1(t),u(l, t) =ω2(t).В частности, если концы струны закреплены, то для любого момента вре- мени t u(0, t) = 0,u(l, t) = 0.2. Граничные условия можно задать следующим образом:∂u(0, t)∂x= h1(t),∂u(l, t)∂x= h2(t).Это соответствует случаю, когда известен закон изменения касательных в концевых точках струны. Если h1(t) = 0 и h2(t) = 0, то в концевых точках для любого t ≥ 0 струна имеет касательные, параллельные оси 0x.3. В случае упругого закрепления концов стержня краевые условия записывают в виде∂u(0, t)∂x= h1(u(0, t) − u0(t)) ,∂u(l, t)∂x= h2(u(l, t) − u l(t)) ,где u0(t) и u l(t) – заданные функции.Все рассмотренные краевые условия являются линейными и описыва- ются уравнениямиR1∂u(0, t)∂x− S1u(0, t) = g1(t),R2∂u(l, t)∂x+ S2u(l, t) = g2(t).При R1= 0, R2= 0 – это краевые условия первого рода, или условияДирихле. При S1= 0, S2= 0 – условия второго рода, или условия Неймана,а при R1S1 6= 0, R2S2 6= 0 – условия третьего рода.83 Заметим, что при решении физических задач на концах x = 0 и x = l могут задаваться краевые условия разного рода.4. Если струна бесконечная, то для функции u(x, t) задают только на- чальные условия, краевые условия на нее обычно не накладываются, но при этом предполагается, что на бесконечности функция u(x, t) ограниче- на.Аналогично задаются начальные и краевые условия для уравнения(3.2). Пусть в области Ω с границей Γ функция u(M, t) (M ∈ Ω и t > 0)удовлетворяет волновому уравнению (3.2). К этому уравнению добавляют- ся начальные условия u(M, 0) =ϕ(M ),∂u(M, 0)∂t=ψ(M ).На границе Γ ставится одно из трех условий:uΓ=µ(M, t),∂u∂nΓ=ν(M, t),∂u∂n+ huΓ=χ(M, t),гдеµ,ν,χ, h – заданные и непрерывные на границе Γ функции; M – точка границы области Γ;∂u∂n– производная по направлению внешней нормали к границе.На разных частях границы Γ могут быть заданы краевые условия раз- ного рода. В таком случае краевые условия называются смешанными.3.3. Колебания бесконечной струны. Метод ДаламбераРассмотрим бесконечную струну, которую в начальный момент вре- мени вывели из положения равновесия. Будем считать, что внешние силы отсутствуют и струна совершает свободные колебания. Для нахождения поперечных колебаний такой струны следует решить волновое уравнение∂2u∂t2= a2∂2u∂x2,−∞ < x < +∞,t > 0,(3.3)при начальных условиях u(x, 0) =ϕ(x),∂u(x, 0)∂t=ψ(x),где функцииϕ(x),ψ(x) описывают начальное положение и начальную ско- рость точек струны.Сформулированная задача называется задачей Коши для бесконечной струны. Решим эту задачу методом, который называется методом Далам- бера, или методом бегущих волн. Покажем сначала, что общее решение уравнения (3.3) имеет вид u(x, t) = P (x − at) + Q(x + at),(3.4)84 где P и Q – произвольные дважды дифференцируемые функции. Введем новые переменныеξ= x − at иη= x + at и запишем волновое урав- нение (3.3) в новых переменных. Используя правило дифференцирования сложной функции, выразим производные функции u(x, t) по x и t через производные поξиη:∂u∂x=∂u∂ξ+∂u∂η,∂2u∂x2=∂2u∂ξ2+ 2∂2u∂η∂ξ+∂2u∂η2,∂u∂t= a ∂u∂η−∂u∂ξ,∂2u∂t2= a2 ∂2u∂η2− 2∂2u∂η∂ξ+∂2u∂ξ2Подставим выражения для∂2u∂x2и∂2u∂t2в уравнение (3.3), приведем подобные слагаемые и сократим на (−4a2). Тогда в новых переменных уравнение запишется в виде∂2u∂η∂ξ= 0,или∂∂η ∂u∂ξ= 0.Из этого следует, что∂u∂ξ= f (ξ),где f (ξ) – некоторая функция. Интегрируя последнее равенство, получим u =Zf (ξ)dξ+ Q(η),где Q(η) – произвольная функция.Обозначим P (ξ) =Zf (ξ)dξи подставим вместоξиηих выражения через x и t. В итоге получим равенство (3.4).Функции P (x − at) и Q(x + at) называются волнами отклонения. На- звание функций связано с их свойствами. Построим графики этих функций при t = 0: y = P (x) и y = Q(x). График функции y = P (x − at) получается параллельным переносом графика функции y = P (x) на at единиц вправо(a > 0). Соответственно, график функции y = Q(x + at) получается парал- лельным переносом графика y = Q(x) на at единиц влево. Таким образом,при непрерывном изменении t происходит перемещение графика функции y = P (x) вправо (рис. 3.2), а графика функции y = Q(x) влево (рис. 3.3).Для того чтобы решить поставленную задачу Коши, следует, поль- зуясь начальными условиями, определить неизвестные функции P и Q.85 =yP ( x −y = P ( )x)atatxy0Рис. 3.2−atyxyy0= Q ( x + at )= Q ( x )Рис. 3.3Продифференцируем функцию u(x, t) по t:∂u∂t= −aP0(x − at) + aQ0(x + at).Подставим выражения для u и∂u∂t в начальные условия, положив t = 0.В результате получим систему уравнений для функций P (x) и Q(x):(P (x) + Q(x) =ϕ(x),−aP0(x) + aQ0(x) =ψ(x).Интегрируя второе равенство в пределах от 0 до x, получим−a(P (x) − P (0)) + a(Q(x) − Q(0)) =xZ0ψ(ξ)dξОтсюда−P (x) + Q(x) =1a xZ0ψ(ξ)dξ+ C,где C = −P (0) + Q(0) – постоянная величина. Решая систему уравнений,найдемP (x) =1 2ϕ(x) −1 2a xZ0ψ(ξ)dξ−C2,Q(x) =1 2ϕ(x) +1 2a xZ0ψ(ξ)dξ+C2Подставляя полученные выражения для P (x) и Q(x) в (3.4), найдем функ- цию u(x, t):u(x, t) =1 2ϕ(x − at) −1 2a x−atZ0ψ(ξ)dξ+1 2ϕ(x + at) +1 2a x+atZ0ψ(ξ)dξ86 Заметим, что−x−atZ0ψ(ξ)dξ+x+atZ0ψ(ξ)dξ=0Zx−atψ(ξ)dξ+x+atZ0ψ(ξ)dξ=x+atZx−atψ(ξ)dξТогда функцию u(x, t), являющуюся решением поставленной задачи, мож- но представить в виде u(x, t) =ϕ(x − at) +ϕ(x + at)2+1 2a x+atZx−atψ(ξ)dξПолученное равенство называется формулой Даламбера решения за- дачи Коши для уравнения колебаний бесконечной струны.Найденное решение представляет собой сумму двух волн P (x − at) иQ(x + at). Одна волна “бежит” вправо, другая – влево. Число a =r Tρв уравнении колебаний струны называется скоростью распространения вол- ны.Функция u(x, t), полученная методом Даламбера, будет решением по- ставленной задачи при условии, что функцияϕ(x) дважды дифференци- руема, а функцияψ(x) дифференцируема один раз. В некоторых задачахϕ(x) иψ(x) не имеют нужных производных. Например, если струна в на- чальный момент времени имеет форму ломаной линии (рис. 3.4). В таких случаях считают, что формула Даламбера также дает решение задачи, хо- тя при этом функция u(x, t) не всюду дважды дифференцируема. Такое решение называют обобщенным решением задачи.Рассмотрим несколько примеров, иллюстрирующих применение мето- да Даламбера.Пример 3.1. Решить задачу Коши для уравнения колебаний неогра- ниченной струны:∂2u∂t2= a2∂2u∂x2,t > 0,−∞ < x < +∞,u(x, 0) = he−x2,∂u(x, 0)∂t= ve−x2Для поставленной задачиϕ(x) = he−x2,ψ(x) = ve−x2. Используя фор- мулу Даламбера, получим u(x, t) =he−(x−at)2+ he−(x+at)2 2+1 2a x+atZx−at ve−ξ2dξ=87 =h2e−(x−at)2+ e−(x+at)2+v√π4a(erf(x + at) − erf(x − at)) ,где erf x =2√πxZ0e−ξ2dξ– функция ошибок.Пример 3.2. Изобразить форму бесконечной струны для моментов времени t0= 0, t1=l2a, t2=l a, t3=2l a, если начальная скорость точек струны равна нулю, а начальная форма – это треугольник на отрезке [−l, l]с максимальным отклонением h (рис. 3.4).1   2   3   4   5   6   7   8   9   10

uhxluhxluhxluhtxluhxlt = tt = tt = t= tt = t0 04 32 10Рис. 3.10Процесс колебаний струны изображен на промежутке времени t ∈0;l a. В следующий промежуток времени t ∈ l a;2l aструна прой- дет все отмеченные положения в обратном порядке и вернется в первона- чальное положение, после чего процесс будет периодически повторяться с периодом T =2l aПример 3.6. Найти закон колебаний конечной струны (0 ≤ x ≤ l),если в начальный момент времени по струне ударили и ее точкам на про- межуткеα≤ x ≤βпридали скорость v.101 В данном случае функция u(x, t), описывающая свободные колебания струны,является решением задачи(3.8)–(3.10),при этомϕ(x) = 0 иψ(x) =(v,x ∈ [α,β],0,x 6∈ [α,β].Найдем для этих функций коэффициенты Фурье. Очевидно, чтоϕk= 0.Для определенияψk вычислим скалярное произведение(ψ(x), y k(x)) =lZ0ψ(x) sinπkx ldx =βZαv sinπkx ldx == −vlπk cosπkx l βα=vlπkcosπkαl− cosπkβlУчитывая, что ky k(x)k2=l2, получим решение поставленной задачи u(x, t) =2vlπ2a+∞Xk=1 1k2cosπkαl− cosπkβlsinπkat lsinπkx l4. УРАВНЕНИЯ ЛАПЛАСА И ПУАССОНА4.1. Определения. Постановка краевой задачиЛинейное дифференциальное уравнение второго порядка вида∆u = f,где ∆ – оператор Лапласа, f – заданная функция, называется уравнениемПуассона. Если в уравнении f = 0, то оно называется уравнением Лапласа.Такие уравнения описывают стационарное распределение тепла в од- нородном изотропном теле, электростатический и гравитационный потен- циалы в точках пространства и другие стационарные явления.Пусть функция u = u(M ) рассматривается в области Ω с кусочно- гладкой границей Γ. Часто ставится следующая задача: найти дважды дифференцируемую функцию u, удовлетворяющую в области Ω уравне- нию Лапласа или Пуассона, а на границе Γ – одному из краевых условий:1)uΓ=µ– краевое условие первого рода (условие Дирихле);2)∂u∂nΓ=ν– краевое условие второго рода (условие Неймана);102 3)u + h∂u∂nΓ=χ– краевое условие третьего рода.Здесьµ,ν,χи h – заданные и непрерывные на границе Γ функции.Название краевой задачи связано с краевым условием, которому удо- влетворяет функция u(M ) на границе Γ.Краевая задача∆u = f,uΓ=µназывается задачей Дирихле для уравнения Пуассона. Если на границе задано условие∂u∂nΓ=ν,то краевая задача называется задачей Неймана, а в случае краевого усло- вия видаu + h∂u∂nΓ=χ– третьей краевой задачей.В некоторых случаях на разных частях границы Γ могут быть заданы краевые условия разного рода. Тогда краевая задача называется смешан- ной.4.2. Гармонические функцииОпределение 4.1. Функция u = u(M ) называется гармонической в области Ω, если она имеет непрерывные частные производные первого и второго порядков в области Ω и удовлетворяет уравнению Лапласа∆u = 0.Приведем примеры различных гармонических фунций. При решении задач обычно используются декартовы, полярные, цилиндрические или сфериче- ские координаты. Вид оператора Лапласа зависит от выбранной системы координат (см. 2.2).Пример 4.1. Рассмотрим уравнение Лапласа, в котором оператор Ла- пласа записан в декартовой системе координат:∂2u∂x2+∂2u∂y2+∂2u∂z2= 0.Очевидно, что любая линейная функция вида u(x, y, z) = Ax + By + Cz + Dявляется гармонической функцией.103 Пример 4.2. Пусть (ρ,ϕ) – полярные координаты точек плоскости.Рассмотрим уравнение Лапласа, записав оператор Лапласа в полярной си- стеме координат:1ρ∂∂ρρ∂u∂ρ+1ρ2∂2u∂ϕ2= 0.Найдем u = u(ρ) – решение этого уравнения, которое не зависит от угловой координатыϕ(осесимметричное решение). Решение, обладающее этим свойством, удовлетворяет уравнению1ρ∂∂ρρ∂u∂ρ= 0⇔ρ∂u∂ρ= C1Преобразуем его к виду∂u∂ρ=C1ρПроинтегрировав, найдем гармоническую функцию u(ρ) = C1ln(ρ) + C2Здесь C1, C2– произвольные константы. Зададим C1= −1, C2= 0 и получим функцию u(ρ) = ln1ρЭту функцию называют фундаментальным решением уравнения Лапласа на плоскости. Она является гармонической функцией на всей плоскости,кроме начала координат. Обозначим r = |M M0| =p(x − x0)2+ (y − y0)2расстояние между точками M (x, y) и M0(x0, y0). Нетрудно заметить, что функция u(r) = ln1r также будет гармонической функцией на всей плоскости, исключая точкуM0Пример 4.3. Пусть (ρ,θ,ϕ) – сферические координаты точек про- странства. Запишем теперь оператор Лапласа в сферической системе ко- ординат:1ρ2∂∂ρρ2∂u∂ρ+1ρ2 1sinθ∂∂θsinθ∂u∂θ+1ρ2sin2θ∂2u∂ϕ2= 0.Будем искать решение u = u(ρ) уравнения Лапласа, которое не зависит от угловых координатθиϕ(центрально-симметричное решение). Очевидно,что u(ρ) удовлетворяет уравнению1ρ2∂∂ρρ2∂u∂ρ= 0⇔ρ2∂u∂ρ= C1 104 Разделим наρ2:∂u∂ρ=C1ρ2После интегрирования получим функцию u(ρ) =−C1ρ+ C2, гармоническую при любых константах C1и C2. Пусть C1= −1, C2= 0. Функция u(ρ) =1ρназывается фундаментальным решением уравнения Лапласа в простран- стве. Она является гармонической функцией на всей плоскости, кроме на- чала координат. Пусть r = |M M0| =p(x − x0)2+ (y − y0)2+ (z − z0)2– расстояние между точками M (x, y, z) и M0(x0, y0, z0). Очевидно, что функ- ция u(r) =1r также является гармонической функцией в пространстве, исключая точкуM0 4.3. Формулы ГринаПусть Ω ∈ R3– область, ограниченная поверхностью Γ.Утверждение 4.1. Если функции u = u(M ) и v = v(M ) имеют непрерывные частные производные второго порядка в области Ω и непре- рывны вместе с частными производными первого порядка на множествеΩS Γ, то справедливы формулы:Z Z ZΩ(grad v, grad u) dV +Z Z ZΩv∆u dV =Z ZΓv∂u∂n ds,(4.1)Z Z ZΩ(v∆u − u∆v) dV =Z ZΓv∂u∂n− u∂v∂nds,(4.2)называемые первой и второй формулами Грина соответственно. Здесь∂u/∂n = (grad u, n) – производная функции u по направлению внешней нормали n к поверхности Γ.Доказательство. Докажем сначала первую формулу Грина (4.1).Воспользуемся формулой Остроградского–ГауссаZ Z ZΩdiv F dV =Z ZΓ( F , n) ds.(4.3)105 Зададим F = v grad u. Поскольку верны следующие равенства:∂∂xv∂u∂x=∂v∂x∂u∂x+ v∂2u∂x2,∂∂yv∂u∂y=∂v∂y∂u∂y+ v∂2u∂y2,∂∂zv∂u∂z=∂v∂z∂u∂z+ v∂2u∂z2,то, складывая их почленно, получим div F = div(v grad u) = (grad v, grad u) + v∆u.При этом( F , n) = (v grad u, n) = v(grad u, n) = v∂u∂nС учетом полученных равенств из формулы Остроградского–Гаусса(4.3) следует первая формула ГринаZ Z ZΩ(grad v, grad u) dV +Z Z ZΩv∆u dV =Z ZΓv∂u∂n ds.Поменяем теперь местами функции u и v в первой формуле Грина:Z Z ZΩ(grad u, grad v) dV +Z Z ZΩu∆v dV =Z ZΓu∂v∂n ds.Вычтем почленно из первой формулы Грина это равенство:Z Z ZΩ(v∆u − u∆v) dV =Z ZΓv∂u∂n− u∂v∂nds.Получилась вторая формула Грина (4.2).4.4. Свойства гармонических функцийИспользуя полученные формулы, выделим некоторые свойства гармо- нических функций.Утверждение 4.2. Если u = u(M ) – функция, гармоническая в об- ласти Ω, то производная этой функции по нормали к поверхности Γ,ограничивающей область Ω, удовлетворяет условиюZ ZΓ∂u∂n ds = 0.(4.4)106 Доказательство. Воспользуемся первой формулой Грина (4.1). По- скольку u – гармоническая в области Ω функция, то ∆u = 0 и выполняется равенствоZ Z ZΩ(grad v, grad u) dV =Z ZΓv∂u∂n ds.Зададим функцию v = 1, тогда, учитывая, что grad 1 = 0, получим дока- зываемое утверждение.Рассмотрим функцию v(M ) =1r, где r = |M M0| =p(x − x0)2+ (y − y0)2+ (z − z0)2– расстояние между точками M (x, y, z) и M0(x0, y0, z0). В 4.2 было показа- но, что она является гармонической функцией, т. е. удовлетворяет уравне- нию Лапласа ∆u = 0 в пространстве, исключая точку M0Пусть M0(x0, y0, z0) – это некоторая внутренняя точка множества Ω.Утверждение 4.3. Если u = u(M ) – функция, гармоническая в обла- сти Ω, непрерывная вместе с частными производными первого порядка на множестве ΩS Γ, то для любой внутренней точки M0(x0, y0, z0) мно- жества Ω справедлива формула u(M0) =1 4πZ ZΓ 1r∂u∂n− u∂∂n 1rds,(4.5)называемая основной формулой теории гармонических функций. Здесь n– внешняя нормаль к поверхности Γ.Доказательство. Воспользуемся второй формулой Грина (4.2). По- скольку функция v(M ) =1r имеет разрыв в точке M0(x0, y0, z0), то сразу применить к функциям u и v вторую формулу Грина нельзя. Окружим точ- ку M0шаром Kε(M0) с центром в точке M0малого радиусаε. Рассмотрим область Ω \ Kε(M0). Граница этой области Γ ∪ Γε, где Γε– сфера, огра- ничивающая шар Kε(M0). Применим вторую формулу Грина к функциям u(M ) и v(M ) =1r в области Ω \ Kε(M0):Z Z ZΩ\Kε(M0) 1r∆u − u∆ 1rdV ==Z ZΓ 1r∂u∂n− u∂∂n 1rds +Z ZΓε1r∂u∂n ds −Z ZΓεu∂∂n 1rds.(4.6)107 Поскольку u(M ) и v(M ) =1r– гармонические в области Ω \ Kε(M0) функ- ции, то левая часть равенства равна нулю.Второй интеграл в правой части вычисляется по сфере Γε, поэтому множитель1r=1εне зависит от M и его можно вынести за знак интеграла.Интеграл преобразуется к виду1εZ ZΓε∂u∂n ds.Функция u(M ) – гармоническая, поэтому, согласно утверждению 4.2, зна- чение этого интеграла равно нулю.Рассмотрим в правой части формулы (4.6) третий интеграл по сфе- ре Γε. Поскольку внешняя нормаль к этой поверхности n направлена по радиусу к центру сферы, то∂∂n 1r Γε= −∂∂r 1r Γε=1ε2не зависит от M .Тогда получится равенствоZ ZΓ 1r∂u∂n− u∂∂n 1rds =1ε2Z ZΓεu ds.Применим теорему о среднем к интегралу, стоящему в правой части.Функция u(M ) непрерывна, следовательно, на сфере найдется такая точ- ка Mε, что интеграл по Γεбудет равен значению u(Mε), умноженному на площадь сферы 4πε2. Тогда1ε2Z ZΓεu ds = 4πu(Mε).Устремимεк нулю. При этом сфера Γε(M0) стягивается в точку M0:limε→0u(Mε) = u(M0).В итоге, получим равенство u(M0) =1 4πZ ZΓ 1r∂u∂n− u∂∂n 1rds.108 Замечание.Полученная формула справедлива для гармонических функций трех переменных.На плоскости решением уравнения Лапласа является функция u(r) = ln1r,где r = |M M0| =p(x − x0)2+ (y − y0)2Дословно повторяя описанные ранее рассуждения, нетрудно получить ос- новную формулу теории гармонических функций для гармонических функ- ций двух переменных:u(M0) =1 2πZ ZΓln1r ∂u∂n− u∂∂nln1rds.Пусть KR(M0) – шар радиуса R с центром в точке M0, а ΓR– сфера,ограничивающая этот шар.Утверждение 4.4. Если u = u(M ) – функция, гармоническая в шареKR(M0) и непрерывная вместе с частными производными первого поряд- ка на множестве KR(M0)S ΓR, то ее значение в центре M0шара есть среднее ее значений на сфере:u(M0) =1 4πR2Z ZΓRu ds.Доказательство. Поскольку u(M ) – гармоническая функция, то со- гласно утверждению 4.3 выполняется равенство u(M0) =1 4πZ ZΓR 1r∂u∂n− u∂∂n 1rds.Поверхность ΓR– это сфера. Для точек, лежащих на ней, r = |M0M | = R.Внешняя нормаль к поверхности сферы направлена по радиусу от центра сферы, поэтому∂∂n 1r ΓR=∂∂r 1r ΓR= −1R2Тогда u(M0) =1 4πZ ZΓR 1R∂u∂n+ u1R2ds =1 4πRZ ZΓR∂u∂n ds +1 4πR2Z ZΓRu ds.Из утверждения 4.2 следует, чтоRRΓR∂u∂n ds = 0, поэтому получаем доказы- ваемое равенство.109 Утверждение 4.5. Функция u = u(M ), гармоническая в ограничен- ной области Ω и непрерывная в области Ω ∪ Γ, достигает наибольшего и наименьшего значений на границе Γ этой области.Это утверждение обычно называют принципом максимума.Доказательство. 1. Если u(M ) = C = const, то теорема очевидна.2. Пусть u(M ) 6= const. Докажем теорему методом рассуждения от противного. Предположим, что во внутренней точке M0области Ω функция u имеет локальный экстремум, например локальный максимум. Окружим точку M0шаром KR(M0) достаточно малого радиуса R с границей ΓR,такого, что шар целиком лежит в области Ω и для всех M ∈ KR(M0) ∪ ΓRвыполняется неравенство u(M0) > u(M ).Поскольку u – гармоническая в шаре KR(M0) функция, то согласно утверждению 4.4 справедливо равенство u(M0) =1 4πR2Z ZΓRu ds.Здесь ΓR– сфера, ограничивающая шар. Применим теорему о среднем к интегралу, стоящему в правой части. Так как функция u(M ) непрерывна,то на сфере найдется такая точка MR, что интеграл по ΓRбудет равен значению u(MR), умноженному на площадь сферы 4πR2. Тогда u(M0) =1 4πR2 4πR2u(MR) = u(MR).Получилось противоречие с условием u(M0) > u(M ). Значит, внутренние точки области Ω не являются точками экстремума функции, и поскольку u(M ) непрерывна в области Ω ∪ Γ, то она достигает наибольшего и наи- меньшего значений на границе Γ области Ω.4.5. Теоремы единственности для уравненийЛапласа и ПуассонаУтверждение 4.6. Решение u = u(M ) задачи Дирихле для уравне- ния Пуассона∆u = f,u|Γ=µединственно.Доказательство. Предположим, что две функции u1и u2являются решениями задачи Дирихле. Рассмотрим функцию u = u1− u2. Эта функ- ция – гармоническая в области Ω. Она непрерывна в области Ω ∪ Γ и равна нулю на границе Γ:∆u = 0,u|Γ= 0.110 Согласно принципу максимума наибольшее и наименьшее значения функ- ции u достигаются на границе и равны нулю. Следовательно, u = u1− u2== 0, т.е. u1= u2Рассмотрим теперь задачу Неймана для уравнения Пуассона:∆u = f,∂u∂nΓ=νРешение этой задачи существует только при выполнении условияZ Z ZΩf dV =Z ZΓνds.Утверждение 4.7. Решение u = u(M ) задачи Неймана для уравне- ния Пуассона единственно с точностью до постоянного слагаемого C.Доказательство.Подставив в первую формулу Грина (4.1) функцию v = u, получим равенствоZ Z ZΩu∆u dV +Z Z ZΩ|grad u|2dV =Z ZΓu∂u∂n ds.(4.7)Предположим, что есть 2 решения u1и u2задачи Неймана. Тогда u == u1− u2является гармонической функцией (∆u = 0), удовлетворяющей краевому условию∂u∂nΓ= 0. Тогда равенство (4.7) преобразуется к видуZ Z ZΩ|grad u|2dV = 0.Следовательно, grad u = 0, т. е. u = C = const. Значит, u1= u2+ C.Утверждение 4.8. Решение u = u(M ) третьей краевой задачи для уравнения Пуассона∆u = f,u + h∂u∂nΓ=χ,h > 0единственно.Доказательство. Предположим, что третья краевая задача имеет2 решения u1и u2. Рассмотрим функцию u = u1− u2. Подставим ее в формулу (4.7). Тогда, учитывая, что u = u1−u2является решением задачи∆u = 0,u + h∂u∂nΓ= 0,111 получим для этой функции равенствоZ Z ZΩ|grad u|2dV + hZ ZΓ ∂u∂n2ds = 0.Поскольку h > 0, то равенство возможно только тогда, когда значения обоих слагаемых в левой части равенства равны нулю. Следовательно,grad u = 0 в области Ω и∂u∂n= 0 на границе Γ. Из условия grad u = 0 сле- дует, что u = C = const в области Ω. Функция u удовлетворяет краевому условиюu + h∂u∂nΓ= 0. И поскольку выполняется условие∂u∂nΓ= 0,то u = C = 0. Значит, u1= u2 4.6.Метод Фурье для уравнений Лапласа и ПуассонаВ некоторых областях Ω решение краевой задачи для уравнения Ла- пласа или Пуассона можно получить в виде ряда Фурье по собственным функциям оператора Штурма–Лиувилля. Алгоритм нахождения решения такой же, как при решении уравнения теплопроводности или волнового уравнения. Рассмотрим пример решения задачи Дирихле для уравненияЛапласа в шаре (осесимметричный по переменнойϕслучай).Пример 4.4. Найти стационарное распределение температуры в шаре радиуса R, верхняя половина поверхности которого поддерживается при температуре T0, а нижняя – при температуре, равной нулю. Определить температуру точек шара, расположенных вдоль радиусаθ= 0.Пусть u – искомая функция, описывающая стационарное распределе- ние температуры в шаре. Она является решением краевой задачи∆u = 0,uΓ= u0,где u0– функция, описывающая распределение температуры на границе шара.В модели будем использовать сферическую систему координат. Усло- вия задачи таковы, что функция u = u(ρ,θ) не зависит от угловой коор- динатыϕ. Она удовлетворяет уравнению Лапласа1ρ2∂∂ρρ2∂u∂ρ+1ρ2sinθ∂∂θsinθ∂u∂θ= 0и краевым условиям: u(ρ,θ) ограничена приθ→ 0 + 0 и приθ→π− 0;u(ρ,θ) ограничена приρ→ 0 + 0, u(R,θ) = u0(θ),где u0(θ) =T0,0 ≤θ≤π2,0,π2<θ≤π112 По переменнойθфункция u = u(ρ,θ) удовлетворяет естественным краевым условиям ограниченности функции. Эти условия являются одно- родными, поэтому решение задачи будем искать в виде ряда Фурье, разло- жив функцию u(ρ,θ) по собственным функциям оператораLθ(y) = −1sinθd dθ(sinθdy dθ).1. Решим задачу Штурма–Лиувилля для этого оператора:−1sinθ(sinθy0)0= λy,0 <θ<π,y(θ)ограничена приθ→ 0 + 0и приθ→π− 0.Решение этой задачи подробно разобрано в пособии [1] (пример 2.7).Собственные числа оператора: λk= k(k + 1), k = 0, 1, 2, ...; собствен- ные функции: y k(θ) = Pk(cosθ),k = 0, 1, 2, ..., где Pk(x) – многочленыЛежандра (Pk(x) =1 2k k!(x2− 1)k(k)); квадраты норм собственных функ- ций: ky k(θ)k2=2 2k + 1. Собственные функции образуют полную ортого- нальную систему в пространстве L2[0,π; sinθ].2. Будем искать функцию u(ρ,θ) в виде ряда Фурье по системе соб- ственных функций:u(ρ,θ) =+∞Xk=0c k(ρ)y k(θ).Функцию u0(θ) из краевого условия разложим в ряд Фурье по той же системе функций:u0(θ) =+∞Xk=0u0k yk(θ),где u0k=(u0, y k)||y k||2. Найдем скалярные произведения(u0, y k) =πZ0u0(θ)y k(θ) sinθdθ=π2Z0T0Pk(cosθ) sinθdθВыполнив замену x = cosθполучим(u0, y k) = T0 1Z0Pk(x)dx.113 Для вычисления интеграла воспользуемся свойствами многочленов Лежанд- ра. Для этих многочленов выполняются рекуррентные соотношения (1.34),(1.35):(k + 1)Pk+1(x) − (2k + 1)xPk(x) + Pk−1(x) = 0,(2k + 1)Pk(x) = P0k+1(x) − P0k−1(x)и справедливы равенства Pk(1) = 1, k = 0, 1, ... . Тогда(u0, y k) =T0 2k + 1 1Z0(P0k+1(x) − P0k−1(x))dx ==T0 2k + 1(Pk−1(0) − Pk+1(0)) =k + 2k + 1Pk−1(0),k = 1, 2, ... .Отдельно найдем(u0, y0) = T0 1Z0P0(x)dx = T0 1Z0 1 dx = T0Получим теперь коэффициенты Фурье u0k:u0 0=T0 2,u0k=T0 2k + 2k + 1Pk−1(0),k = 1, 2, ... .Найдем значения коэффициентов c k(ρ), подставив ряды Фурье функ- ций u(ρ,θ) и u0(θ) в уравнение Лапласа и в краевые условия, заданные по переменнойρ. Используя равенства −1sinθ(sinθ· y k(θ)0)0= λy k(θ) и свой- ство единственности разложения функции в ряд Фурье, получим краевые задачи для коэффициентов c k(r):((ρ2c0k(ρ))0= k(k + 1)c k(ρ),c k(ρ) ограничена приρ→ 0 + 0,c k(R) = u0kУравнение (ρ2c0k(ρ))0= k(k + 1)c k(ρ) ⇔ρ2c00k(ρ) + 2ρc0k(ρ) − k(k + 1)c k(ρ) = 0является уравнением Эйлера. Используем подстановкуρ= e t. Тогда спра- ведлива цепочка равенств c k(ρ) = c k(e t) = z k(t) = z k(lnρ). Нетрудно полу- чить уравнение, которому удовлетворяет функция z k(t):z00k(t) + z0k(t) − k(k + 1)z k(t) = 0.Запишем его характеристическое уравнение: λ2+ λ − k(k + 1) = 0. Найдем корни этого уравнения: λ1= −(k + 1), λ2= k. Общим решением диффе- ренциального уравнения будет функция z k(t) = Ak e−(k+1)+ Bk ek. Значит,c k(ρ) = Akρ−(k+1)+ Bkρk,k = 1, 2, ... .114 Используя первое краевое условие, получим Ak= 0 (k = 1, 2, ...). Из вто- рого условия следует, что BkRk= u0k, тогда Bk=u0kRk. В итоге ck(ρ) =T0 2k + 2k + 1Pk−1(0)ρRk,k = 1, 2, ... .Найдем теперь значение коэффициента c0(ρ). Рассмотрим при k = 0 диф- ференциальное уравнение(ρ2c0 0(ρ))0= 0⇔ρ2c0 0(ρ) = A0Тогда c0 0(ρ) =A0ρ2, c0(ρ) =ZA0ρ2dr = −A0ρ+ B0. Используя краевые усло- вия, получим A0= 0, B0=T0 2. Значит,c0(ρ) =T0 2Все коэффициенты найдены. Искомая функция u(ρ,θ) представима в виде ряда Фурье:u(ρ,θ) =T0 2+T0 2+∞Xk=1k + 2k + 1Pk−1(0)ρRkPk(cosθ).Известно [5], чтоPk(0) =(−1)m4m(2m)!(m!)2,k = 2m,0,k = 2m + 1.Используя это свойство многочленов Лежандра, искомую функцию можно преобразовать к виду u(ρ,θ) =T0 2+T0 2+∞Xm=0(−1)m(2m)!(2m + 3)4m(m!)2(2m + 2)ρR2m+1P2m+1(cosθ).Вычислим теперь температуру точек шара, расположенных вдоль радиусаθ= 0. Для этого подставим в последнюю формулуθ= 0. Учитывая, чтоPk(1) = 1, получим искомую температуру u(r, 0) =T0 2+T0 2+∞Xm=0(−1)m(2m)!(2m + 3)4m(m!)2(2m + 2)ρR2m+1(4.8)115 4.7. Задача на собственные значения для оператора ЛапласаВ некоторых случаях краевую задачу для уравнения Пуассона или уравнения Лапласа удается свести к задаче с однородными краевыми усло- виями на всей границе. Решение такой задачи целесообразно искать в виде ряда Фурье, разложив искомую функцию по собственным функциям опе- ратора Лапласа.Пусть Ω – некоторая ограниченная область, Γ – граница этой области.Определение 4.2. Числа λ, для которых краевая задача−∆u = λu,R∂u∂n+ SuΓ= 0(|R| + |S| 6= 0)имеет ненулевые решения u = u(M ) в области Ω, называются собствен- ными числами оператора Лапласа. При этом ненулевые решения u на- зываются собственными функциями оператора Лапласа, соответству- ющими числам λ.Задача нахождения собственных чисел и соответствующих им соб- ственных функций называется задачей на собственные значения для опе- ратора Лапласа.Однородное краевое условие R∂u∂n+ SuΓ= 0 может быть первого,второго или третьего рода. Также можно задать смешанное однородное краевое условие. Задача на собственные значения оператора Лапласа яв- ляется многомерным аналогом задачи Штурма–Лиувилля.Пусть функции f = f (M ) и g = g(M ) определены в области Ω. Введем скалярное произведение этих функций по правилу(f, g) =Z Z ZΩf gρdV,гдеρ=ρ(M ) – непрерывная на множестве Ω функция,ρ(M ) ≥ρ0> 0.Как и ранее,ρ(M ) называется весовой функцией.Норму, порожденную скалярным произведением, определим по фор- муле kf (M )k =v uu tZ Z ZΩf2ρdV .Множество функций, для которых норма конечна (kf (M )k < +∞),образуют пространство L2[Ω;ρ(M )]. Еслиρ(M ) = 1, то пространство обо- значим L2[Ω].Сформулируем без доказательства некоторые свойства собственных чисел и собственных функций оператора Лапласа:116 1) собственные числа оператора Лапласа вещественные;2) множество собственных чисел счетно, т. е. представляет собой по- следовательность {λn}+∞n=1;3) числа λn≥ 0 и lim n→+∞λn= +∞;4) собственные функции оператора Лапласа, соответствующие различ- ным собственным числам, ортогональны в пространстве L2[Ω;ρ(M )];5) собственные функции оператора Лапласа образуют полную, замкну- тую в пространстве L2[Ω;ρ(M )] систему функций. Если f ∈ L2[Ω;ρ(M )],то функция f может быть разложена в ряд Фурье по системе собствен- ных функций оператора и этот ряд будет сходиться к функции в норме пространства L2[Ω;ρ(M )].Рассмотрим несколько задач на собственные значения оператора Ла- пласа для разных областей и разных типов краевых условий.Пример 4.5. Найти собственные числа и собственные функции опе- ратора Лапласа в прямоугольнике Ω = {(x, y) : 0 ≤ x ≤ A, 0 ≤ y ≤ B},если на границе области функция u удовлетворяет условиям Дирихле.Решим задачу− ∂2u∂x2+∂2u∂y2= λu(x, y),0 < x < A, 0 < y < B,u(0, y) = 0,u(A, y) = 0,u(x, 0) = 0,u(x, B) = 0.Искомую функцию представим в виде произведения двух функций u(x, y) = X(x)Y (y), предполагая, чтоX(0) = 0,X(A) = 0,Y (0) = 0,Y (B) = 0.Тогда u(x, y) будет удовлетворять заданным краевым условиям. Под- ставим u(x, y) = X(x)Y (y) в дифференциальное уравнение:−(X00(x)Y (y) + X(x)Y00(y)) = λX(x)Y (y).Разделим обе части этого равенства на X(x)Y (y):− X00(x)X(x)+Y00(y)Y (y)= λ.Это равенство выполняется только в том случае, если−X00(x)X(x)= λ1и−Y00(y)Y (y)= λ2,где λ1и λ2– константы. При этом λ = λ1+ λ2 117 Значит, функции X(x) и Y (y) являются решениями двух задач Штурма–Лиувилля:(−X00= λ1X,0 < x < A,X(0) = 0,X(A) = 0и(−Y00= λ2Y,0 < y < B,Y (0) = 0,Y (B) = 0.Эти задачи решены в 1.3:λ1,k=πkA2,Xk(x) = sinπkxA,k = 1, 2, ...;λ2,m=πmB2,Ym(y) = sinπmyB,m = 1, 2, ... .Тогда числа λkm= λ1,k+ λ2,m и функции Ukm(x, y) = sinπkxAsinπmyB,k = 1, 2, ...;m = 1, 2, ..., будут собственными числами и собственными функциями оператора Лапласа для задачи Дирихле в прямоугольнике Ω.Найденные собственные функции образуют ортогональную, полную в про- странстве L2[Ω] систему функций.Пример 4.6. Найти собственные числа и собственные функции опе- ратора Лапласа в круге радиуса R (Ω = {(ρ,ϕ) : 0 <ρ≤ R, 0 ≤ϕ< 2π}),если на границе кругаρ= R функция u удовлетворяет условиям Дирихле.Запишем оператор Лапласа в полярной системе координат и решим следующую задачу:− 1ρ∂u∂ρρ∂u∂ρ+1ρ2∂2u∂ϕ2= λu,0 <ρ≤ R,0 ≤ϕ< 2π;u(ρ,ϕ) ограничена приρ→ 0 + 0,u(R,ϕ) = 0;u(ρ, 0) = u(ρ, 2π),∂u(ρ, 0)∂ϕ=∂u(ρ, 2π)∂ϕПо переменнойϕзадаются естественные для задачи периодические краевые условия.Функцию u, как и в предыдущем примере, будем искать в виде про- изведения двух функций u(ρ,ϕ) = Ψ(ρ)Φ(ϕ), считая, что выполняются условия:Ψ(ρ) ограничена приρ→ 0 + 0,Ψ(R) = 0,Φ(0) = Φ(2π),Φ0(0) = Φ0(2π).Тогда краевые условия для функции u будут выполнены. Подставим u(ρ,ϕ) == Ψ(ρ)Φ(ϕ) в дифференциальное уравнение:− 1ρ(ρΨ0(ρ))0Φ(ϕ) +1ρ2Ψ(ρ)Φ00(ϕ)= λΨ(ρ)Φ(ϕ).118 Пусть функция Φ(ϕ) является решением задачи Штурма–Лиувилля:−Φ00(ϕ) =νΦ(ϕ),0 ≤ϕ< 2π,Φ(ϕ+ 2π) = Φ(ϕ),Φ0(ϕ+ 2π) = Φ0(ϕ).Тогда функция Ψ(ρ) есть решение другой задачи Штурма–Лиувилля:−1ρ(ρΨ0(ρ))0+νρ2Ψ(ρ) = λΨ(ρ),0 <ρ≤ R,Ψ(ρ) ограничена приρ→ 0 + 0,Ψ(R) = 0.Решение задачи для функции Φ(ϕ) было получено в пособии [1] (при- мер 2.4). Собственные числаνk= k2, k = 0, 1, ... . Собственному числуν0= 0 соответствует одна собственная функция Φ0(ϕ) = 1, а каждо- му собственному числуνk= k2, k = 1, 2, ... – две собственные функцииΦk,1(ϕ) = cos(kϕ), Φk,2(ϕ) = sin(kϕ).Каждому числуνk будет соответствовать функция Ψk(ρ). Для функ- ции Ψk(ρ) задача Штурма–Лиувилля примет вид−1ρ(ρΨ0k(ρ))0+k2ρ2Ψk(ρ) = λkΨk(ρ),0 <ρ≤ R,Ψk(ρ) ограничена приρ→ 0 + 0,Ψk(R) = 0(k = 0, 1, ...).Оператор в правой части уравнения есть оператор Бесселя. В 1.6 най- дено, что ограниченным приρ→ 0 + 0 решением уравнения будет функ- ция Бесселя Ψk(ρ) = Jk(√λkρ). Подставив эту функцию в краевое условиеΨk(R) = 0, получим:Jk(pλkR) = 0.Обозначимγ=√λkR и запишем уравнение в виде Jk(γ) = 0. Это уравнение имеет множество простых решенийγk,m(m = 1, 2, ...). Числаγk,m– это корни функции Бесселя Jk(γ) (см. 1.5). Учитывая это, получим множество собственных чиселλk,m=γk,mR2,которым соответствуют функции Ψk,m(ρ) = Jkγk,mρR(k = 0, 1, ...,m == 1, 2, ...).Полученные числа λk,m являются также собственными числами опера- тора Лапласа. Соответствующие им собственные функции, как отмечалось ранее, будем искать в виде u(ρ,ϕ) = Ψ(ρ)Φ(ϕ). Используя полученные решения задач Штурма–Лиувилля найдем эти функции.119 Собственным числам оператора Лапласа λ0,m=γ0,mR2(m = 1, 2, ...)соответствуют собственные функции Ψ0,m(ρ) = J0γ0,mρRКаждому собственному числу λk,m=γk,mR2(k = 1, 2, ...,m == 1, 2, ...) соответствуют две собственные функцииΨk,m(ρ,ϕ) = cos(kϕ)Jkγk,mρR,Φk,m(ρ,ϕ) = sin(kϕ)Jkγk,mρRЗдесьγk,m(k = 0, 1, ...,m = 1, 2, ...) – корни уравнения Jk(γ) = 0.Найденные собственные функции оператора Лапласа образуют орто- гональную, полную в пространстве L2[Ω;ρ] систему функций.Пример 4.7. Найти собственные числа и собственные функции опе- ратора Лапласа на сфере радиуса R = 1 (Ω = {(θ,ϕ) : 0 ≤θ≤π, 0 ≤ϕ<< 2π}).Запишем оператор Лапласа в сферической системе координат и поло- жимρ= 1:∆Y =1sinθ∂∂θsinθ∂Y∂θ+1sin2θ∂2Y∂ϕ2Полученный оператор называется оператором Бельтрами. Поставим зада- чу на собственные значения этого оператора:−1sinθ∂∂θsinθ∂Y∂θ+1sin2θ∂2Y∂ϕ2= λY,0 <θ<π,0 <ϕ< 2π;Y (θ,ϕ) ограничена приθ→ 0 + 0 и приθ→π− 0;Y (θ, 0) = Y (θ, 2π),∂Y (θ, 0)∂ϕ=∂Y (θ, 2π)∂ϕПо обеим переменным заданы естественные краевые условия, следующие из свойств сферических координат.Решение будем искать в виде произведения двух функций Y (ρ,ϕ) == Ψ(θ)Φ(ϕ), предположив, что выполняются условия:Ψ(θ) ограничена приρ→ 0 + 0 и приθ→π− 0,Φ(0) = Φ(2π),Φ0(0) = Φ0(2π).Тогда для функции Y будут выполнены краевые условия исходной задачи.Подставим Y (θ,ϕ) = Ψ(θ)Φ(ϕ) в дифференциальное уравнение:−1sinθ(sinθ· Ψ0(θ))0Φ(ϕ) +1sin2θΨ(θ)Φ00(ϕ)= λΨ(θ)Φ(ϕ).120 Как и в предыдущем примере, для функции Φ(ϕ) решим задачу Штурма–Лиувилля:−Φ00(ϕ) =νΦ(ϕ),0 ≤ϕ< 2π,Φ(ϕ+ 2π) = Φ(ϕ),Φ0(ϕ+ 2π) = Φ0(ϕ).Тогда функция Ψ(θ) – решение другой задачи Штурма–Лиувилля:−1sinθ(sinθ· Ψ0(θ))0+νsin2θΨ(θ) = λΨ(θ),Ψ(θ) ограничена приρ→ 0 + 0 и приθ→π− 0.Первая задача Штурма–Лиувилля рассмотрена в предыдущем приме- ре. Собственные числаνk= k2, k = 0, 1, ... . Собственному числуν0= 0 со- ответствует одна собственная функция Φ0(ϕ) = 1, а каждому собственному числуνk= k2, k = 1, 2, ... – две собственные функции Φk,1(ϕ) = cos(kϕ),Φk,2(ϕ) = sin(kϕ).Каждому числуνk будет соответствовать функция Ψk(θ). При этом задача Штурма–Лиувилля для функции Ψk(θ) примет вид−1sinθ(sinθ· Ψ0k(θ))0+k2sin2θΨk(θ) = λΨk(θ),(4.9)Ψk(θ) ограничена приρ→ 0 + 0 и приθ→π− 0.Преобразуем уравнение относительно функции Ψk(θ) к виду−1sinθd dθ− sin2θ(−1sinθ)dΨk(θ)dθ+k2sin2θΨk(θ) = λΨk(θ).Введем новую переменную t = cosθ(−1 < t < 1). Будет выполняться цепочка равенствΨk(θ) = Ψk(arccos t) = v k(t) = v k(cosθ)(4.10)и v0k(t) = Ψ0k(θ)−1√1 − t2= −1sinθΨ0k(θ). Учитывая эту взаимосвязь между производными функций v k(t) и Ψk(θ), а также равенство sin2θ= 1 − t2,получим уравнение для функции v k(t)−d dt(1 − t2)dv kdt+k2 1 − t2v k= λv kПри этом функция v k(t) будет удовлетворять однородным краевым условиям:v k(t) ограничена при t → −1 + 0 и при t → 1 − 0.121 Получилась задача Штурма–Лиувилля для оператораLk(v) = − (1 − t2)v00+k2 1 − t2v.Эта задача была рассмотрена в 1.8. Собственные числа оператораLk(v) – это числа λn= n(n + 1) (n = k, k + 1, ...). Соответствующие им собственные функции оператора Lk(v) – это присоединенные функции Ле- жандра:v k,n(t) = Pk n(t) = (1 − t2)k/2d kdt kPn(t),где Pn(t) – многочлены Лежандра.Возвращаясь к переменнойθ, учитывая равенства (4.10), получим ре- шения задачи Штурма–Лиувилля (4.9). Искомые собственные числа – это числа λn= n(n + 1) (n = k, k + 1, ...). Им соответствуют собственные функции Ψk,n(θ) = Pk n(cosθ).Очевидно, что найденные собственные числа λn являются также соб- ственными числами оператора Бельтрами. Ранее указывалось, что соб- ственные функции этого оператора будем искать в виде Y (ρ,ϕ) = Ψ(θ)Φ(ϕ).Функции Ψ(θ) и Φ(ϕ) найдены. Учитывая эти функции, заключаем, что каждому собственному числу λn= n(n + 1) (n = 0, 1, ...) соответствует2n + 1 собственная функция вида:Y0n(θ,ϕ) = P0n(cosθ) = Pn(cosθ),Yk n(θ,ϕ) = Pk n(cosθ) sin(kϕ),k = 1, ..., n,Y−k n(θ,ϕ) = Pk n(cosθ) cos(kϕ),k = 1, ..., n.Найденные собственные функции оператора Бельтрами называются сферическими функциями. Они образуют ортогональную, полную в про- странстве L2[Ω; sinθ] систему функций. Доказательство приводится, на- пример, в [8].4.8. Решение задач разложением по собственным функциям оператора ЛапласаРассмотрим несколько примеров, в которых целесообразно искать ре- шение в виде ряда Фурье, разложив искомую функцию по собственным функциям оператора Лапласа.Пример 4.8. Найти стационарное распределение температуры в длин- ном металлическом стержне с прямоугольным сечением (0 ≤ x ≤ A,0 ≤ y ≤ B) при условии, что в стержне выделяется тепло с объемной плот- ностью распределения тепловой энергии Q, а его грани поддерживаются122 при температуре T0. Температура в каждом поперечном сечении стержня считается одинаковой.Функция u(x, y), описывающая стационарное распределение темпера- туры в поперечном сечении стержня, не зависит от времени (∂u∂t= 0),поэтому она удовлетворяет уравнению∂2u∂x2+∂2u∂y2= −Q0(0 < x < A,0 < y < B),краевым условиям:u(0, y) = T0,u(A, y) = T0,u(x, 0) = T0,u(x, B) = T0,где Q0=QK, K – коэффициент теплопроводности.Искомая функция u(x, y) является решением задачи Дирихле для урав- нения Пуассона в прямоугольнике Ω = {(x, y) : 0 ≤ x ≤ A, 0 ≤ y ≤ B}:∆u = −Q0,uΓ= T0Сначала сведем эту краевую задачу к задаче с однородными краевыми условиями. Используем замену u(x, y) = v(x, y) + w(x, y). Очевидно, что функция w(x, y) = T0удовлетворяет тем же краевым условиям, что и функция u(x, y). Тогда функция v(x, y) будет удовлетворять однородным краевым условиям. Учитывая, что u(x, y) = v(x, y) + T0, запишем краевую задачу относительно функции v(x, y):∂2v∂x2+∂2v∂y2= −Q0(0 < x < A,0 < y < B),v(0, y) = 0,v(A, y) = 0,v(x, 0) = 0,v(x, B) = 0.(4.11)Функция v(x, y) удовлетворяет однородным краевым условиям на всей границе.Полученную краевую задачу будем решать методом Фурье, выполнив разложение в ряд по системе собственных функций оператора Лапласа.Задача на собственные значения оператора Лапласа в прямоугольнике при условии, что функция на границе удовлетворяет однородным краевым условиям Дирихле, была решена в примере 4.5. Собственные числа опера- тора:λkm=πkA2+πmB2,соответствующие им собственные функции:Vkm(x, y) = sinπkxAsinπmyB,k = 1, 2, ...;m = 1, 2, ... .123 Собственные функции образуют ортогональную, полную в пространствеL2[Ω] систему функций.Функцию v(x, y) представим в виде ряда Фурье по найденной системе собственных функций оператора Лапласа:v(x, y) =+∞Xk=1+∞Xm=1c kmVkm(x, y).Функцию Q0разложим в ряд Фурье по той же системе функций:Q0=+∞Xk=1+∞Xm=1q kmVkm(x, y),где q km=(Q0, Vkm(x, y))kVkm(x, y)k2. Вычислим скалярные произведения:(Q0, Vkm(x, y)) =AZ0BZ0Q0Vkm(x, y)dx dy = Q0AZ0sinπkxAdxBZ0sinπmyBdy == Q0AB((−1)k− 1)((−1)m− 1)π2km,kVkm(x, y)k2=AZ0BZ0V2km(x, y)dx dy =AZ0sin2πkxAdxBZ0sin2πmyBdy =AB4Тогда qkm=4Q0((−1)k− 1)((−1)m− 1)π2kmПодставим функции v(x, y) и Q0, представленные в виде рядов, в урав- нение краевой задачи (4.11):+∞Xk=1+∞Xm=1c km ∂2Vkm∂x2+∂2Vkm∂y2= −+∞Xk=1+∞Xm=1q kmVkm(x, y).Поскольку разложение в ряд Фурье единственно и для функций Vkm(x, y)выполняются равенства− ∂2Vkm∂x2+∂2Vkm∂y2= λkmVkm,то c km=q kmλkm для всех k = 1, 2, ...;m = 1, 2, ... и решение краевой задачи(4.11) представляется в виде v(x, y) =+∞Xk=1+∞Xm=1q kmλkmVkm(x, y) =124 =4Q0π2+∞Xk=1+∞Xm=1((−1)k− 1)((−1)m− 1)km (πkA)2+ (πmB)2sinπkxAsinπmyBУчитывая, что u(x, y) = v(x, y) + T0, получим решение исходной задачи u(x, y) = T0+4Q0π2+∞Xk=1+∞Xm=1((−1)k− 1)((−1)m− 1)km (πkA)2+ (πmB)2sinπkxAsinπmyBПример 4.9. Найти стационарное распределение температуры в длин- ном металлическом стержне с круглым сечением (0 ≤ρ≤ R,0 ≤ϕ≤≤ 2π), если в стержне выделяется тепло с объемной плотностью распреде- ления тепловой энергии Q, а поверхность стержняρ= R поддерживается при температуре T0. Температура в каждом поперечном сечении стержня считается одинаковой.Пусть функция u описывает стационарное распределение температуры в поперечном сечении стержня. Как и в предыдущем примере, она является решением задачи Дирихле для уравнения Пуассона∆u = −Q0,uΓ= T0Поскольку в этом примере сечение стержня – круг, то в модели исполь- зуем полярную систему координат. Функция u = u(ρ,ϕ) удовлетворяет стационарному уравнению теплопроводности1ρ∂∂ρρ∂u∂ρ+1ρ2∂2u∂ϕ2= −Q0(0 <ρ< R,0 ≤ϕ< 2π)и краевым условиям:u(ρ,ϕ) ограничена приρ→ 0 + 0,u(R,ϕ) = T0,u(ρ, 0) = u(ρ, 2π),∂u(ρ, 0)∂ϕ=∂u(ρ, 2π)∂ϕЗдесь Q0=QK, K – коэффициент теплопроводности.Сведем эту краевую задачу к задаче с однородными краевыми усло- виями. Аналогично предыдущему примеру, используем замену u(ρ,ϕ) == v(ρ,ϕ) + T0. Функция v(ρ,ϕ) будет удовлетворять однородным краевым условиям и являться решением следующей краевой задачи:1ρ∂∂ρρ∂v∂ρ+1ρ2∂2v∂ϕ2= −Q0(0 <ρ< R,0 ≤ϕ< 2π),(4.12)краевые условия:v(ρ,ϕ) ограничена приρ→ 0 + 0,v(R,ϕ) = 0,125 v(ρ, 0) = v(ρ, 2π),∂v(ρ, 0)∂ϕ=∂v(ρ, 2π)∂ϕРешение этой задачи будем искать в виде ряда Фурье по системе соб- ственных функций оператора Лапласа в круге.Задача на собственные значения оператора Лапласа в круге для кра- евых условий Дирихле была решена в примере 4.6. Собственным числам оператора Лапласа λ0m=γ0mR2(m = 1, 2, ...) соответствуют собственные функции Ψ0m(ρ) = J0γ0mρRКаждому собственному числу λkm=γkmR2(k = 1, 2, ...,m == 1, 2, ...) соответствуют две собственные функцииΨkm= cos(kϕ)JkγkmρR,Φkm= sin(kϕ)JkγkmρRЗдесьγkm(k = 0, 1, ...,m = 1, 2, ...) – корни уравнения Jk(γ) = 0.Найденные собственные функции оператора Лапласа образуют орто- гональную, полную в пространстве L2[Ω;ρ] систему функций.Представим функцию v(ρ,ϕ) в виде ряда Фурье по найденной системе собственных функций оператора Лапласа:v(ρ,ϕ) =+∞Xm=1a0mΨ0m(ρ,ϕ) ++∞Xk=1+∞Xm=1(a kmΨkm(ρ,ϕ) + b kmΦkm(ρ,ϕ)).Разложим функцию Q0в ряд по той же системе функций:Q0=+∞Xm=1α0mΨ0m(ρ,ϕ) ++∞Xk=1+∞Xm=1(αkmΨkm(ρ,ϕ) +βkmΦkm(ρ,ϕ)),гдеα0m=(Q0, Ψ0m(ρ,ϕ))kΨ0m(ρ,ϕ)k2,αkm=(Q0, Ψkm(ρ,ϕ))kΨkm(ρ,ϕ)k2,βkm=(Q0, Φkm(ρ,ϕ))kΦkm(ρ,ϕ)k2Вычислим(Q0, Ψkm(ρ,ϕ)) =2πZ0RZ0Q0Ψkm(ρ,ϕ)ρdϕdρ== Q0 2πZ0cos(kϕ)dϕRZ0JkγkmρRρdρПоскольку2πZ0cos(kϕ)dϕ= 0, то (Q0, Ψkm(ρ,ϕ)) = 0. Аналогично,126 учитывая, что2πZ0sin(kϕ)dϕ= 0, получим (Q0, Φkm(ρ,ϕ)) = 0. Значит,αkm= 0,βkm= 0(k = 1, 2, ...,m = 1, 2, ...).Найдемα0m. Используя формулы (1.28), (1.29), свойство J0 0(x) = −J1(x)и равенство J0(γ0m) = 0 получим:kΨ0m(ρ,ϕ)k2=2πZ0RZ0Ψ2 0m(ρ,ϕ)ρdϕdρ=2πZ0dϕRZ0J2 0γ0mρRρdρ== 2πR2 2J2 1(γ0m),(Q0, Ψ0m(ρ,ϕ)) =2πZ0RZ0Q0Ψ0m(ρ,ϕ)ρdϕdρ= Q0 2πZ0dϕRZ0J0γ0mρRρdρ== Q0 2πR2γ0mJ1(γ0m).Тогдаα0m=2Q0γ0mJ1(γ0m)Подставим функции v(x, y) и Q0, представленные в виде рядов, в урав- нение (4.12) краевой задачи:+∞Xm=1a0m∆Ψ0m(ρ,ϕ) ++∞Xk=1+∞Xm=1(a km∆Ψkm(ρ,ϕ) + b km∆Φkm(ρ,ϕ)) == −+∞Xm=1α0mΨ0m(ρ,ϕ) −+∞Xk=1+∞Xm=1(αkmΨkm(ρ,ϕ) +βkmΦkm(ρ,ϕ)).Поскольку разложение в ряд Фурье единственно и для функций Ψ0m(ρ,ϕ),Ψkm(ρ,ϕ) и Φkm(ρ,ϕ) выполняются равенства−∆Ψ0m= λkmΨ0m,−∆Ψkm= λkmΨkm,−∆Φkm= λkmΦkm,то a0m=α0mλ0m, a km=αkmλkm, b km=βkmλkm для всех k = 1, 2, ...;m = 1, 2, ... и решение краевой задачи представляется в виде v(ρ,ϕ) =+∞Xm=1α0mλ0mΨ0m(ρ,ϕ) =+∞Xm=1 2Q0R2γ3 0mJ1(γ0m)J0γ0mρR127 Учитывая, что u(x, y) = v(x, y)+T0, получим решение исходной задачи u(x, y) = T0+ 2Q0R2+∞Xm=1 1γ3 0mJ1(γ0m)J0γ0mρR4.9. Метод функции ГринаВ ряде случаев требуется найти аналитическое решение краевой зада- чи не в виде ряда, а в виде интеграла. Одним из методов нахождения та- кого решения является метод функции Грина. Идея метода состоит в том,что сначала внутри области Ω строится некоторая гармоническая функцияG1. Затем с помощью этой функции находится решение заданной краевой задачи.Пусть Ω ∈ R3– ограниченная область, M0– некоторая фиксированная точка, лежащая внутри Ω. Обозначим через G1(M0, M ) функцию, удовле- творяющую следующим условиям:1) G1(M0, M ) непрерывна на множестве Ω ∪ Γ;2) G1(M0, M ) является гармонической функцией внутри области Ω: ∆G1== 0;3) G1(M0, M ) на границе Γ удовлетворяет условию: G1(M0, M )Γ=1rΓ,где r = |M0M | =p(x − x0)2+ (y − y0)2+ (z − z0)2Если область Ω имеет гладкую или кусочно-гладкую границу, то такая функция существует. Из утверждения 4.6 следует, что функция G1(M0, M )единственна.Прежде чем привести примеры построения функции G1(M0, M ), по- кажем как с ее помощью можно найти значения другой гармонической функции в области Ω, если известны значения этой функции на границе Γ.Пусть u = u(M ) – гармоническая в области Ω функция, непрерывная на множестве ΩS Γ. Тогда для любой внутренней точки M0(x0, y0, z0) мно- жества Ω справедлива основная формула теории гармонических функций(4.5):u(M0) =1 4πZ ZΓ 1r∂u∂n− u∂∂n 1rds,(4.13)где r = |M0M | =p(x − x0)2+ (y − y0)2+ (z − z0)2Применим к функциям u и G1вторую формулу Грина (4.2), задав v = G1. Поскольку ∆u = 0 и ∆G1= 0, то0 =Z ZΓu∂G1∂n− G1∂u∂nds.128 Умножив обе части этого равенства на1 4πи сложив с (4.13), получим u(M0) =1 4πZ ZΓu ∂G1∂n−∂∂n 1r−∂u∂nG1−1rds.Из условия G1(M0, M )Γ=1rΓследует равенство u(M0) =1 4πZ ZΓu∂∂nG1−1rds.(4.14)Определение 4.3. Функция G(M0, M ) =1r− G1(M0, M ) называется функцией Грина в области Ω.Используя функцию Грина, формулу (4.14) можно записать в виде u(M0) = −1 4πZ ZΓu∂∂nG(M0, M )ds.(4.15)Если известны значения функции u на границе Γ и известна функцияГрина G(M0, M ) в области Ω, то формула (4.15) дает возможность вычис- лять значения гармонической функции u в любой внутренней точке M0области Ω.Отметим, что построение функции Грина не является тривиальной задачей, но для некоторых областей ее удается найти.Пример 4.10. Построить функцию Грина в шаре KR(O) радиуса R.Найдем сначала функцию G1(M0, M ) – непрерывное в шаре решение задачи Дирихле∆G1= 0,G1(M0, M )ΓR=1rΓR1. Пусть точка M0(x0, y0, z0) – центр шара. Тогда очевидно, что функ- ция G1(M0, M ) =1R. Действительно, эта функция-константа является непре- рывной и гармонической в шаре. Если M лежит на границе ΓR, то |M0M | == R и выполняется краевое условие G1(M0, M )ΓR=1RΓR2. Пусть теперь M0не является центром шара. На луче [OM0) найдем такую точку M1(x1, y1, z1), что |OM0| · |OM1| = R2. Точка M1находится вне шара.Рассмотрим произвольную точку M (x, y, z), находящуюся внутри или на поверхности шара. Соединим ее с точками M0и M1. Обозначим129 r = |M M0| =p(x − x0)2+ (y − y0)2+ (z − z0)2,r1= |M M1| =p(x − x1)2+ (y − y1)2+ (z − z1)2,а такжеρ= |OM0| (рис 4.1).1   2   3   4   5   6   7   8   9   10

MMOM0 1ρwrγ1rРис. 4.1Покажем, что функцияG1(M0, M ) =Rρ1r1является искомой функцией.Посколькуρи R являются постоян- ными, то функцияG1(M0, M ) ==Rρ1p(x − x1)2+ (y − y1)2+ (z − z1)2является гармонической в пространствеR3, исключая точку M1(см. пример 4.3).Значит, она является непрерывной и гармонической в шаре. Покажем, что на границе шара выполняется условие G1(M0, M )ΓR=1rΓRСоединим теперь точку M с точкой O (рис. 4.1). Точки O, M0, M1ле- жат на одной прямой, поэтому все 4 точки O, M0, M1и M лежат в одной плоскости. Обозначим w = |OM | иγ=6M OM1. Рассмотрим треуголь- ники OM0M и OM M1. Учитывая, что |OM0| · |OM1| = R2, по теореме косинусов получим r =p w2+ρ2− 2wρcosγ,r1=s w2+R2ρ2− 2wRρcosγ(4.16)Если точка M лежит на границе ΓR, то |OM | = w = R и справедливы равенства:rΓR=pR2+ρ2− 2Rρcosγ,r1ΓR=sR2+R2ρ2− 2R2ρcosγ=Rρpρ2+ R2− 2Rρcosγ=RρrΓRИз этого следует, что краевое условие G1(M0, M )ΓR=Rρ1r1ΓR=1rΓRвыполняется.Из теоремы единственности решения задачи Дирихле (утверждение4.6) следует, что других функций, удовлетворяющих тем же условиям, не существует.130 Тогда искомая функция Грина имеет видG(M0, M ) =1r−1R,M0– центр шара,1r−Rρ1r1,M0– не центр шара.Покажем, как найденная функция используется для решения задачи.Пример 4.11. Найти гармоническую в шаре радиуса R функцию,принимающую на границе шара ΓRзначения u0Пусть u = u(M ) – искомая функция. Область Ω = KR(O) – шар радиуса R с центром в точке O. Функция u является решением задачиДирихле:∆u = 0,uΓR= u0Требуется найти u(M0), где M0лежит внутри шара.1. Пусть точка M0(x0, y0, z0) – центр шара. Воспользуемся формулой(4.15) и подставим в нее найденную функцию Грина:u(M0) = −1 4πZ ZΓRu∂∂nG(M0, M )ds = −1 4πZ ZΓRu∂∂n 1r−1Rds.Поскольку M0– центр шара, то∂∂n 1r−1R=∂∂r 1r−1R= −1r2Тогда u(M0) =1 4πZ ZΓRu1r2ds =1 4πZ ZΓRu0 1R2ds.В итоге получили результат утверждения 4.4:u(M0) =1 4πR2Z ZΓRu0ds.(4.17)2. Пусть M0(x0, y0, z0) не является центром шара. Для нее применим формулу (4.15):u(M0) = −1 4πZ ZΓRu0∂∂n 1r−Rρ1r1ds.Поскольку w = |OM |, то нетрудно заметить, что∂∂n=∂∂w. Тогда,используя формулы (4.16), получим∂∂n 1r−Rρ1r1= −w −ρcosγ(p w2+ρ2− 2wρcosγ)3+131 +Rρw −R2ρcosγ s w2+R2ρ2− 2wRρcosγ!3== −w −ρcosγr3+Rρw −R2ρcosγr3 1Для точек, лежащих на поверхности ΓR, w = |OM | = R и выполняется равенство r1=Rρr, поэтому∂∂n 1r−Rρ1r1= −R −ρcosγr3+RρR −R2ρcosγR3ρ3r3== −R −ρ2Rr3= −R2−ρ2Rr3Тогда u(M0) =R2−ρ2 4πRZ ZΓRu0r3ds.(4.18)Этот интеграл называется интегралом Пуассона для шара.Учитывая, что для точек, расположенных на сфере ΓR, выполняется равенство r2= R2+ρ2− 2Rρcosγ, преобразуем полученный интеграл к виду u(M0) =R2−ρ2 4πRZ ZΓRu0(R2+ρ2− 2Rρcosγ)3 2ds.Перейдем теперь к сферическим координатам с началом в точке O:x = R sinθcosϕ,y = R sinθsinϕ,z = R cosθДля сферы ds = R2sinθdθdϕ, поэтому u(M0) =R(R2−ρ2)4π2πZ0dϕπZ0u0(θ,ϕ)sinθ(R2+ρ2− 2Rρcosγ)3 2dθ(4.19)132 Этот интеграл, как и интеграл (4.18), называется интегралом Пуассона для шара.Приведем пример использования полученной формулы.Пример 4.12. Верхняя половина поверхности шара радиуса R под- держивается при температуре T0, а нижняя – при температуре, равной ну- лю. Найти стационарное распределение температуры вдоль радиуса шараθ= 0.ρθOMM0Рис. 4.2В задаче удобно использовать сфери- ческую систему координат. Опишем шарKR(O) ∪ ΓR= {(ρ,θ,ϕ) : 0 ≤ρ≤≤ R,0 ≤ Θ ≤π,0 ≤ϕ< 2π}. Ис- комая функция u является решением за- дачи Дирихле:∆u = 0,uΓR= u0,где u0(θ) =T0,0 ≤θ≤π2,0,π2<θ≤π– функция,описывающая температуру на поверхно- сти шара.Согласно условиям задачи функция u не зависит от угловой коорди- натыϕЭта задача была решена в примере 4.4. Решение было получено в виде ряда Фурье (4.9). Найдем теперь функцию u другим способом.Пусть M0(ρ0,θ0,ϕ0) – произвольная точка шара. Найдем температуру шара в этой точке, используя интеграл Пуассона (4.19):u(M0) =T0R(R2−ρ2)4π2πZ0dϕπ2Z0sinθ(R2+ρ2− 2Rρcosγ)3 2dθ==T0R(R2−ρ2)2π2Z0sinθ(R2+ρ2− 2Rρcosγ)3 2dθВычислить температуру во всех точках шара с помощью этого инте- грала сложно.Рассмотрим точки M0(ρ, 0,ϕ), расположенные вдоль радиусаθ= 0(рис. 4.2). Пусть M – произвольная точка поверхности шара. Тогдаγ=133 =6M0OM =θ, следовательно,u(M0) =T0R(R2−ρ2)2π2Z0sinθ(R2+ρ2− 2Rρcosθ)3 2dθ==T0R(R2−ρ2)4Rρπ2Z0d(R2+ρ2− 2Rρcosθ)(R2+ρ2− 2Rρcosθ)3 2==T0(R2−ρ2)4ρ−2(R2+ρ2− 2Rρcosθ)1 2π2 0Подставив пределы интегрирования, получим u(M0) =T0 21 +Rρ−R2−ρ2ρpR2+ρ2!Эту формулу можно использовать для всех точек радиуса, исключая центр шара (ρ= 0). Температуру в центре шара найдем по формуле (4.17)u(M0) =1 4πR2Z ZΓRu0ds =T0 4πR2Z ZΓ1ds.Здесь Γ1– верхняя половина сферы. ПосколькуRRΓ1ds = 2πR2(площадь половины поверхности сферы), то u(M0) =T0 4πR2 2πR2=T0 24.10. Метод сеток для уравнения ПуассонаПокажем, как строится разностная схема для уравнения Пуассона в прямоугольнике Ω с краевыми условиями Дирихле на границе Γ. Пусть требуется найти функцию u(x, y), удовлетворяющую уравнению Пуассона∂2u∂x2+∂2u∂y2= f (x, y)(4.20)в области 0 < x < a, 0 < y < b и краевым условиям u(0, y) =ϕ1(y), u(a, y) =ϕ2(y),u(x, 0) =ψ1(x), u(x, b) =ψ2(x).(4.21)Покроем прямоугольник Ω = [0; a] × [0; b] сеткойωh xh y= {(x i, y j)}, i == 0, ..., n, j = 0, ..., m, разбив отрезок [0; a] на n равных частей точками134 xi= ih x, i = 0, ..., n, а отрезок [0; b] – на m равных частей точками y j= jh y,j = 0, ..., m (h x=a n, h y=b m– шаги сетки по осям 0x и 0y соответственно)(рис. 4.3).O xxxyyhh1i0yyx0a xnbym==jyx1Рис. 4.3Все функции, заданные в краевой задаче, будем рассматривать только в узлах сетки. Обозначим u ji= u(x i, y j).Для внутренних узлов сетки запишем уравнение Пуассона, заменив производные в узлах разностными отношениями (1.46):u ji−1− 2u ji+ u ji+1h2x+ O(h2x) +u j−1i− 2u ji+ u j+1i h2y+ O(h2y) = f (x i, y j),i = 1, ..., n − 1,j = 1, ..., m − 1.В граничных узлах функция u(x, y) будет удовлетворять условиям uj0=ϕ1(y j),u jn=ϕ2(y j),j = 0, ..., m,u0i=ψ1(x i),u mi=ψ2(x i),i = 0, ..., n.Отбросим в уравнениях погрешности аппроксимации производных. В ре- зультате получим уравнения относительно сеточной функции ˜u:˜u ji−1− 2˜u ji+ ˜u ji+1h2x+˜u j−1i− 2˜u ji+ ˜u j+1i h2y= f (x i, y j),i = 1, ..., n − 1, j = 1, ..., m − 1.(4.22)Для этой функции запишем краевые условия˜u j0=ϕ1(y j), ˜u jn=ϕ2(y j),j = 0, ..., m,˜u0i=ψ1(x i), ˜u mi=ψ2(x i),i = 0, ..., n.(4.23)135 i( x,yj +1)( xi−1(xi, yj)( xi+1( xi, yj −1), yj), yj)Рис. 4.4Уравнения (4.22) с условиями (4.23)называются разностной схемой для крае- вой задачи (4.20), (4.21).На рис. 4.4 изображен шаблон схе- мы, т. е. указаны узлы сетки, которые ис- пользовались при аппроксимации произ- водных в узле (x i, y j). Построенная схема для уравнения Пуассона является пяти- точечной разностной схемой.В уравнениях (4.22), (4.23) ˜u ji– при- ближенные значения функции u(x, y) в узлах сетки. Так как значения ˜u в граничных узлах известны, необходимо решить только уравнения (4.22).Умножая на h2x h2y, преобразуем их к виду h2x˜u j−1i+ h2y˜u ji−1− 2(h2x+ h2y)˜u ji+ h2y˜u ji+1++ h2x˜u j+1i= h2x h2y f (x i, y j),i = 1, ..., n − 1, j = 1, ..., m − 1.Получилась система, состоящая из (n−1)(m−1) уравнений с (n−1)(m−1)неизвестными ˜u ji. Каждое уравнение в системе связывает 5 неизвестных.Для решения систему удобно записать в матричном виде. Вид матрицы коэффициентов зависит от способа нумерации неизвестных.1 20 12xxxxyyyyyxi1 23n −1n +1 n +2n +3 2n −2n2−1n2n2+1n3−3m(n −1)y(−1)xn0 0=ambj=Рис. 4.5Если неизвестные ˜u ji, i = 1, ..., n − 1, j = 1, ..., m − 1, перенумеровать так, как это указано на рис. 4.5, то получится уравнениеA˜u = F ,(4.24)где ˜u = [˜u1 1...˜u1n−1˜u2 1...˜u2n−1...˜u m−1 1...˜u m−1n−1]t– матрица-столбец неизвестных.Число неизвестных в системе N = (n − 1)(m − 1); A – матрица коэффи- циентов, имеющая 5 отличных от нуля диагоналей. Структура матрицы136 показана на рис. 4.6. Крестиками отмечены отличные от нуля элементы матрицы. При указанном способе нумерации узлов сетки A – это ленточ- ная матрица. Ширина ленты определяется числом узлов сетки по оси 0x.Кроме того, A – это матрично-трехдиагональная матрица. Размер каждой подматрицы – (n − 1) × (n − 1). Число подматриц по вертикали и гори- зонтали – (m − 1). Естественно, что структура матрицы учитывается при выборе алгоритма решения системы; F – матрица-столбец правых частей.При i = 1, i = n − 1, j = 1, j = m − 1 в матрице F учитываются краевые условия.x x x x x x x x x x x x x xx xx x x x x xx xx x x x x x x xx xx x x x x x x x x x x x x x x x x xx xx xx xx x x x x x x x x x x x x xx xx xx xx xx xA =Рис. 4.6При решении прикладных задач обычно задают n и m ∼ 50–100. Пусть n=101, m=101, тогда число уравнений в системе (4.24)N = (n − 1)(m − 1) = 10 4. Для решения таких больших систем уравнений следует использовать алгоритмы, учитывающие разряженность матрицы коэффициентов A.Отметим некоторые методы, применяемые для решения систем урав- нений подобного вида: 1) прямые методы – метод Гаусса для ленточных матриц, метод матричной прогонки, продольно-поперечный метод; 2) ите- рационные методы – метод Зейделя, метод последовательной верхней ре- лаксации. В настоящее время имеются стандартные программы, реализу- ющие все перечисленные методы.Покажем, что система (4.24) однозначно разрешима. Если f (x, y) = 0,то для построенной разностной схемы справедлив принцип максимума.Утверждение 4.9. Если f (x i, y i) = 0, i = 1, ..., n − 1, j = 1, ..., m −137 − 1, то решение разностной задачи (4.22), (4.23) достигает наибольшего(наименьшего) значения в граничных узлах сетки.Доказательство. Предположим, напротив, что наибольшее значение функция ˜u принимает во внутреннем узле сетки. Пусть (x p, y q) – внутрен- ний узел, имеющий наибольшую абсциссу среди тех узлов сетки, в которых˜u достигает наибольшего значения. Для этого узла ˜u qp+1− ˜u qp< 0. Запишем для него уравнение (4.22), учитывая, что f (x i, y j) = 0:(˜u qp−1− ˜u qp) + (˜u qp+1− ˜u qp)h2x+(˜u q−1p− ˜u qp) + (˜u q+1p− ˜u qp)h2y= 0.Поскольку ˜u qp+1− ˜u qp< 0 и ˜u qp−1− ˜u qp≤ 0, ˜u q−1p− ˜u qp≤ 0, ˜u q+1p− ˜u qp≤ 0.Последнее равенство невозможно, получили противоречие.Из принципа максимума следует, что разностная задача (4.22), (4.23),равносильная уравнению A˜u = 0 при условии, что f (x, y) = 0,ϕ1(y) = 0,ϕ2(y) = 0,ψ1(x) = 0,ψ2(x) = 0, имеет единственное нулевое решение.Значит, определитель матрицы A отличен от нуля и система (4.24) имеет единственное решение.Построенная разностная схема устойчива по правой части, т. е. для функции ˜u в сеточной норме справедливо неравенство k˜uk hx hy≤ C kf k hx hy+ kϕ1k hy+ kϕ2k hy+ kψ1k hx+ kψ2k hx,где C не зависит от h xи h y. Доказательство приводится, например, в [8].Очевидно, что эта разностная схема аппроксимирует краевую задачу с по- грешностью O(h2x) + O(h2y).Из устойчивости и аппроксимации разностной схемы следует ее сходи- мость. Доказательство выполняется аналогично доказательству сходимо- сти явной разностной схемы.4.11. Классификация уравнений второго порядкаЗапишем в самом общем виде линейное дифференциальное уравнение второго порядка nXk=1nXi=1a ki∂2u∂x k∂x i+nXi=1b i∂u∂x i+ cu = f,(4.25)где a ki, b i, c, f – функции n независимых переменных x1, ..., x n, заданных в некоторой области Ω ⊂ Rn. Считаем, что a ki= a ik. Зафиксируем некото- рую точку M0(x0 1...x0n) и рассмотрим матрицу коэффициентов при вторых138 производныхA0= A(x0 1, ..., x0n) =a0 11... a0 1n a0n1... a0nnОбозначим через λ1, λ2, ..., λn собственные числа матрицы A0. Так как A0симметрична, все собственные числа действительны.Определение 4.4. Если все собственные числа λ1, λ2, ..., λn имеют одинаковые знаки, то уравнение (4.25) называется уравнением эллипти- ческого типа в точке M0. Если λ1, ..., λn одинакового знака во всех точках области Ω, то уравнение называется уравнением эллиптического типа в области Ω.Определение 4.5. Если одно собственное число имеет один знак, а все остальные – другой, то уравнение называется уравнением гиперболи- ческого типа в точке M0(или во всей области).Определение 4.6. Если одно собственное число, например λ1, равно нулю, а все остальные одного знака, то уравнение называется уравнением параболического типа в точке M0(или в области Ω).Наконец, если в одних точках области Ω уравнение принадлежит к одному типу, а в других – к другому, то (4.25) называется уравнением смешанного типа.Рассмотрим согласно данной классификации уравнения, изученные ранее. Поскольку использовались уравнения с постоянными коэффициен- тами и матрица A с ее собственными числами не менялись от точки к точке,все уравнения имели определенный тип в области.Пример 4.13. Рассмотрим уравнение Лапласа∂2u∂x2+∂2u∂y2+∂2u∂z2= 0.Матрица A имеет видA =1 0 0 0 1 0 0 0 1,собственные числа матрицы λ1= λ2= λ3= 1. Следовательно, уравнениеЛапласа – это уравнение эллиптического типа.Пример 4.14. Волновое уравнение в двумерном случае имеет вид∂2u∂t2= a2 ∂2u∂x2+∂2u∂y2+ f (x, y).139 Матрица коэффициентов при вторых производныхA =1 00 0 −a2 00 0−a2Собственные числа λ1= 1, λ2= λ3= −a2. Уравнение имеет гиперболиче- ский тип.Пример 4.15. Рассмотрим уравнение теплопроводности∂u∂t= a2 ∂2u∂x2+∂2u∂y2+∂2u∂z2+ f (x, y, z, t).В этом уравнении можно считать, что коэффициент при∂2u∂t2равен нулю.Матрица A будетA =0 00 00 −a2 00 00−a2 00 00−a2Собственные числа матрицы A: λ1= 0, λ2= λ3= λ4= −a2. Уравнение теплопроводности – это уравнение параболического типа.Из разобранных примеров видно, что наиболее просто тип уравнения определяется в случае, когда матрица A диагональная, т. е. когда уравне- ние не содержит смешанных производных. Такой вид уравнения принято называть каноническим.Уравнение с постоянными коэффициентами всегда можно привести к каноническому виду с помощью линейной замены переменных, например используя ту же методику, что и для приведения квадратичной формы к каноническому виду. Проиллюстрируем эту методику на примере.Пример 4.16. Привести уравнение2∂2u∂x2+ 3∂2u∂y2+ 4∂2u∂z2− 4∂2u∂x∂y− 4∂2u∂y∂z− 2∂u∂x− 2∂u∂y−∂u∂z= f (x, y, z) (4.26)к каноническому виду.Выпишем симметричную матрицу коэффициентов при вторых произ- водных:A =2 −2 0−2 3 −2 0 −2 4140 Нетрудно проверить, что собственными числами матрицы A будут λ1= 0,λ2= 3, λ3= 6, а соответствующие нормированные собственные векторыe1=2/3 2/3 1/3,e2=−2/3 1/3 2/3,e3=1/3−2/3 2/3Матрица перехода, составленная из собственных векторов, и даст необ- ходимое преобразование координатT =2/3 −2/3 1/3 2/3 1/3−2/3 1/3 2/3 2/3Введем новые координаты x0, y0, z0следующим образом:x0y0z0= T−1x yz,или, учитывая, что T−1= Tt,x0=2 3x +2 3y +1 3z,y0= −2 3x +1 3y +2 3z,z0=1 3x −2 3y +2 3z.По правилу дифференцирования сложной функции получим∂∂x=2 3∂∂x0−2 3∂∂y0+1 3∂∂z0,∂∂y=2 3∂∂x0+1 3∂∂y0−2 3∂∂z0,∂∂z=1 3∂∂x0+2 3∂∂y0+2 3∂∂z0Перепишем уравнение (4.26) в новых координатах2 4 9∂2u∂x02+4 9∂2u∂y02+1 9∂2u∂z02−8 9∂2u∂x0∂y0+4 9∂2u∂x0∂z0−4 9∂2u∂y0∂z0++ 3 4 9∂2u∂x02+1 9∂2u∂y02+4 9∂2u∂z02+4 9∂2u∂x0∂y0−8 9∂2u∂x0∂z0−4 9∂2u∂y0∂z0+141 + 4 1 9∂2u∂x02+4 9∂2u∂y02+4 9∂2u∂z02+4 9∂2u∂x0∂y0+4 9∂2u∂x0∂z0+8 9∂2u∂y0∂z0−− 4 4 9∂2u∂x02−2 9∂2u∂y02−2 9∂2u∂z02−2 9∂2u∂x0∂y0−2 9∂2u∂x0∂z0+5 9∂2u∂y0∂z0−− 4 2 9∂2u∂x02+2 9∂2u∂y02−4 9∂2u∂z02+5 9∂2u∂x0∂y0+2 9∂2u∂x0∂z0−2 9∂2u∂y0∂z0−− 2 2 3∂u∂x0−2 3∂u∂y0+1 3∂u∂z0− 2 2 3∂u∂x0+1 3∂u∂y0−2 3∂u∂z0−− 1 3∂u∂x0+2 3∂u∂y0+2 3∂u∂z0= g(x0, y0, z0),где g(x0, y0, z0) = f 2 3x0−2 3y0+1 3z0,2 3x0+1 3y0−2 3z0,1 3x0+2 3y0+2 3z0После приведения подобных получим3∂2u∂y02+ 6∂2u∂z02− 3∂u∂x0= g(x0, y0, z0).Полученное уравнение – это уравнение параболического типа в канониче- ской форме.Подводя итоги, отметим, что рассмотренные ранее уравнения принад- лежат разным типам. Для этих уравнений по-разному ставились допол- нительные условия. Правила их постановки целиком определяются типом уравнения.Для уравнений эллиптического типа ставятся краевые условия на всей границе области определения уравнения так, как это делалось для уравне- ния Лапласа.Для уравнений гиперболического типа выделяют одну переменную(обычно она соответствует времени) и ставят по ней 2 начальных усло- вия, а по остальным переменным ставятся краевые условия так, как это делалось для волнового уравнения.Для уравнений параболического типа со временем связывают перемен- ную, по которой отсутствует вторая производная. По ней ставится началь- ное условие, а по остальным – краевые, как для уравнения теплопроводно- сти.4.12. Корректные и некорректные задачи для уравнения ЛапласаВ 4.11 были сформулированы условия, которые следует накладывать на уравнения разных типов. Возникает естественный вопрос: что будет,142 если условия поставлены по-другому, например, для уравнения Лапласа поставлена начально-краевая задача?Введем понятие корректности задачи. Задача называется корректной,если выполнены 3 условия:1) решение задачи существует;2) решение задачи единственно;3) решение задачи устойчиво, т. е. малые изменения данных (коэффи- циентов уравнения, краевых условий и др.) приводят к малым изменениям решения.Вопросы существования и единственности решения рассматривались ранее. Теоремы единственности формулировались и доказывались. Теоре- мы существования в явной форме не формулировались, однако для ряда задач решение было построено и тем самым доказано его существование.Теперь приведем примеры устойчивой и неустойчивой задач.Пример 4.17. Пусть область Ω – прямоугольник [a, b] × [c, d], Γ –граница Ω. Рассмотрим в Ω уравнение Лапласа с краевыми условиями Ди- рихле(∆u = 0,uΓ=ϕ(4.27)Пусть u1и u2– решение задачи (4.27) с краевыми условиямиϕ1иϕ2соответственно. Тогда их разность u = u1− u2является решением задачи(∆u = 0,uΓ=ϕ1−ϕ2Так как решение этой задачи – гармоническая функция, то в силу принципа максимума |u(x, y)| ≤ maxS|ϕ1−ϕ2| для всех (x, y) ∈ Ω. Следовательно,решение устойчиво по начальным данным и задача (4.27) корректна.Пример 4.18. В области Ω = [0,π] × [0, a] рассмотрим следующую задачу:∆u = 0,u(0, y) = u(π, y) = 0,u(x, 0) =1n sin(nx),∂u∂y(x, 0) = 0,(4.28)где n – натуральное число.В этой задаче на границе y = 0 поставлены 2 условия, а на границе y = a условий нет, т. е. поставлена начально-краевая задача так, как это делается для волнового уравнения. Будем решать задачу (4.28) методом143 Фурье. Вспомогательная задача Штурма–Лиувилля(−y00(x) = λy(x),y(0) = y(π) = 0,ее решение λk= k2, y k(x) = sin(kx), k = 1, 2, ... . Решение задачи (4.28)ищем в виде u(x, y) =+∞Xk=1c k(y) sin(kx). После подстановки ряда в задачу(4.28) для коэффициентов c k(y) получаемc00k(y) − k2c k(y) = 0,c k(0) =αk,c0k(0) = 0,(4.29)гдеαk– коэффициент Фурье начального условия1n sin(nx). Заметим, что функция sin(nx) входит в ортогональную систему {sin(kx)}+∞k=1, следова- тельно, ее ряд Фурье совпадает с самой функцией, т. е.αk= 0, если k 6= n иαn=1n. Отсюда c k(y) = 0, если k 6= n, так как задача Коши для одно- родного уравнения с нулевыми начальными условиями имеет только три- виальное решение. Для k = n имеем c00n− n2c n= 0,c n(0) =1n, c0n(0) = 0.Отсюда общее решение cn(y) = Ae ny+ Be−nyИз начальных условий получаемA + B =1n,A − B = 0.Следовательно, A = B =1 2n, c n=1 2n(e ny+ e−ny) и решение задачи (4.28)u(x, y) =1 2n(e ny+ e−ny) sin(nx).Таким образом, решение задачи (4.28) существует и можно доказать, что оно единственно. Однако если n устремить к бесконечности, то lim n→+∞1n sin(nx) = 0 для всех x ∈ [0,π], а ∀y > 0 найдутся точки, где144 решение u(x, y) неограниченно возрастает, так как lim n→+∞e ny2n= +∞, т. е.устойчивости по начальным данным нет и начально-краевая задача (4.28)некорректна.5. ЭЛЕМЕНТЫ ВАРИАЦИОННОГО ИСЧИСЛЕНИЯ5.1. Задача о минимуме функционалаИнтегральным функционалом или просто функционалом назовем отоб- ражение видаJ (y) =bZaF (x, y, y0)dx,(5.1)где F – заданная функция трех переменных. В общем случае функциона- лом называют отображение некоторого класса функций одной или несколь- ких переменных на множество действительных чисел. Областью определе- ния функционала (5.1) назовем множество функций, дифференцируемых на промежутке [a, b] и удовлетворяющих краевым условиям ДирихлеD(J ) =n y ∈ C1[a,b], y(a) = t1, y(b) = t2oПоставим задачу нахождения минимума функционала (5.1), т. е. на- хождения такой функции y0∈ D(J), которая удовлетворяет следующему условию: ∃ε> 0, такое, что ∀y ∈ D(J )y6=y0и max x∈[a,b]|y(x) − y0(x)| <εвыполнено неравенство J (y) > J (y0).Предположим, что на функции y0функционал J (y) достигает мини- мума. Пусть v(x) – произвольная дифференцируемая функция, удовлетво- ряющая условию v(a) = v(b) = 0. Рассмотрим следующую функцию:ϕ(t) = J (y0+ tv) =bZaF (x, y0+ tv, y0 0+ tv0)dx.(5.2)Если функционал J (y) имеет минимум в y0, тоϕ(t) будет иметь минимум при t = 0 и, следовательно,ϕ0(0) = 0. Найдем производнуюϕ(t), не обос- новывая возможность дифференцирования интеграла (5.2) по параметруϕ0(0) =bZa ∂F∂y v +∂F∂y0v0dx,145 и проинтегрируем второе слагаемое по частям:ϕ0(0) =bZa∂F∂y v dx +∂F∂y0v ba−bZa ddx ∂F∂y0v dx.В силу условия v(a) = v(b) = 0 получаемϕ0(0) =bZa ∂F∂y−d dx ∂F∂y0v dx = 0.Так как функция v(x) произвольна, последнее равенство возможно только при∂F∂y−d dx ∂F∂y0= 0.(5.3)Уравнение (5.3) называется уравнением Эйлера для функционала (5.1).Пример 5.1. Найти линию наименьшей длины, соединяющую точкиA(a, t1) и B(b, t2).Как известно, длина кривой y(x):J (y) =bZa p1 + (y0)2dx.Найдем уравнение Эйлера для F (x, y, y0) =p1 + (y0)2. Имеем:∂F∂y= 0,∂F∂y0=y0p1 + (y0)2,d dx ∂F∂y0=y00(1 + (y0)2)3/2Тогда, согласно (5.3), уравнение Эйлера y00= 0.Его общее решение y = c1x + c2– прямая линия. Если подставить сюда координаты точек A и B, получим значения констант c1и c2Пример 5.2. Задача о наименьшей поверхности вращения: опреде- лить кривую с заданными граничными точками, от вращения которой во- круг оси абсцисс образуется поверхность наименьшей площади.Как известно, площадь поверхности вращенияS = 2πbZa yp1 + (y0)2dx.146 Найдем минимум функционалаJ (y) =bZa yp1 + (y0)2dy.Имеем:F (x, y, y0) = y p1 + (y0)2,∂F∂y=p1 + (y0)2,∂F∂y0=yy0p1 + (y0)2;d dx ∂F∂y0=(y0)2p1 + (y0)2+yy00p(1 + (y0)2)3Уравнение Эйлера будет p1 + (y0)2−(y0)2p1 + (y0)2−yy00p(1 + (y0)2)3= 0.Приведем уравнение к общему знаменателю(1 + (y0)2)2− (y0)2(1 + (y0)2) − yy00p(1 + (y0)2)3= 0и получим yy00= (y0)2+ 1.Данное уравнение не содержит в явном виде x, поэтому его можно свести к уравнению первого порядка заменой y0= p(y), y00=dp dy y0=dp dy p. Получаем ydp dy p = p2+ 1.Разделяем переменные и интегрируемZp p2+ 1dp =Zdy y,или1 2ln(p2+ 1) = ln |y| + ln |c1|.Отсюда p2+ 1 = c2 1y2и получаем уравнение y0= ±q c2 1y2− 1.147 Разделяем переменные и интегрируемZdy±pc2 1y2− 1=Zdx,или1c1lnc1y +q c2 1y2− 1 = x + c2Поскольку c1, c2– произвольные постоянные, знаки ± опущены. Разрешим полученное уравнение относительно y:c1y +q c2 1y2− 1 = e c1x+c3,где c3= c1c2. Перенесем c1y в правую часть и возведем в квадрат c2 1y2− 1 = e c1x+c3− c1y2Отсюда −1 = (e c1x+c3)2− 2c1ye c1x+c3и окончательно y =1 2c1e c1x+c3+ e−(c1x+c3)=1c1ch(c1x + c3).Получилось так называемое уравнение цепной линии. Значения констант c1и c3определяются из граничных условий.Пример 5.3. Рассмотрим функционалJ (y) =bZap(x)(y0)2+ q(x)y2− 2f (x)y dx.(5.4)Найдем производные подынтегральной функцииF (x, y, y0) = p(x)(y0)2+ q(x)y2− 2f (x)y :∂F∂y= 2q(x)y − 2f (x),d dx ∂F∂y0= (2p(x)y0)0Отсюда получаем уравнение Эйлера−(p(x)y0)0+ q(x)y = f (x).Получилось хорошо знакомое из разд. 1 уравнение с симметричным диф- ференциальным оператором.Пример 5.4. Рассмотрим дифференциальное уравнениеL(y) ≡ −(p(x)y0)0+ q(x)y = f (x)148 с краевыми условиями y(0) = 0, y(b) = 0. Построим функционалJ (y) = (L(y), y) − 2(f, y) =bZa−(p(x)y0)0y + q(x)y2 dx − 2bZa f (x)y dx.Проинтегрируем первое слагаемое по частям−bZa(p(x)y0)0y dx = −p(x)y0y ba+bZa p(x)y0y0dx.С учетом краевых условий получаемJ (y) =bZap(x)(y0)2+ q(x)y2− 2f (x)y dx,т. е. функционал вида (5.4).Пусть L(y) – линейный симметричный положительно-определенный дифференциальный оператор. Справедлива следующая теорема.Теорема 5.1. Задача решения операторного уравненияL(y) = f (x)и задача поиска минимума функционалаJ (y) = (L(y), y) − 2(f, y)эквивалентны.Доказательство. Возьмем линейную комбинацию функций y +αv и распишем функционал J (y +αv), используя свойства скалярного произве- дения и линейность L:J (y +αv) = (L(y +αv), y +αv) − 2(f, y +αv) == (L(y), y) + (L(αv), y) + (L(y),αv) + (L(αv),αv) − 2(f, y) − 2(f,αv).В силу симметричности оператора L имеем(L(αv), y) = (αv, L(y)) =α(L(y), v).ПоэтомуJ (y +αv) = (L(y), y) − 2(f, y) + 2α(L(y), v) − 2α(f, v) +α2(L(v), v),илиJ (y +αv) = J (y) + 2α(L(y) − f, v) +α2(L(v), v).149 Пусть y(x) – решение уравнения L(y) = f , тогдаJ (y +αv) = J (y) +α2(L(v), v).Оператор L – положительно-определенный, т. е. (L(v), v) > 0 для любой функции v(x) 6= 0 (kvk > 0). Следовательно, J (y +αv) > J (y), т. е. на функции y функционал J достигает минимума.Обратное утверждение: пусть J (y) – минимальное значение функцио- нала. Тогда при любомαдолжно выполняться неравенство2α(L(y) − f, v) +α2(L(v), v) > 0.Но при изменении знакаαлевая часть неравенства будет менять знак,если только (L(y) − f, v) 6= 0. Следовательно, для произвольной функции v скалярное произведение (L(y) − f, v) = 0. Это возможно только, еслиL(y) = f .К доказательству теоремы следует сделать одно принципиальное за- мечание. Дело в том, что области определения дифференциального опера- тора L(y) и функционала J (y) не совпадают. Функционал J (y) определен на функциях, имеющих первую производную и удовлетворяющих краевым условиям. Решение же дифференциального уравнения L(y) = f должно быть дважды дифференцируемым. В случае, когда функция y(x), мини- мизирующая функционал, дважды дифференцируема, она является и ре- шением уравнения L(y) = f в обычном смысле, и такое решение принято называть классическим решением уравнения.Если же у функции y(x) непрерывна только первая производная, y(x)называют обобщенным решением уравнения L(y) = f .5.2. Задача о минимуме функционала в многомерном случаеПусть u(x, y) – дифференцируемая функция двух переменных, удо- влетворяющая на границе области Ω краевым условиям Дирихле. Введем в рассмотрение функционалJ (u) =Z ZΩF (x, y, u, u0x, u0y)dx dy.(5.5)Требуется найти функцию u0(x, y), на которой функционал (5.5) достигает минимального значения. По аналогии с одномерной задачей напишем для функционала (5.5) уравнение Эйлера∂F∂u−d dx ∂F∂u0x−d dy ∂F∂u0y= 0.(5.6)Аккуратный вывод уравнения (5.6) см., например, в [10].150 Решением уравнения Эйлера и будет функция, обеспечивающая мини- мум функционала (5.5).Пример 5.5. Рассмотрим функционалJ (u) =Z ZΩ" ∂u∂x2+ ∂u∂y2+ qu2− 2f u#dx dy(5.7)и найдем для него уравнение Эйлера:F (x, y, u, u0x, u0y) = (u0x)2+ (u0y)2+ qu2− 2f u,∂F∂u= 2qu − 2f,∂F∂u0x= 2u0x,d dx ∂F∂u0x= 2∂2u∂x2,∂F∂u0y= 2u0y,d dy ∂F∂u0y= 2∂2u∂y2Таким образом, уравнение Эйлера будет−∂2u∂x2−∂2u∂y2+ qu = f.(5.8)Как и в одномерном случае, если L(u) – симметричный положительно- определенный оператор, то решение операторного уравнения L(u) = f эк- вивалентно задаче оминимуме функционалаJ (u)== (L(u), u) − 2(f, u).Пример 5.6. Построим функционал для уравнения (5.8). Если функ- ция q(x, y) ≥ q0> 0, то оператор L(u) = −∂2u∂x2−∂2u∂y2+ q(x, y)u будет симметричным и положительно-определенным. К уравнению (5.8) добавим однородные краевые условия Дирихле uΓ= 0, где Γ – граница области Ω.Запишем функционалJ (u) = (Lu, u)−2(f, u) =Z ZΩ−∂2u∂x2−∂2u∂y2+ quu dx dy −2Z ZΩf u dx dy.Первые слагаемые в интеграле преобразуем по формуле Грина−Z ZΩ ∂2u∂x2+∂2u∂y2u dx dy =Z ZΩ(grad u, grad u)dx dy −ZΓu∂u∂n ds,151 где n – внешняя нормаль к границе Γ. С учетом краевых условий имеем−Z ZΩ ∂2u∂x2+∂2u∂y2u dx dy =Z ZΩ" ∂u∂x2+ ∂u∂y2#dx dy и, следовательно,J (u) =Z ZΩ" ∂u∂x2+ ∂u∂y2+ qu2− 2f u#dx dy,т. е. получили тот же функционал (5.7).5.3. Решение задачи о минимуме функционала методом РитцаВ предыдущих параграфах было показано, что задача нахождения ми- нимума функционала эквивалентна решению краевой задачи для диффе- ренциального уравнения Эйлера. Однако существуют методы нахождения минимума, не использующие уравнение Эйлера. Их принято называть пря- мыми. Более того, в ряде случаев от исходной краевой задачи переходят к задаче поиска минимума функционала. Именно такую задачу и рассмот- рим.Пусть L(u) – симметричный положительно-определенный оператор.Поставим краевую задачу ДирихлеL(u) = f ; uΓ= 0,(5.9)где Γ – граница области Ω. Решение этой задачи будем искать как минимум функционалаJ (u) = (L(u), u) − 2(f, u).(5.10)Область определения функционала D(J ) состоит из дифференцируемых вΩ функций, удовлетворяющих краевым условиям uΓ= 0. Возьмем линей- но независимый набор функций {ϕk}n k=1;ϕk∈ D(J), и рассмотрим множе- ство линейных комбинаций {α1ϕ1+α2ϕ2+ ... +αnϕn} = Dn. МножествоDn является подпространством линейного пространства D(J ), а функцииϕ1, ...,ϕn будут базисом этого подпространства. Любой элемент u ∈ Dn од- нозначно определяется набором коэффициентовα1, ...,αn– координатами u в базисеϕ1, ...,ϕnНайдем минимум функционала J (u) на множестве Dn, т. е. найдем та-152 кой набор коэффициентовα1, ...,αn, который обеспечит минимум функцииΦ(α1,α2, ...,αn) = J (α1ϕ1+α2ϕ2+ ... +αnϕn) == (L(α1ϕ1+ ... +αnϕn),α1ϕ1+ ... +αnϕn) − 2(f,α1ϕ1+ ... +αnϕn) ==α2 1(L(ϕ1),ϕ1) +α1α2(L(ϕ1),ϕ2) + ... +α1αn(L(ϕ1),ϕn) ++α2α1(L(ϕ2),ϕ1) +α2 2(L(ϕ2),ϕ2) + ... +α2αn(L(ϕ2),ϕn) +− 2α1(f,ϕ1) − 2α2(f,ϕ2) − ... − 2αn(f,ϕn).Найдем все частные производные∂Φ∂αk, k = 1, 2, ..., n, и приравня- ем их нулю. При этом учитываем, что оператор L симметричный, т. е.(L(ϕk),ϕm) = (ϕk, L(ϕm)). В итоге получится система линейных урав- ненийα1(L(ϕ1),ϕ1) +α2(L(ϕ2),ϕ1) + ... +αn(L(ϕn),ϕ1) = (f,ϕ1),α1(L(ϕ1),ϕ2) +α2(L(ϕ2),ϕ2) + ... +αn(L(ϕn),ϕ2) = (f,ϕ2),α1(L(ϕ1),ϕn) +α2(L(ϕ2),ϕn) + ... +αn(L(ϕn),ϕn) = (f,ϕn).(5.11)Покажем, что система (5.11) однозначно разрешима. Для этого рас- смотрим соответствующую однородную систему и предположим, что она имеет нетривиальное решениеα∗1,α∗2, ...,α∗n, т. е. выполняются равенстваα∗1(L(ϕ1),ϕ1) +α∗2(L(ϕ2),ϕ1) + ... +α∗n(L(ϕn),ϕ1) = 0,α∗1(L(ϕ1),ϕ2) +α∗2(L(ϕ2),ϕ2) + ... +α∗n(L(ϕn),ϕ2) = 0,α∗1(L(ϕ1),ϕn) +α∗2(L(ϕ2),ϕn) + ... +α∗n(L(ϕn),ϕn) = 0.В силу линейности скалярного произведения и оператора L имеем(L(α∗1ϕ1+α∗2ϕ2+ ... +α∗nϕn),ϕ1) = 0,(L(α∗1ϕ1+α∗2ϕ2+ ... +α∗nϕn),ϕ2) = 0,(L(α∗1ϕ1+α∗2ϕ2+ ... +α∗nϕn),ϕn) = 0.Умножим первое равенство наα∗1, второе – наα∗2и т. д., затем сложим и получим(L(α∗1ϕ1+α∗2ϕ2+ ... +α∗nϕn),α∗1ϕ1+α∗2ϕ2+ ... +α∗nϕn) = 0.Так как оператор L – положительно-определенный (т. е. ∀u, u 6= 0,(L(u), u) > 0), тоα∗1ϕ1+α∗2ϕ2+ ... +α∗nϕn= 0. Полученное равенство противоречит условию линейной независимостиϕ1,ϕ2, ...,ϕn. Следова- тельно, однородная система имеет только тривиальное решение, а значит,153 система (5.11) разрешима единственным образом. Решив систему (5.11),получим элемент u n=α1ϕ1+α2ϕ2+ ... +αnϕn, который можно рассмат- ривать как приближенное решение задачи о минимуме функционала (5.10),а также как приближение обобщенного решения краевой задачи (5.9). Ес- ли расширить систему базисных функцийϕ1,ϕ2, ...,ϕn,ϕn+1, ...,ϕn1и найти элемент u n1, обеспечивающий минимум функционала на множествеDn1= {α1ϕ1+ ... +αn1ϕn1}, то J(u n1) ≤ J (u n). Минимум функционала не может увеличиться, потому что u n∈ Dn1. Затем можно еще добавить ба- зисных функций, получить u n2и т. д. Получится последовательность при- ближенных решений u n, u n1, u n2, ... .Доказательство сходимости последовательности к точному решению задачи выходит за рамки настоящего пособия, его можно найти, например,в [11].Описанный метод построения последовательности приближенных ре- шений называется методом Ритца.Отметим, что в методе Ритца от базисных функцийϕk требуется на- личие только первых производных.Пусть оператор L имеет видL(y) = −(p(x)y0)0+ q(x)y,функция y(x) задана на промежутке [a, b] и удовлетворяет краевым усло- виям y(a) = y(b) = 0. Тогда скалярное произведение(L(ϕk),ϕm) =bZa[−(p(x)ϕ0k)0+ q(x)ϕk]ϕm dx == −p(x)ϕ0kϕm ba+bZa(p(x)ϕ0kϕ0m+ q(x)ϕkϕm)dx.С учетом краевых условий(L(ϕk),ϕm) =bZa(p(x)ϕ0kϕ0m+ q(x)ϕkϕm)dx.В двумерном случае, например для оператораL(u) = −∆u + q(x, y)u,скалярные произведения будут выглядеть следующим образом:(L(ϕk),ϕm) =Z ZΩ ∂ϕk∂x∂ϕm∂x+∂ϕk∂y∂ϕm∂y+ q(x, y)ϕkϕmdx dy.154 5.4. Понятие о методе ГалеркинаПусть L(u) – дифференциальный оператор 2-го порядка (не обязатель- но симметричный). Поставим краевую задачуL(u) = f,(5.12)где u удовлетворяет однородным краевым условиям Дирихле. Возьмем произвольную дифференцируемую функцию v, удовлетворяющую одно- родным краевым условиям Дирихле, и скалярно умножим уравнение (5.12)на v:(L(u), v) = (f, v),(5.13)или, что то же самое, (L(u) − f, v) = 0. Поскольку скалярное произведе- ние – это интеграл, соотношение (5.13) принято называть интегральным тождеством.Функцию u, удовлетворяющую интегральному тождеству с произволь- ной функцией v, будем называть обобщенным решением краевой задачи(5.12).Очевидно, что если краевая задача (5.12) имеет обычное решение, то это решение будет и обобщенным. Однако обобщенное решение может быть дифференцируемо только один раз, тогда как решение задачи (5.12) долж- но быть дифференцируемо 2 раза. Такой случай уже встречался при рас- смотрении задачи о минимуме функционала. Более того, в 5.1 доказатель- ство эквивалентности задачи о минимуме функционала и краевой задачи основывалось на выполнении соотношения (L(u)−f, v) = 0, т. е. интеграль- ного тождества (5.13). Таким образом, для симметричного оператора зада- ча нахождения минимума функционала и задача поиска функции, удовле- творяющей интегральному тождеству, эквивалентны. Однако соотношение(5.13) можно написать и для несимметричного оператора и не связывать его с каким-либо функционалом.Метод Галеркина состоит в построении решения задачи на конечно- мерном подпространстве так же, как и в методе Ритца. Возьмем систему из n линейно независимых дифференцируемых функций, удовлетворяю- щих однородным краевым условиямϕ1,ϕ2, ...,ϕn. Рассмотрим множество всех линейных комбинаций этих функций Dn= {α1ϕ1+α2ϕ2+ ... +αnϕn}.Будем искать функцию u n=α1ϕ1+α2ϕ2+ ... +αnϕn, удовлетворяющую соотношению (5.13) для произвольной функции v ∈ Dn:(L(α1ϕ1+α2ϕ2+ ... +αnϕn), v) = (f, v).(5.14)Если интегральное тождество будет выполнено для всех базисных функций пространства Dn, то в силу линейности скалярного произведения оно будет выполнено и для любой функции v ∈ Dn. Последовательно подставляя в155 (5.14) вместо v базисные функцииϕ1,ϕ2, ...,ϕn, получим систему линей- ных уравненийα1(L(ϕ1),ϕ1) +α2(L(ϕ2),ϕ1) + ... +αn(L(ϕn),ϕ1) = (f,ϕ1),α1(L(ϕ1),ϕ2) +α2(L(ϕ2),ϕ2) + ... +αn(L(ϕn),ϕ2) = (f,ϕ2),α1(L(ϕ1),ϕn) +α2(L(ϕ2),ϕn) + ... +αn(L(ϕn),ϕn) = (f,ϕn).(5.15)В случае симметричного оператора L система (5.15) совпадает с системой,полученной методом Ритца. Решая систему (5.15), получим решение задачи un=α1ϕ1+α2ϕ2+ ... +αnϕn∈ Dn, которое можно принять за приближен- ное решение задачи (5.13) и за приближенное обобщенное решение задачи(5.12). Затем набор базисных функций расширяется, строится новое при- ближенное решение u n1и т. д. В итоге, как и в методе Ритца, получаем последовательность приближенных решений u n, u n1, ... . Вопросы сходимо- сти приближенного решения к точному, как и разрешимость системы (5.15)в общем случае, обсуждаются в известных монографиях, например в [11].Покажем, как строится интегральное тождество для двумерной крае- вой задачи. Пусть в области Ω задано уравнение эллиптического типа−∂2u∂x2−∂2u∂y2+ b1(x, y)∂u∂x+ b2(x, y)∂u∂y+ q(x, y)u = f (x, y).(5.16)На границе области Γ заданы однородные условия Дирихле uΓ= 0.Заметим, что оператор задачи (5.16) не симметричен и поэтому метод Рит- ца к ней не применим. Умножим уравнение (5.16) на произвольную, удовле- творяющую однородным краевым условиям функцию v и проинтегрируем по Ω:−Z ZΩ ∂2u∂x2+∂2v∂y2v dx dy +Z ZΩb1(x, y)∂u∂x+ b2(x, y)∂u∂yv dx dy ++Z ZΩq(x, y)uv dx dy =Z ZΩf (x, y)v dx dy.(5.17)Первый интеграл преобразуем по формуле интегрирования по частям:−Z ZΩ ∂2u∂x2+∂2v∂y2v dx dy =156 =Z ZΩ ∂u∂x∂v∂y+∂u∂y∂v∂ydx dy −ZΓ∂u∂n v ds,где n – нормаль к границе Γ. Так как функция v на границе обращается в нуль, то и последний интеграл тоже равен нулю, и формула (5.17) примет видZ ZΩ ∂u∂x∂v∂x+∂u∂y∂v∂ydx dy +Z ZΩb1(x, y)∂u∂x+ b2(x, y)∂u∂yv dx dy ++Z ZΩq(x, y)uv dx dy =Z ZΩf (x, y)v dx dy.Полученное соотношение и является интегральным тождеством для урав- нения (5.16).5.5. Пример решения уравнения Пуассона методом конечных элементовВ изложении метода конечных элементов будем следовать монографии[12]. Рассмотрим задачу Дирихле для уравнения Пуассона в прямоуголь- нике Ω = [0, a] × [0, b]−∂2u∂x2−∂2u∂y2= f (x, y),uΓ= 0.(5.18)Пусть v – произвольная функция, удовлетворяющая однородным краевым условиям. Для задачи (5.18) запишем интегральное тождество метода Га- леркинаZ ZΩ ∂u∂x∂v∂x+∂u∂y∂v∂ydx dy =Z ZΩf v dx dy.(5.19)Далее в методе Галеркина необходимо выбрать систему базисных функ- ций. Для этого зададим в области Ω сетку. Пусть заданы шаги сетки h x=a nи h y=b m. Обозначим x k= kh x, y k= ih y, k = 0, 1, ..., n, i = 0, 1, ..., m.В каждом прямоугольнике построенной сетки проведем диагональ из левого нижнего угла. Таким образом, область Ω оказалась разбита на тре- угольники (рис. 5.1). Такой вид сетки принято называть триангуляцией.Введем в рассмотрение функцииϕki(x, y), обладающие следующими свойствами: функцияϕki равна единице в узле (x k, y i) и нулю во всех157 остальных узлах сетки; функция линейна во всех треугольниках сетки(рис. 5.2). Такие функции называют функциями-палатками.1   2   3   4   5   6   7   8   9   10

a
x
x
x
y
x
x
0 1
2 3
0 1
2 3
y
y
y
0
b
y
Рис. 5.1
Рис. 5.2
h
x
h
x
y
y
y
h
h
x
k
y
y
x
k
i
x
i
i
V
IV
VI
I
II
III
k
+1
−1
−1
+1
Рис. 5.3
Понятно, что функция
ϕ
ki будет отлична от нуля только в тех треуголь- никах сетки, которые имеют вершину в точке (x k
, y i
). Занумеруем треугольни- ки вокруг узла (x k
, y i
) как показано на рис. 5.3.
В таблице в явном виде приведе- ны функция
ϕ
ki
(x, y) и ее производные для каждого треугольника сетки.
Система функций {
ϕ
ki
}
k=1,...,n−1
i=1,...,m−1
линейно-независима, и ее можно взять в качестве базисных функций ме- тода Галеркина. Решение интегрального тождества (5.19) будем искать в виде
˜
u(x, y) =
n−1
X
k=1
m−1
X
i=1
u ki
ϕ
ki
(x, y).
Так как в каждом узле сетки отлична от нуля только одна базисная функ- ция, то ˜
u(x k
, y i
) = u ki
. Т. е., если найти коэффициенты решения задачи u ki
,
то фактически получим приближенные значения решения в узлах сетки.
Для того чтобы найти u ki
, запишем систему линейных уравнений, анало- гичную (5.15), учитывая при этом, что оператор задачи (5.18) есть L = −∆:
158

Номер треугольника
ϕ
ki
(x, y)
∂ϕ
ki
∂x
∂ϕ
ki
∂y
I
1 −
y − y i
h y
0
−
1
h y
II
1 −
x − x k
h x
−
1
h x
0
III
1 −
x − x k
h x
+
y − y i
h y
−
1
h x
1
h y
IV
1 +
y − y i
h y
0 1
h y
V
1 +
x − x k
h x
1
h x
0
VI
1 +
x − x k
h x
−
y − y i
h y
1
h x
−
1
h y













































u
11
(−∆
ϕ
11
,
ϕ
11
) + u
21
(−∆
ϕ
21
,
ϕ
11
) + ...
... + u n−1 m−1
(−∆
ϕ
n−1 m−1
,
ϕ
11
) = (f,
ϕ
11
),
u
11
(−∆
ϕ
11
,
ϕ
21
) + u
21
(−∆
ϕ
21
,
ϕ
21
) + ...
... + u n−1 m−1
(−∆
ϕ
n−1 m−1
,
ϕ
21
) = (f,
ϕ
21
),
u
11
(−∆
ϕ
11
,
ϕ
n−1 1
) + u
21
(−∆
ϕ
21
,
ϕ
n−1 1
) + ...
... + u n−1 m−1
(−∆
ϕ
n−1 m−1
,
ϕ
n−1 1
) = (f,
ϕ
n−1 1
),
u
11
(−∆
ϕ
11
,
ϕ
12
) + u
21
(−∆
ϕ
21
,
ϕ
12
) + ...
... + u n−1 m−1
(−∆
ϕ
n−1 m−1
,
ϕ
12
) = (f,
ϕ
12
),
u
11
(−∆
ϕ
11
,
ϕ
n−1 m−1
) + u
21
(−∆
ϕ
21
,
ϕ
n−1 m−1
) + ...
... + u n−1 m−1
(−∆
ϕ
n−1 m−1
,
ϕ
n−1 m−1
) = (f,
ϕ
n−1 m−1
).
(5.20)
Возьмем в системе (5.20) строчку с номером (k, i) и вычислим все вхо- дящие в нее скалярные произведения, которые с учетом (5.19) имеют вид
(−∆
ϕ
lj
,
ϕ
ki
) =
Z Z
Ω
∂
ϕ
lj
∂x
∂
ϕ
ki
∂x
+
∂
ϕ
lj
∂y
∂
ϕ
ki
∂y

dx dy.
(5.21)
Так как функции
ϕ
ki и
ϕ
lj отличны от нуля только в шести треуголь- никах сетки, то интеграл (5.21) равен нулю, если шестерки треугольников,
примыкающих к узлам (ki) и (lj), не пересекаются. Выпишем все случаи,
159
когда пересечение примыкающих треугольников есть. В последующих фор- мулах нумерация треугольников выбрана согласно рис. 5.2. Очевидно, что площади всех треугольников равны
1 2
h x
h y
:
(−∆
ϕ
ki
,
ϕ
ki
) =
Z Z
Ω
∂
ϕ
ki
∂x
∂
ϕ
ki
∂x
+
∂
ϕ
ki
∂y
∂
ϕ
ki
∂y

dx dy =
=
Z Z
∆
I

−
1
h y

2
dx dy +
+
Z Z
∆
II

−
1
h x

2
dx dy +
Z Z
∆
III
"

−
1
h x

2
+
1
h y

2
#
dx dy +
+
Z Z
∆
IV
1
h y

2
dx dy +
Z Z
∆
V
1
h x

2
dx dy +
Z Z
∆
VI
"
1
h x

2
+

−
1
h y

2
#
dx dy =
= 2
h x
h y
+
h y
h x

;
(−∆
ϕ
k i+1
,
ϕ
ki
) =
Z Z
Ω
∂
ϕ
k i+1
∂x
∂
ϕ
ki
∂x
+
∂
ϕ
k i+1
∂y
∂
ϕ
ki
∂y

dx dy =
=
Z Z
∆
I

−
1
h x
0 +
1
h y
−1
h y

dx dy +
Z Z
∆
VI

0 1
h x
+
1
h y
−1
h y

dx dy = −
h x
h y
;
(−∆
ϕ
k+1 i+1
,
ϕ
ki
) =
Z Z
Ω
∂
ϕ
k+1 i+1
∂x
∂
ϕ
ki
∂x
+
∂
ϕ
k+1 i+1
∂y
∂
ϕ
ki
∂y

dx dy =
=
Z Z
∆
I
1
h x
0 + 0
−1
h y

dx dy +
Z Z
∆
II

0
−1
h x
+
1
h y
0

dx dy = 0.
Аналогично вычисляем
(−∆
ϕ
k+1 i
,
ϕ
ki
) = −
h y
h x
;
(−∆
ϕ
k i−1
,
ϕ
ki
) = −
h x
h y
;
(−∆
ϕ
k−1 i−1
,
ϕ
ki
) = 0;
(−∆
ϕ
k−1 i
,
ϕ
ki
) = −
h y
h x
Интеграл в правой части уравнения обозначим f
ki
= (f,
ϕ
ki
) =
Z Z
P ∆
f
ϕ
ki dx dy.
160

Таким образом, уравнение с номером (ki) примет вид
−
h x
h y
u k i−1
−
h y
h x
u k−1 i
+ 2
h x
h y
+
h y
h x

u ki
−
h y
h x
u k+1 i
−
h x
h y
u k i+1
= f ki
Тогда систему линейных уравнений (5.20) можно переписать в виде



















2
h x
h y
+
h y
h x

u
11
−
h y
h x
u
21
− ... −
h x
h y
u
12
− ... = f
11
,
−
h y
h x
u
11
+ 2
h x
h y
+
h y
h x

u
21
−
h y
h x
u
31
− ... −
h x
h y
u
22
... = f
21
,
−
h x
h y
u n−1 m−2
−...−
h y
h x
u n−2 m−1
+2
h x
h y
+
h y
h x

u n−1 m−1
= f n−1 m−1
(5.22)
Матрица системы (5.22) имеет такую же блочно-трехдиагональную струк- туру, как в системе метода конечных разностей. Более того, если все урав- нения системы (5.22) умножить на h x
h y
, то матрицы систем просто сов- падут. Только в отличие от метода конечных разностей в правых частях уравнения стоит не значение функции f (x, y) в узле сетки, а некоторое усреднение по шести треугольникам, примыкающим к узлу. Основные пре- имущества метода конечных элементов по сравнению с методом конечных разностей проявляются при решении задач в областях сложной формы.
Если необходимо использовать нерегулярную сетку (значения h x
и h y
ме- няются от шага к шагу), то формулы аппроксимации производных значи- тельно усложняются, если же необходима непрямоугольная сетка, то по- строение разностных отношений просто невозможно. В методе конечных элементов коэффициенты линейных уравнений есть интегралы по некото- рым треугольникам, и при любой форме треугольника их вычисление не вызовет принципиальных трудностей.

Список литературы
1. Меркулов А. Л., Трегуб В. Л., Червинская Н. М. Задачи и упраж- нения по математической физике: учеб. пособие. СПб.: Изд-во СПбГЭТУ
„ЛЭТИ“, 2014.
2. Очан Ю. С. Методы математической физики. М.: Высш. шк., 1965.
3. Бодунов Н. А., Пилюгин С. Ю. Дифференциальные уравнения:
учеб. пособие. СПб.: Изд-во СПбГЭТУ „ЛЭТИ“, 2011.
4. Боревич Е. З., Фролова Е. В., Челкак С. И. Ряды Фурье: учеб. посо- бие. СПб.: Изд-во СПбГЭТУ „ЛЭТИ“, 2010.
5. Справочник по специальным функциям / под ред. М. Абрамовица и И. Стиган. М.: Наука, 1979.
6. Ватсон Г. Н. Теория бесселевых функций: в 2 ч. Ч. 1. М.: Изд-во иностр. лит., 1949.
7. Арсенин В. Я. Методы математической физики и специальные функ- ции. М.: Наука, 1984.
8. Тихонов А. Н., Самарский А. А. Уравнения математической физики.
М.: Наука, 2004.
9. Араманович И. Г., Левин В. И. Уравнения математической физики.
М.: Наука, 1969.
10. Эльсгольц Л. Э. Дифференциальные уравнения и вариационное ис- числение. М.: Эдиториал УРСС, 2002.
11. Канторович Л. В., Крылов В. И. Приближенные методы высшего анализа. М.: Физматлит, 1962.
12. Оганесян Л. А., Руховец Л. А. Вариационно-разностные методы ре- шения эллиптических уравнений. Ереван: Изд-во АН АрмССР, 1979.

Содержание
Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1. КРАЕВЫЕ ЗАДАЧИ ДЛЯ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ВТОРОГО ПОРЯДКА . . . . . . . 4 1.1. Постановка краевой задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2. Оператор Штурма–Лиувилля . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3. Задача Штурма–Лиувилля для оператора −y
00
. . . . . . . . . . . . . . . . 13 1.4. Решение краевой задачи для обыкновенного дифференциального уравнения методом Фурье . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.5. Уравнение Бесселя. Функция Бесселя . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.6. Задача на собственные значения для оператора Бесселя . . . . . . . 25 1.7. Оператор Лежандра. Многочлены Лежандра . . . . . . . . . . . . . . . . . . 29 1.8. Присоединенные функции Лежандра . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.9. Решение краевой задачи методом конечных разностей . . . . . . . . . 36 2. УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ . . . . . . . . . . . . . . . . . . . . . . . . 43 2.1. Вывод одномерного уравнения теплопроводности . . . . . . . . . . . . . . 43 2.2. Постановка начально-краевых задач . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.3. Решение краевой задачи для уравнения теплопроводности методом Фурье . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.4. Принцип максимума . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.5. Решение уравнения теплопроводности для бесконечного однородного стержня . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.6. Метод сеток . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3. УРАВНЕНИЕ КОЛЕБАНИЙ СТРУНЫ. ВОЛНОВОЕ
УРАВНЕНИЕ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .79 3.1. Уравнение колебаний струны . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 3.2. Постановка начальных и краевых условий . . . . . . . . . . . . . . . . . . . . . 82 3.3. Колебания бесконечной струны. Метод Даламбера . . . . . . . . . . . . 84 3.4. Колебания полубесконечной струны . . . . . . . . . . . . . . . . . . . . . . . . . . . .92 3.5. Метод Фурье . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4. УРАВНЕНИЯ ЛАПЛАСА И ПУАССОНА . . . . . . . . . . . . . . . . . . . . . 102 4.1 Определения. Постановка краевой задачи. . . . . . . . . . . . . . . . . . . . . 102 4.2. Гармонические функции . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .103 4.3. Формулы Грина . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.4. Свойства гармонических функций . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.5. Теоремы единственности для уравнений Лапласа и Пуассона 110 4.6. Метод Фурье для уравнений Лапласа и Пуассона . . . . . . . . . . . . 112 4.7. Задача на собственные значения для оператора Лапласа . . . . . 116 4.8. Решение задач разложением по собственным функциям оператора Лапласа . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.9. Метод функции Грина . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 163

4.10. Метод сеток для уравнения Пуассона . . . . . . . . . . . . . . . . . . . . . . . . 134 4.11. Классификация уравнений второго порядка . . . . . . . . . . . . . . . . . 138 4.12. Корректные и некорректные задачи для уравнения
Лапласа . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5. ЭЛЕМЕНТЫ ВАРИАЦИОННОГО ИСЧИСЛЕНИЯ . . . . . . . . . . . 145 5.1. Задача о минимуме функционала . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 5.2. Задача о минимуме функционала в многомерном случае . . . . . 150 5.3. Решение задачи о минимуме функционала методом Ритца . . . 152 5.4. Понятие о методе Галеркина . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 5.5. Пример решения уравнения Пуассона методом конечных элементов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
Меркулов Александр Львович,
Трегуб Вера Леонидовна,
Червинская Нина Михайловна
МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
Учебное пособие
Редактор Э. К. Долгатов
Подписано в печать 12.07.16. Формат 60×84 1/16. Бумага офсетная.
Печать цифровая. Гарнитура “Times New Romen”. Печ. л. 10,25.
Тираж 143 экз. Заказ
Издательство СПбГЭТУ “ЛЭТИ”
197376, С.-Петербург, ул. Проф. Попова, 5

Смотрите также файлы