Пушкарь Е.А. - Дифференциальные уравнения, страница 8

Начальное значение можно задать точно и считать, что z0 = 0. Второе слагаемое обусловленоотброшенными членами в формуле Тейлора (4.20) при выводесхемы ломаных (4.22). Оценим это слагаемое сверху; заменяявсе функции под интегралами их модулями и вынося max h(x)за знак интеграла, получим|zm | M(xn ) max hn = O(max hn ),0nmгдеM(xm ) =12xmx0(4.26)xmdτ |y (τ )| exp  |fy |dµ τM4 M3 |xm−x0 |e· |xm − x0 |.2M3(4.27)так как |y (τ )| M4 и экспонента не превышает eM3 (xm −τ ).Следовательно, при h → 0 приближенное решение сходится к точному равномерно (на ограниченном интервале|x − x0 | a) с первым порядком точности.4.5 Метод Рунге-Кутта59Замечание 1.

Оценка погрешности (4.26) — мажорантная.Для функций со знакопеременными производными она можетбыть сильно завышена по сравнению с асимптотической оценкой (4.25).Замечание 2. Первый (экспоненциальный) член в оценке(4.25) характеризует расхождение интегральных кривых (см.рис. 4.1); если он очень велик, то исходная задача Коши плохообусловлена.4.5. Метод Рунге-КуттаЭтот метод позволяет строить схемы различного порядкаточности. Построим семейство схем второго порядка и на егопримере разберем основные особенности метода.В качестве исходного выражения возьмем ряд Тейлора (4.20),оставив в нем член второго порядка точности, порядок которого равен предполагаемому порядку точности схемы1yn+1 = yn + hnyn + h2n yn .2(4.28)Чтобы избежать дифференцирования, запишем вторую производную в виде:y =df (x̃, ỹ) − f (x, y)f (x, y) =,dx∆x(4.29)где нужно подобрать x̃, ỹ и ∆x так, чтобы обеспечить максимально возможную точность.

Возьмем, например, x̃ = xn + γh,ỹ = un + δh.После такой замены, объединив второй и третий члены в(4.28) и заменяя yn сеточной функцией un , приведем это выражение к видуun+1 = un + h [βf (xn , un ) + αf (xn + γh, un + δh)] ,где α, β, γ, δ — параметры, которые нужно определить.(4.30)604Общая теорияРассматривая правую часть (4.30) как функцию h, разложим ее в ряд по степеням h. Получим:f (xn + γh, un + δh) = f (xn , un) + (γfx + δfy )n h + .

. .un+1 = un + h(α + β) f (xn , un) +αh2 (γfx + δfy )n + . . . y1/2[fx +f fy ]Выберем α, β, γ, δ так, чтобы это разложение было наиболееблизко к ряду (4.28). Чтобы правильно передать два первыхчлена формулы Тейлора, необходимо положить11α + β = 1; αγ = ; αδ = f (xn , un ).22Для определения четырех параметров мы получили только три уравнения, так что один параметр остается свободным.Выразим через α все остальные параметрыβ = 1 − α,1,γ=2αf δ = n,2αи подставим их в (4.30). Получим однопараметрическое семейство двучленных схем Рунге-Куттаhhun+1 = un + h (1 − α)f (xn , un ) + αf (xn +, un +fn ) ,2α2α(4.31)где 0 < α 1.

Выбрать параметр α так, чтобы схема (4.31)правильно передавала третий член формулы Тейлора (4.20),невозможно.Погрешность можно исследовать аналогично методу ломаных. При этом доказывается следующий результат.Теорема 4.5. Если f (x, y) непрерывна и ограничена вместесо своими вторыми производными, то решение, полученное посхеме Рунге-Кутта (4.31), равномерно сходится к точному4.5 Метод Рунге-Кутта61решению с погрешностью O(max h2n ), т.е. двучленная схемаnРунге-Кутта имеет второй порядок точности.Формула (4.31) имеет неплохую точность и нередко используется в численных расчетах.

