close

Вход

Забыли?

вход по аккаунту

?

7 Дифф уравнения (2)

код для вставкиСкачать
 7. РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
7.1. Задача Коши
Обыкновенное дифференциальное уравнение порядка n имеет вид:
(7.1)где )- производные первого, второго, ... , n -го порядков от искомой функции y. Его решением является семейство функций
y = y(x,a1, a2, ... , an),
где a1, a2, ... , an - произвольные константы.
Например, простейшее дифференциальное уравнение имеет решение y = aex (рис.7.1). Каждому какому угодно значению параметра a соответствует своя функция, и все эти функции удовлетворяют исходному уравнению.
Рис.7.1. Семейство кривых - решение дифференциального уравнения y' = y
Рис.7.2. Результат численного решения задачи Коши Если в дополнение к уравнению (7.1) задать конкретные значения
для некоторого значения x0 в виде
; ; ;...;,
(7.1')то тем самым определяется конкретный набор a1, a2, ... , an и, следовательно, единственная конкретная функция y(x,a1, a2, ... , an) из всего семейства решений.
Условия (7.1') называются начальными условиями, а вся задача, включающая дифференциальное уравнение (7.1) и начальные условия (7.1'), называется "задачей Коши".
К сожалению, класс дифференциальных уравнений, позволяющих аналитическими методами получить решение, довольно узок. Например, уравнение y' = x2 + y2 не имеет аналитического решения. В большинстве практических задач функция F или коэффициенты, входящие в нее, могут содержать существенные нелинейности или даже задаваться в виде таблиц экспериментальных данных, и тогда аналитическое решение задачи Коши становится невозможным.
При численном решении задачи Коши необходимо задаваться границами xнач, xкон изменения аргумента x и величиной h, являющейся шагом его изменения, который определяет дискретность вычисления значений функции y = y(x).
Решение, полученное численным методом, есть таблица соответствующих значений (xi,yi), i = 0,1,2,...,n, где x0=xнач, xn= xкон; xi+1 = xi + h (см.рис.7.2).
Поскольку численное решение задачи Коши широко применяется в различных областях науки и техники, то оно в течение многих лет было объектом пристального внимания и число разработанных для него методов очень велико. Здесь мы остановимся на следующих двух группах методов решения задачи Коши.
1. Одношаговые методы, в которых для нахождения каждой новой точки на кривой y=y(x) требуется информация лишь об одном предыдущем шаге. Одношаговыми являются метод Эйлера и методы Рунге-Кутта.
2. Методы прогноза и коррекции (многошаговые методы), в которых для отыскания каждой следующей точки кривой y=y(x) требуется информация более чем об одной из предыдущих точек. К этой группе относятся методы Адамса, Гира и т.д.
7.2. Погрешности
Прежде чем перейти к обсуждению конкретных методов численного решения дифференциальных уравнений, остановимся на источниках погрешностей, связанных с численной аппроксимацией. Таких источников три:
1. Погрешность округления обусловлена ограничениями на представление чисел в используемой ЭВМ, так как для любой из них число значащих цифр, запоминаемых и используемых в вычислениях, ограниченно.
2. Погрешность ограничения (усечения) связана с тем, что для аппроксимации функции вместо бесконечных рядов часто используются лишь несколько первых ее членов. Это обычный для численных методов прием, являющийся источником погрешностей, целиком обусловленных применяемым методом и не зависящих от характеристик самой ЭВМ.
3. Погрешность распространения является результатом накопления погрешностей, появившихся на предыдущих этапах счета. Так как ни один приближенный метод не может дать совершенно точных результатов, то любая возникшая в ходе вычислений погрешность сохраняется и на последующих стадиях счета.
Указанные три источника погрешностей являются причиной наблюдаемых ошибок двух типов:
1. Локальная ошибка - сумма погрешностей, вносимых в вычислительный процесс на каждом шаге вычислений.
2. Глобальная ошибка - разность между вычисленным и точным значением величины на каждом этапе реализации численного алгоритма, определяющая суммарную погрешность, накопившуюся с момента начала вычислений.
7.3. Одношаговые методы решения задачи Коши
Одношаговые методы рассмотрим на примере решения обыкновенного дифференциального уравнения первого порядка вида
y' = f (x,y) ,(7.2)при начальном условии
y(x0) = y0.(7.2') С помощью этих методов вычисляют последовательные значения y, соответствующие дискретным значениям независимой переменной x.
7.3.1. Решение с помощью рядов Тейлора
Предположим, что нами уже найдено приближенное решение уравнения (7.2) для точек x0, x1, x2, . . . , xm. При этом последовательные значения xi расположены на расстоянии h друг от друга, т.е. xi+1 = xi + h.
Разложим искомую функцию y(x) в ряд Тейлора в окрестности точки xm:
,
где - значение j -й производной от функции y(x), вычисленное в точке x = xm.
Найдем приближенное значение ym+1 для точки xm+1, подставив в это разложение вместо x величину xm+1 :
(7.3)Чем больше членов этого ряда мы возьмем для вычисления, тем точнее будет решение.
Из (7.2) имеем: . Дифференцируя обе части (7.2) по x и учитывая, что y есть функция от x, получаем:
или для сокращения записи: , где fx, fy - частные производные от функции по x и y соответственно.
Тогда выражение (7.3) приобретает вид:
(7.4)где O(h3) означает, что в следующие (отброшенные) члены ряда значение h входит в степени не ниже третьей.
Таким образом, если для решения уравнения (7.2) будет использована формула (7.4), то погрешность усечения будет приблизительно равна Ch3, где C - некоторая постоянная, не зависящая от h.
Решение дифференциального уравнения данным способом является одноступенчатым, так как для вычисления каждого ym+1 требуется информация только об одной предыдущей точке (xm, ym).
С практической точки зрения трудность использования этого метода заключается в необходимости нахождения и вычисления частных производных fx, fy, что в некоторых ситуациях бывает просто невозможно. Кроме того, если попытаться получить лучшее приближение, то необходимо вычислять уже третью производную:
.
Производные более высоких порядков становятся еще более сложными. На практике этот метод не используется, а здесь он приведен как основа для вывода других методов и оценки их погрешностей.
7.3.2. Метод Эйлера
Это простейший метод решения задачи Коши, позволяющий интегрировать дифференциальные уравнения первого порядка. Его точность невелика, и поэтому на практике им пользуются сравнительно редко. Однако на основе этого метода легче понять алгоритмы других, более эффективных методов.
Итак, решается задача Коши (7.2, 7.2'). Запишем разложение (7.3) для m=0, отбросим в нем члены, содержащие h во второй и более высоких степенях, и получим:
.
(7.5)Величину находим из дифференциального уравнения (7.2), подставив в него начальное условие: . Таким образом можно получить приближенное значение зависимой переменной при малом смещении h от начальной точки.
Этот процесс можно продолжить, используя соотношение
(7.6)и делая сколь угодно много шагов. Графически метод Эйлера показан на рис.7.3. Хотя тангенс угла наклона касательной к истинной кривой в исходной точке известен и равен y'(x0), он изменяется в соответствии с изменением независимой переменной. Поэтому в точке x0+h наклон касательной уже не таков, каким он был в точке x0. Следовательно, при сохранении начального наклона касательной на всем интервале [x0,x1] в результаты вносится погрешность. Ошибка метода имеет порядок h2, а сам метод является методом первого порядка, так как в его вычислительной формуле (7.6) параметр h имеет максимальную степень -1.
Рис.7.3. Геометрическая интерпретация метода Эйлера
Рис.7.4. Ошибка метода Эйлера на m-м шаге На рис.7.5 представлена блок-схема решения задачи Коши по методу Эйлера
Рис.7.5. Алгоритм метода Эйлера
7.3.3. Усовершенствованный метод Эйлера Точность метода Эйлера можно существенно повысить, улучшив аппроксимацию производной. Используются два способа такой аппроксимации. Первый называется "усовершенствованным" методом Эйлера, второй - "модифицированным" методом Эйлера.
Геометрическая интерпретация усовершенствованного метода приведена на рис.7.6.
Рис.7.6. Усовершенствованный метод Эйлера
Прямая L1 есть касательная к истинной кривой y=y(x) в точке (xm, ym). Ее наклон к оси OX равен углу , для которого
,
или в силу (7.2):
.
Прямая L2 есть касательная к решению уравнения (7.2) в точке , являющейся пересечением L1 c прямой x = xm+1. Наклон L2 равен углу , для которого
. Прямая проходит через точку , а ее угол наклона равен , для которого
или .
Прямая L параллельна и проходит через точку (xm, ym), а ее пересечение с x = xm+1 как раз и определяет окончательное значение ym+1.
Данное геометрическое построение в аналитическом виде выглядит следующим образом: сначала вычисляется значение функции y(x) в точке xm+1 по методу Эйлера:
,(7.7)а затем оно используется для вычисления , т.е. приближенного значения производной в конце интервала (xm,xm+1):
.
Вычислив среднее между этим значением производной и ее значением f(xm, ym) в начале интервала, найдем более точное значение ym+1:
.
(7.8) Формула (7.8) представляет собой вычислительный алгоритм усовершенствованного метода Эйлера.
Принцип, на котором основан усовершенствованный метод Эйлера, можно пояснить
и иначе. Для этого в разложении (7.3) функции в ряд Тейлора в окрестности точки xm сохраним член, содержащий h2, и отбросим члены более высоких порядков:
(7.9)Однако, чтобы сохранить член с h2 надо знать вторую производную y"(xm). Ее можно аппроксимировать конечной разностью :
(7.10)Подставив это выражение в ряд Тейлора (7.9) и заменяя его приближением по методу Эйлера, найдем
(7.11)что совпадает с полученным в (7.8) выражением.
Этот метод является методом второго порядка, так как в нем используется член ряда Тейлора, содержащий h2. Ошибка метода на каждом шаге имеет порядок h3.
На рис.7.7 приведена блок-схема алгоритма усовершенствованного метода Эйлера.
Рис.7.7. Алгоритм усовершенствованного метода Эйлера
7.3.4. Модифицированный метод Эйлера
В усовершенствованном методе Эйлера усреднялись наклоны касательных, т.е. производные от искомой функции. В модифицированном методе происходит усреднение точек (рис.7.8).
Прямая L1 есть касательная к истинной кривой y=y(x) в точке (xm, ym). Ее наклон к оси OX равен углу , для которого
,
или в силу (7.2):
.
Рис.7.8. Модифицированный метод Эйлера
Прямая L2 есть касательная к решению уравнения (7.2) в точке , являющейся пересечением L1 c прямой x = xm+h/2. Наклон L2 равен углу , для которого
. Прямая L параллельна прямой L2 и проходит через точку (xm, ym), а ее пересечение c прямой x = xm+1 и определяет окончательное значение ym+1 решения уравнения в тоске xm+1. Уравнение прямой L можно записать в виде
,
где .
Поэтому
(7.12)Выражение (7.12) есть вычислительная формула модифицированного метода Эйлера. Он также согласуется с разложением в ряд Тейлора с точностью до h2. Блок-схема этого алгоритма аналогична предыдущей и отличается лишь формулой в блоке "ордината следующей точки".
7.3.5. Метод Рунге-Кутта
При использовании "исправленного" или "модифицированного" методов Эйлера по сравнению с методом Эйлера для получения каждой новой точки приходится вычислять
значение функции f(x,y) уже дважды - в точках () и () или () и (), т.е. за повышение точности приходится расплачиваться дополнительными затратами машинного времени.
Более высокая точность может быть достигнута, если пользователь готов потратить дополнительное машинное время на еще лучшую аппроксимацию производной путем сохранения большего числа членов ряда Тейлора. Эта идея лежит в основе методов Рунге-Кутта. Собственно, усовершенствованный и модифицированный методы Эйлера являются методами Рунге-Кутта 2-го порядка.
Чтобы удержать в ряде Тейлора член n-го порядка, необходимо каким-то образом вычислить n -ую производную зависимой переменной. При использовании модификаций метода Эйлера для получения второй производной в конечно-разностной форме достаточно было знать наклоны кривой на концах рассматриваемого интервала. Чтобы вычислить третью производную в конечно-разностном виде, необходимо знать значение второй производной по меньшей мере в двух точках. Для этого необходимо дополнительно определить наклон кривой в некоторой промежуточной точке интервала h, т.е. между xm и ym. Oчевидно, чем выше порядок вычисляемой производной, тем больше дополнительных вычислений потребуется внутри интервала. Метод Рунге-Кутта дает набор формул для расчета координат внутренних точек, требуемых для реализации этой идеи. Так как существует несколько способов расположения внутренних точек и выбора относительных весов для найденных производных, то метод Рунге-Кутта в сущности объединяет целое семейство методов решения дифференциальных уравнений первого порядка.
Наиболее распространенным из них является метод, при котором удерживаются все члены, включая h4. Это метод четвертого порядка точности, для которого ошибка на шаге имеет порядок h5. Расчеты при использовании этого классического метода производятся по формуле
,(7.13)
гдеK0 = h f(xm, ym), K1 = h f(xm+0.5h, ym+0.5K0), K2 = h f(xm+0.5h, ym+0.5K1),K3 = h f(xm+h, ym+K2). Блок-схема алгоритма метода Рунге-Кутта приведена на рис.7.10.
Рис.7.10. Блок-схема метода Рунге-Кутта 4-го порядка
По сравнению с методами Эйлера и его модификациями метод Рунге-Кутта четвертого порядка имеет важное преимущество, так как обеспечивает более высокую точность, которая с лихвой оправдывает дополнительное увеличение объема вычислений.
В таблице 7.1. приведены результаты решения дифференциального уравнения
при начальных условиях y(1.0)=1.0. Вторая графа содержит точки полученного табулированием точного, аналитического, решения
Следующие четыре графы содержат результаты численного решения методами, соответственно, Эйлера, усовершенствованного - Эйлера, модифицированного - Эйлера, Рунге-Кутта четвертого порядка. В последней графе помещены ошибки метода Рунге-Кутта четвертого порядка по сравнению с точным решением.
Результаты решения дифференциального уравнения численными методами
Таблица 7.1.
x
Точное
решение
Метод
ЭйлераУсовершенствованный
ЭйлераМодифици-
рованный
ЭйлераМетод
Рунге-
КуттаОшибка
метода
Рунге-Кутта1.01.00001.00001.00001.00001.00000.00000001.10.99610.99270.99600.99580.99610.00000101.20.99800.99190.99780.99770.99800.00000111.31.00490.99631.00461.00471.00490.00000121.41.01641.00541.01611.01661.01640.00000131.51.03251.01891.03221.03321.03250.00000131.61.05321.03691.05281.05451.05320.00000141.71.07841.05921.07781.08041.07830.00000151.81.10771.08571.10711.11061.10770.00000151.91.14061.11611.13991.14461.14060.00000162.01.17641.14961.17551.18161.17640.00000162.11.21381.18531.21271.22041.21380.00000172.21.25151.22201.25021.25951.25150.00000172.31.28771.25811.28621.29731.28770.00000182.41.32051.29211.31891.33181.32050.00000182.51.34821.32191.34651.36111.34820.00000192.61.36891.34591.36701.38311.36890.00000192.71.38101.36221.37901.39611.38100.00000192.81.38331.36951.38121.39891.38330.00000192.91.37511.36661.37301.39041.37510.00000183.01.35631.35321.35431.37061.35630.0000018
Из таблицы хорошо видно, что метод Рунге-Кутта четвертого порядка дает самую высокую точность.
Более высокая точность метода Рунге-Кутта часто позволяет увеличить шаг интегрирования h. Допустимая погрешность на шаге определяет его максимальную величину. Чтобы обеспечить высокую эффективность вычислительного процесса, величину h следует выбирать именно из соображений максимальной допустимой ошибки на шаге. Такой выбор часто осуществляется автоматически и включается как составная часть в алгоритм, построенный по методу Рунге-Кутта.
7.3.6. Метод Кутта-Мерсона
Мерсон предложил модификацию метода Рунге-Кутта четвертого порядка, позволяющую оценивать погрешность на каждом шаге и принимать решение об изменении величины шага. Схема Мерсона выглядит следующим образом:
(7.14)
гдеK1 = h3 f(xm, ym), h3=h/3,K2 = h3 f(xm+h3, ym+K1), K3 = h f(xm+h3, ym+(K1+K2)/2),K4 = K1+ h3 f(xm+h/2, ym+0,375(K1+K3)),K5 = h3 f(xm+h, ym+1,5(K4 - K3)).
Эта схема требует на каждом шаге вычислять правую часть дифференциального уравнения в пяти точках, но она позволяет на каждом шаге определять погрешность решения R по формуле
R = 0,1(2K4 - 3K3 - K5).(7.15) Для автоматического изменения шага интегрирования рекомендуется следующий критерий. Если абсолютное значение величины R, вычисленное по формуле (7.15), на (m+1)-м шаге окажется больше допустимой заранее заданной погрешности , т.е. , то шаг h уменьшается вдвое и вычисления по схеме (7.14) повторяются с точки (xm,ym). При выполнении условия 32 шаг h можно удвоить начиная с точки (xm+1,ym+1).
Следует обратить внимание, что, если по условиям задачи требуется сохранять в памяти ЭВМ все вычисленные точки до конца решения, то, по сравнению с другими методами, здесь необходимо организовывать массив и для абсцисс точек, т.к. шаг изменения по оси OX - переменный.
На рис. 7.11. в виде блок-схемы представлен алгоритм решения задачи Коши для ОДУ первого порядка методом Кутта-Мерсона.
Рис.7.11. Блок-схема метода Кутта-Мерсона
7.3.7. Применение одношаговых методов для решения обыкновенных
дифференциальных уравнений высоких порядков
Преобразуем дифференциальное уравнение (7.1) n-го порядка к системе n дифференциальных уравнений 1-го порядка
(7.16)Запись уравнения (7.1) в виде системы (7.16) называется формой Коши. Начальные условия (7.1') при таком преобразовании выглядят так:
z1(x0)=z1,0 ; z2(x0)=z2,0 ; z3(x0)=z3,0 ; . . . ; zn(x0)=zn,0 ;(7.16')
Для примера преобразуем к форме Коши уравнение Бесселя :
.
Обозначим искомую функцию y(x) через z1(x), а ее первую производную - z2(x). Тогда получим систему уравнений первого порядка, эквивалентную исходному уравнению: Вычислительный алгоритм "усовершенствованного" метода Эйлера для задачи Коши (7.16,7.16') выглядит аналогично (7.8):
(7.17)где i=1,2,...,n ; ; k=1,2,...,n .
Вычислительный алгоритм "модифицированного" метода Эйлера для задачи Коши (7.16,7.16') выглядит аналогично (7.12):
(7.18)где i=1,2,...,n ; ; k=1,2,...,n .
Вычислительная схема метода Рунге-Кутта четвертого порядка для задачи Коши (7.16,7.16') имеет вид:
,(7.19)
гдеK i,0 = h fi(xm, z 1,m, z 2,m,..., z n,m), K i,1 = h fi(xm+0.5h, z 1,m+0.5K1,0, z 2,m+0.5K2,0,..., z n,m+0.5Kn,0), K i,2 = h fi(xm+0.5h, z 1,m+0.5K1,1, z 2,m+0.5K2,1,..., z n,m+0.5Kn,1),K i,3 = h fi(xm+h, z 1,m+K1,2, z 2,m+K2,2,..., z n,m+Kn,2);i=1,2,...,n. Аналогичным образом обобщается на случай системы уравнений и схема (7.14) Кутта-Мерсона.
7.3.8. Общая характеристика одношаговых методов
Всем одношаговым методам присущи определенные общие черты:
1) Чтобы получить информацию в новой точке, надо иметь данные лишь в одной предыдущей точке. Это свойство можно назвать "самостартованием", поэтому одношаговые методы в литературе иногда называются самостартующими.
2) В основе всех одношаговых методов лежит разложение функции в ряд Тейлора, в котором сохраняются члены, содержащие h в степени до k включительно. Целое число k называется порядком метода. Погрешность на шаге имеет порядок k+1.
3) Все одношаговые методы не требуют действительного вычисления производных - вычисляется лишь сама функция, однако могут понадобиться ее значения в нескольких промежуточных точках. Это влечет за собой дополнительные затраты времени и усилий.
4) Свойство "самостартования" позволяет легко менять величину шага интегрирования h.
7.4. Методы прогноза и коррекции
Эти методы в отличие от рассмотренных выше являются многоступенчатыми. Алгоритмы таких методов основываются на аппроксимации интерполяционными полиномами либо правой части обыкновенного дифференциального уравнения f(x,y), либо интегральной кривой y=y(x). Ниже рассматриваются два таких метода: метод Адамса и метод Гира.
7.4.1. Метод Адамса
Суть метода Адамса состоит в следующем: одним из одношаговых методов, например, методом Рунге-Кутта четвертого порядка строится искомая интегральная кривая y=y(x) в нескольких точках x1,x2,...,xn, а затем к ним применяются методы аппроксимации для получения решений в новых точках, лежащих как внутри промежутка [x1,xn] - (интерполяция), так и за его пределами (экстраполяция).
Рассмотрим четырехточечный вариант метода Адамса для задачи Коши (7.2),(7.2').
С помощью любого из одношаговых методов вычислим решения y1,y2,y3 заданного уравнения в точках x1,x2,x3. Правая часть (7.2) - функция f(x,y) на интегральной кривой, соответствующей начальному условию (x0,y0), будет, очевидно, функцией только одного аргумента x: f(x,y) = f(x, y(x)) = f(x), значения которой в рассматриваемых точках обозначим f0, f1, f2, f3 (f0 - значение f(x,y) в точке (x0,y0), заданной начальными условиями (7.2')). В окрестности узлов x0,x1,x2,x3 функцию f(x) приближенно заменим интерполяционным полиномом Ньютона
f(x)= f0+ f01(x-x0)+f012(x-x0)(x-x1)+f0123(x-x0)(x-x1)(x-x2),(7.20)где f01, f012, f0123 - разделенные разности.
Представим искомое решение y4 в точке x4=x3+h в виде тейлоровского разложения около точки x3:
,
(7.21)где - производные по x от правой части исходного уравнения в точке x3.
Дифференцируя полином (7.20) получим выражения для производных:
(x)=
f01+f012(x-x0+x-x1)+f0123[(x-x0)(x-x1)+(x-x0)(x-x2)+ (x-x1) (x-x2)]=
= f01+f012(2x-x0-x1)+f0123[3x2-2x(x0+x1+x2)+ x0x1+x0x2+x1x2];(x)=
2f012+2f0123(3x-x0-x1-x2);(x)=
6f0123.Эти соотношения при x=x3 в случае равноотстоящих узлов принимяют вид:
=(2f0+9f1-18f2+119f3) / h ;=(-f0+4f1-5f2+2f3) / h2 ;
(7.22)=(-f0+3f1-3f2+f3) / h3 . Подставляя производные (7.22) в разложение (7.21) получим экстраполяционную формулу Адамса:
(7.23)имеющую четвертый порядок точности.
Изменяя количество членов, учитываемых в разложении (7.21), можно получить формулы Адамса различных порядков.
Остаточный член формулы (7.23) равен
.
Значительная величина коэффициента ( ) в этом остаточном члене обусловлена тем, что точка x4 лежит вне интервала расположения узлов, по которым построен полином Ньютона. Т.е. здесь мы имеем дело с экстраполяцией, погрешность которой всегда больше, чем погрешность интерполяции. С целью уменьшения погрешности можно получить аналогичным способом, но для узлов x1,x2,x3,x4 интерполяционную формулу Адамса :
(7.24)Эта формула является неявной, т.к. искомая величина y4 необходима для вычисления значения функции f4=f(x4,y4) правой части дифференциального уравнения. Выражение (7.24) можно рассматривать как нелинейное уравнение относительно неизвестной y4 и решать его одним из методов решения трансцендентных уравнений. Обычно здесь используется метод простых итераций, т.к. уравнение (7.24) уже имеет необходимую для этого метода форму. При этом в качестве начального приближения берется значение y4, определенное по экстраполяционной формуле (7.23).
Формулу (7.23) называют формулой прогноза, а формулу (7.24) - формулой коррекции. Эти две формулы без труда переносятся на дифференциальные уравнения высоких порядков, записанные в форме Коши.
7.4.2. Метод Гира
Одним из одношаговых методов получим решения y1,y2,y3 задачи Коши (7.2,7.2') в точках x1,x2,x3. В окрестности узлов x0,x1,x2,x3,x4 искомое решение y(x) приближенно заменим интерполяционным полиномом Ньютона четвертой степени:
y(x)=
+y0+y01(x-x0)+y012(x-x0)(x-x1)+
y0123(x-x0)(x-x1)(x-x2)+y01234(x-x0)(x-x1)(x-x2)(x-x3),
(7.25)где y01, y012, y0123, y01234 - разделенные разности порядков с первого по четвертый.
Левую часть уравнения (7.2), т.е. производную y'(x), приближенно найдем путем дифференцирования по x полинома (7.25):
y'(x)=
+
-y01+y012(2x-x0-x1)+y0123[3x2-2x(x0+x1+x2)+ x0x1+x0x2+x1x2]+
y01234[4x3-3x2(x0+x1+x2)+2x(x0x1+x0x2+x0x3+x1x2+x1x3+x2x3)-
x0x1x2-x0x1x3-x0x2x3-x1x2x3];
(7.26) Разделенные разности для равноотстоящих узлов выражаются через узловые значения аппроксимируемой функции:
y01=(y1-y0) / h,y012=(y2-2y1+y0) / (2h2),(7.27)y0123=(y3-3y2+3y1- y0) / (6h3),y01234=(y4-4y3+6y2-4y1+y0) / (24h4). Полагая в (7.26) x=x4 и учитывая (7.27), получим:
y'(x4)=(3y0-16y1+36y2-48y3+25y4) / (12h).(7.28) C другой стороны, исходное дифференциальное уравнение (7.2) при x=x4 принимает вид:
y'(x4)=f(x4,y4).(7.29) Приравнивая правые части (7.28) и (7.29), находим:
y4=[3(4hf(x4,y4)-y0)+16y1-36y2+48y3] / 25.(7.30) Формула (7.30) представляет собой неявную схему Гира четвертого порядка для решения задачи Коши (7.2,7.2'). Выражение (7.30) есть уравнение относительно y4, для решения которого можно применить метод простых итераций. Начальное приближение к y4 можно получить из следующих соображений. Полагая в выражении (7.26) x=x3, имеем:
y'(x3)=(-y0+6y1-18y2+10y3-y4) / (12h).(7.31)Приравнивая правые части (7.2) при x=x3 и выражения (7.31), получим так называемую схему прогноза
y4=4hf(x3,y3)+(y0-10y3)/3-2y1+6y2,(7.32)которую и можно использовать в качестве начального приближения для решения уравнения (7.30).
7.4.3. Общая характеристика методов прогноза и коррекции
По сравнению с одношаговыми методами методы прогноза и коррекции имеют ряд особенностей:
1. Для реализации методов прогноза и коррекции необходимо иметь информацию о нескольких предыдущих точках: другими словами, они не относятся к числу "самостартующих" методов. Для получения исходной информации приходится прибегать к какому-либо одношаговому методу. Если в процессе решения дифференциальных уравнений методом прогноза и коррекции изменяется шаг, то обычно приходиться временно переходить на одношаговый метод.
2. Поскольку для методов прогноза и коррекции требуются данные о предыдущих точках, то соответственно предъявляются и повышенные требования к объему памяти ЭВМ.
3. Одношаговые методы и методы прогноза и коррекции обеспечивают примерно одинаковую точность результатов. Однако вторые в отличие от первых позволяют легко оценить погрешость на шаге. По этой причине, пользуясь одношаговыми методами, величину шага h обычно выбирают несколько меньше, чем это, строго говоря, необходимо, и поэтому методы прогноза и коррекции оказываются более эффективными.
4. При применении метода Рунге-Кутта четвертого порядка на каждом шаге приходиться вычислять четыре значения функции, в то время как для обеспечения сходимости метода прогноза и коррекции того же порядка точности часто достаточно двух значений функции. Поэтому методы прогноза и коррекции требуют почти вдвое меньше машинного времени, чем методы Рунге-Кутта сравнимой точности. Это обстоятельство может оказывать сильное влияние на выбор алгоритма, так как стоимость машинного времени может быть весьма высокой.
7. Решение обыкновенных дифференциальных уравнений
45
Документ
Категория
Без категории
Просмотров
110
Размер файла
4 068 Кб
Теги
уравнения, дифф
1/--страниц
Пожаловаться на содержимое документа