При вычислениях обычно полагают либо α = 1, либо α = 12 .Рассмотрим случай α = 1 (модифицированный метод Эйлера):hhun+1 = un + hf (xn + , un + fn ).22Ее смысл поясняет рис. 4.3.uun + 1u n + 1/2unxnx n + 1/2xn + 1xРис. 4.3. Модифицированный метод ЭйлераВычисления соответствуют следующим геометрическим движениям:(1) Совершаем полушаг h/2 по методу Эйлера, находим коhhординаты “полуцелой” точки xn + , un + fn. В най22денной точке определяем наклон интегральной кривойf (xn+1/2 , un+1/2 ) = un+1/2 ;(2) По найденному значению un+1/2 вычисляем un+1 , используя модифицированное значение наклона u :un+1 = un + hun+1/2 .624Общая теорияКаждое из этих движений происходит по методу Эйлера.Геометрическая интерпретация второго случая (α = 1/2):un+1 = un +h[f (xn , un) + f (xn + h, un + hfn)]2следующая (рис.

4.4):(1) Производится шаг по схеме Эйлера и грубо приближенно находится значение функции ūn+1 = un + hfn и наклон интегральной кривой ūn+1 = f (xn+1 , ūn+1 ) в новойточке.(2) Определяется средний наклон на шаге как полусумманачального и предсказанного конечного значений наклона un+1/2 = 12 (un + ūn+1 ) и по нему уточняется un+1 .Схемы подобного типа нередко называют “предиктор-корректор”.ul2un + 1l1unРис.l 1 || l 2xn4.4. Геометрическаяпредиктор-корректорxn + 1интерпретацияxсхемыМетодом Рунге-Кутта можно построить схемы различного порядка точности. Например, схема ломаных есть схема4.5 Метод Рунге-Кутта63Рунге-Кутта первого порядка точности.

Наиболее употребительны схемы четвертого порядка, образующие семейство четырехчленных схем. Приведем без вывода ту из них, котораязаписана в большинстве стандартных программ ЭВМ:hun+1 = un + (k1 + 2k2 + 2k3 + k4 );6hhk1 = f (xn , un); k2 = f (xn + , un + k1 );22hhk3 = f (xn + , un + k2 ); k4 = f (xn + h, un + hk3 ).22Схемы Рунге-Кутта имеют ряд важных достоинств:(1) Все они (кроме схемы ломаных) имеют хорошую точность.(2) Они являются явными, т.е.

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

Сначала проводят расчет с шагом h, затемс шагом h/2 и погрешность решения с шагом h/2 оцениваютпо формуле: (h)(h/2) (h/2) ≈ u − u·[y]n − u1,2p − 1где p есть порядок аппроксимации. Эта формула может бытьполучена из принципа Рунге.645Уравнения, не разрешенные относительно производной5. Уравнения, не разрешенныеотносительно производной5.1. Основная теорема о решении уравнения,не разрешенного относительно производнойРассмотрим вначале следующий простой пример.Пример.

Пусть уравнение имеет вид: 2dy− 1 = 0.(5.1)dxОчевидно, оно равносильно совокупности двух независимыхуравнений dydx = +1,(5.2)dy=−1.dxПервое из этих уравнений задает поле направлений, наклоненное к оси OX под углом 45◦ , а второе — поле направлений,наклоненное к оси OX под углом 135◦. Уравнению же (5.1) соответствует поле направлений, полученное наложением полейуравнений (5.2). Через каждую точку плоскости (x, y) проходит одна и только одна интегральная линия первого из уравнений (5.2) — прямая, наклоненная к оси OX под углом 45◦ ,yxРис.

5.1. Интегральные линии уравнения (5.1)5.1 Основная теорема об уравнении не разрешенном относительно производной 65одна и только одна интегральная линия второго из уравнений(5.2) — прямая, наклоненная к оси OX под углом 135◦ . Значит, через каждую точку плоскости (x, y) проходят ровно двеинтегральных линии уравнения (5.1) (рис.

5.1).Имеет место следующая общая теорема.Теорема 5.1. Пусть дано уравнениеF (x, y, y ) = 0,(5.3)где функция F (x, y, y ) обладает тремя свойствами:(1) F (x, y, y ) определена в замкнутой и ограниченной области Ḡ пространства (x, y, y ), где она непрерывна.(2) Для некоторой точки (x0 , y0 ) в плоскости (x, y) число различных решений уравнения (5.3) относительноy конечно и равно m.

Пусть этими решениями являются числаb1 , b2 , ..., bm(m > 0).(3) Каждая из точек (x0 , y0 , bi ), (i = 1, ..., m) лежит внутри области G и в некоторой окрестности Ri каждойиз этих точек функция F (x, y, y ) имеет непрерывнуюпроизводную по y и непрерывную производную по y , которая в Ri превосходит по абсолютной величине некоторое положительное число C > 0.Тогда существует окрестность U точки (x0 , y0 ), расположенная в плоскости (x, y), причем через любую точку U проходит m и только m решений уравнения (5.3).Доказательство теоремы.Рассмотрим случай m = 1.Из сделанных предположений следует, что для уравненияF (x, y, y ) = 0 существует точка (x0 , y0 , y0 ) такая, чтоF (x0 , y0 , y0 ) = 0,Fy (x0 , y0 , y0 ) = 0,665Уравнения, не разрешенные относительно производнойтогда по теореме о неявной функции существует окрестностьэтой точки U (x0 , y0 , y0 ), в которой уравнение (5.3) можно разрешить относительно y : y = f1 (x, y), где функция f1 — непрерывная и дифференцируемая по y, следовательно эта функцияудовлетворяет условию Липшица.Действительно, производная∂F∂f1∂y ∂y== − ∂F∂y∂y∂y непрерывна и ограничена в ограниченной замкнутой области.Следовательно, условие Липшица для функции f1 (x, y) выполнено, значит выполнено условие теоремы существованияи единственности и существует единственная функция y(x):y(x0 ) = y0 , y (x0 ) = f1 (x0 , y0 ), такая что F (x0 , y0 , f1 (x0 , y0 )) = 0.5.2.

Решение дифференциальных уравненийв параметрической формеРассмотрим в качестве примера решения дифференциальных уравнений, не разрешенных относительно производной,уравнение видаy = f (y ).(5.4)Если f (y0 ) = 0, то уравнение можно разрешить относительно y . Однако это не всегда возможно сделать аналитически, иимеется ряд существенных трудностей при приведении уравнения (5.4) к виду y = f −1 (y).Тем не менее, покажем, что уравнение (5.4) всегда можнопроинтегрировать (точнее, свести к квадратурам). Для этогоищем решение в виде:y = f (p),(5.5)x = ϕ(p),5.2 Решение дифференциальных уравнений в параметрической форме67то есть в параметрической форме. Так как f (p) известно, то дляdy= p.нахождения решения нужно найти функцию ϕ(p), гдеdxdyфункции y(x),С другой стороны, найдем производнуюdxзаданной параметрически соотношениями (5.5)fp (p)dpdy≡ ,dxϕp (p)dpи приравняем ее параметру p:fp (p)= p.ϕp (p)(5.6)Разрешив (5.6) относительно ϕ (p), находимf (p),ϕ (p) =pзначитϕ(p) =f (p)dp,pи окончательно получаем y = f(p),f (p)dp + Cx=pрешение уравнения y = f (y ) в параметрической форме в видеквадратур.Пример.Рассмотрим уравнение y = y + ln(y ).Будем искать решение в параметрической форме:y = p + ln(p).685Уравнения, не разрешенные относительно производнойdy= p и дифференцируя предыдущееС учетом того, чтоdxравенства, получим дифференциальное уравнение для нахождения x(p):dy11dp dpdx ==1+dp =+ 2,ppppp1+ C.pТаким образом, решение данного уравнения имеет вид1x = ln(p) − + C,p y = p + ln(p).следовательно, x = ln(p) −5.3.

Особые точки и особые линииОпределение. Точку области G, в которой рассматривается дифференциальное уравнение (2.2) (или (2.5)), будем называть обыкновенной точкой уравнения (2.2) или (2.5), если существует такая окрестность этой точки, что через каждую точкуэтой окрестности проходит ровно одна интегральная кривая и,кроме того, по крайней мере одна из правых частей уравнений(2.2) или (2.5) непрерывна в этой окрестности.Определение. Особой точкой дифференциального уравнения называется точка, не являющаяся обыкновенной точкой.Определение. Линия, любая точка которой — особая, называется особой линией.Если особая линия есть интегральная кривая, она называется особой интегральной кривой.Пример.dyy · ln2 |y|, если y = 0;=0,если y = 0;dx– это интегральная кривая неединственности.5.4 Особое решение69В конкретных примерах обычно наибольший интерес представляет разыскание интегральных линий неединственности,так как их знание помогает представить картину интегральных кривых в целом.5.4.