Оценка ошибки численных методов для решения дифференциальных уравнений второго порядка тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Введение

1 Оценка глобальной ошибки метода Нумерова.

1.1 Случай линейного дифференциального уравнения.

1.2 Случай нелинейного дифференциального уравнения.

2 Оценка глобальной ошибки метода Штермера

2.1 Оценка для явного метода Штермера.

2.1.1 Оценивание вспомогательных величин.

2.2 Оценка для неявного метода Штермера.

3 Метод эллипсоидов для оценки глобальной ошибки метода Штермера

3.1 Оценка для явного метода Штермера.

3.2 Оценка для неявного метода Штермера.

1. Постановка задачи. Данная работа посвящена получению строгой оценки ошибки численного решения задачи у" = 1(х,у), 0 < х < Х0, 2/(0) = 2/0, у'т = Уо, (0-1) где у, / 6 В.п. Задачи вида (0.1), где правая часть дифференциального уравнения не зависит от первой производной, встречаются в механике и в квантовой теории рассеяния. Случай, когда задача (0.1) имеет колеблющееся решение (при п = 1 это имеет место, когда < 0), представляет собой наибольший интерес, так как тогда при использовании известных методов оценки ошибки численного решения возникает так называемый эффект раскрутки (эффект Мура [12]), то есть оценка ошибки растет экспоненциально с ростом длины отрезка, в то время как сама ошибка растет значительно медленнее. Предложенные в работе методы позволяют одновременно с получением приближенного решения строить его окрестность, содержащую точное решение, то есть указывать границы, в которых заключено точное решение.

2. О различных методах решения дифференциальных уравнений. Задачу Коши (0.1) можно решать с помощью формулы Тейлора ( [1], стр.358), то есть в качестве численного решения ут, апроксимирующего точное решение у(х) в точке хт, брать г=0 где хт е [0, Х0] при т — 1,2,., М. Значения у(3)(0), у(4)(0) и т.д. можно получать последовательно диффериенцируя /(ж,у(ж)) соответствующее число раз и подставляя значения в точке (0,2/о) предыдущих производных.

Если Хо больше радиуса сходимости ряда Тейлора, то погрешность метода не стремится к нулю при п стремящемся к бесконечности и предложенный метод неприменим. Но его можно применять для нахождения первых нескольких значений ут при использовании многошаговых численных методов.

Кроме того, можно применять этот метод на каждом шаге. То есть находить приближенное значение в точке хт следующим образом: п „(О

ЕУт-1 / \

• | 1)) г=0 где Ут-1 ~ значения производных в точке решения исходного дифференциального уравнения, проходящего через точку (жтх, ут На практике этот метод применяется относительно редко, поскольку при его использовании нужно вычислять значения всех производных функции /(х,у(х)) до (п — 2)-го порядка включительно.

Следующий способ решения задачи (0.1) основан на сведении уравнения второго порядка к системе уравнений первого порядка У для численного решения которой наиболее часто употребимыми являются конечно-разностные методы: к к = (0.2) г=0 1=0

Здесь У = (у,у)Т, Е{х,У) = (г;, /(ж,у))Т (Т - знак транспонирования). Предполагается, что а0 ^ 0. В случае Ь0 = 0 метод является явным и значение Ут находится из (0.2) через значения Ут-1, Ут~2,—Хт-к• Если же 60 ф 0, то Ут находится из (0.2) с помощью итераций. Начальное приближение для Ут определяется из какой-либо явной формулы. Метод (0.2) начинает работать при т > к, поэтому начальные значения Ух,.надо находить од-ношаговым методом с меньшим шагом к или с помощью формулы Тейлора.

Порядком численного метода (0.2) называется такое число р > 1, что невязка к к Ь[У(х)] = ^бЬ{У(хт{) г=0 г=0 является величиной 0(кр+1). Это накладывает следующие ограничения на коэффициенты метода (см. [1],стр.377): ^г^ + б^г*-1^ = 0, 5 = 1,2,., р. (0.3) г'=0 г=0 г=0

Вопрос о сходимости численного решения к точному при уменьшении шага и увеличении точности начальных значений рассматривался в [20]. А в 1956 году Дальквистом [21] было доказано, что необходимым и достаточным условием сходимости разностного метода (0.2) порядка р > 1 (т.е. (0.3) верно при р > 1) является устойчивость оператора Ь\У(х)]. Это значит, что многочлен

PÍO = S a«'£m * не имеет корней вне единичного круга, причем на единичной окружности может лежать только один корень.

На практике из всех сходящихся методов вида (0.2) употреб ляются, как правило, только методы Адамса:

Недостатком сведения к системе уравнений первого порядка является то, что при этом удваивается количество уравнений, что увеличивает количество вычислений. Кроме того, метод Адамса имеет сравнительно большую локальную ошибку. Поэтому методы, приспособленные специально для решения уравнений второго порядка, в правую часть которых не входит первая производная, часто более эффективны.

Традиционно наиболее распостраненными методами интегрирования уравнений вида у" = f(xfy) являются (см. [32]) явный к к к

0.4)

0.5)

1=1 и неявный к-1

0.6) г=0 конечно-разностные методы. Здесь т > к , Х{ = ih.

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

Nm = у{хт) - 2y(xmi) + у(хт-2) - Ь? 2/(®m-t)) г была величиной как можно более высокого порядка по h. Здесь у{х) - точное решение уравнения (1), а величина Nm называется локальной ошибкой метода на m-ом шаге.

Явный метод вида (0.5) принято называть методом Штермера. Одним из наиболее употребительных неявных методов (0.6) является метод четвертого порядка, называемый методом Нумерова.

На практике равенства (0.5) и (0.6) выполняются.с некоторой погрешностью wm, которая в случае (0.5) включает в себя только ошибку округления, а в случае (0.6) еще и ошибку от обрыва процесса итераций. О различных приемах, с промощью которых можно уменьшить ошибку округления, рассказывается в [33], гл.1 и 3. В частности, формулу явного метода Штермера (0.5) лучше использовать в следующем виде: к-1

Угп = Ут-1 + Vym-1 + Л2/(®т-1,Ут-1 г=2 а формулу (0.6) в виде: к-1

Ут = Ут-1 + Vут-г + /г2/(жт-1, Ут-l) + h2 AV7(хт, ут). i=2

Причем, начинать суммирование надо с членов большего порядка малости. Значения коэффициентов fa приведены в [32] и [33].

В статье [22] приводятся результаты работы исследовательского проекта LONGSTOP ( LONg term Gravitational Stability Test for the Outer Planets ). Проект посвящен вопросу накопления ошибок на больших временных интервалах при численном решении задач небесной механики. Для подсчета численного решения использовались методы Штермера, а в качестве критерия точности результатов использовался АГ2е - критерий. То есть результат считался в достаточной степени достоверным, если произведение квадрата числа шагов на машинную точность не являлось величиной много больше единицы.

Кроме того существуют параметрические семейства гибридных методов типа Нумерова, которые используются в случаях, когда решение задачи (0.1) является периодическим.

Например, в статье [17] приводится семейство М4(а,/3) неявных четырехступенчатых методов четвертого порядка: Уп+2 = 2Уп+1 ~Уп + й2/п+1, У*п = Уп + №(Гп+2- 2/п+1+/п),

Уп+1 = Уп+1 + а^2(/га+2 " 2/п+1 + /*), Уп+2 = 2уп+1 -Уп+ у^Ь2{Гп+2 + ю/*+1 + /п), где/* = /(х{,у*).

При а = 0 и (3 = 0 мы получаем метод Нумерова, для которого начальное приближение найдено методом Штермера второго порядка.

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

У" = ~™2У, (0.7) где IV - константа. Для данного численного метода оптимальными являются а = 1/300, (3 = 1/56. Такой выбор значений параметров увеличивает точность метода на два порядка при решении уравнения (0.7).

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

Например, рассмотрим задачу двух тел х х2 + у2* У х2 + у2 с точным решением x{t) = cost, y(t) = sint, когда одно тело движется вокруг другого по окружности. Здесь численное значение правой части системы в точке tm близко к (—x(tm), —y(tm))T, з поскольку \fx{tm)2 + y{tm)2 близко к единице. Но, хотя система и близка к линейной, численные решения, полученные методами М4(1/300,1/56) и М4(0,0) имеют один порядок точности.

3. Методы оценки ошибки. Первые грубые оценки ошибки разностных методов появились в 30-ых годах. В 1923 году в [23] была доказана сходимость метода Адамса. Уравнение для ошибки zm+1 = y(xm+i) — ут+1 получалось вычитанием уравнения для численного решения из уравнения, которому удовлетворяет точное решение. Откуда, в результате взятия всех членов по абсолютной величине, следует неравенство jfe-i

Zm+l\ < \zm\ + MhY^\zm-i\ + MXhk+1. i=1

Здесь константа М выбрана так, что < М, где - коэффициенты метода Адамса, а величина М\ зависит от ошибки метода и ошибки округления на одном шаге.

Таким образом, оценку ошибки > можно получать последовательно по формуле т+1 = ¿т + Мкк8т + Мхкк+1.

В качестве начальных значений нужно взять <£0 > 1^0 (>•••> > \zk-il, где < ¿1 < . < 6к~ 1

Было показано, что при г > к справедлива оценка

Эта оценка ошибки имеет экспоненциальный рост, так как на отрезке при х = %К

Она является сильно завышенной и не достигается ни на каких задачах.

Более точные оценки ошибки, полученные поздже в [25], [30] для явного метода Адамса и в [26] для неявного, имеют экспоненциальный рост с меньшим показателем степени. В 1938 году в [28] была приведена оценка ошибки с показателем экспоненты равным произведению длины отрезка интегрирования на константу Липшица, доказательство которой была дополнено в [29] в 1951 году. Эта оценка достигается на классе уравнений у' = /(х, у) с > 0. В 1957 году в [31] подобный результат был доказан для системы дифференциальных уравнений первого порядка.

М1

1 + Мкку-к+1 ~ (1 +

Мгкк.1 —:—) ~ е г

МИгк

В 1934 году в [27] была получена грубая оценка ошибки метода Штермера. Предлагалось на каждом шаге находить оценку ошибки через оценки, полученные на предыдущих шагах: к = 2^-1 " + Ы + Ят, г=1 где Ь - константа Липшица правой части дифференциального уравнения по 2/, а включает в себя оценку ошибки метода и ошибки округления на одном шаге. Эта оценка ошибки метода Штермера имеет экспоненциальный рост и является сильно завышенной в случае < 0.

В 1955 году Н.С.Бахвалов [7] получил для одного дифференциального уравнения оценку ошибки, использующую оценку щ не по абсолютной величине. В случае, когда < 0 эта оценка ошибки имеет линейный рост или является ограниченной.

Для метода Адамса в случае одного уравнения в [15] А.Ф.Филипповым описан метод построения верхней и нижней границы точного решения. Для получения оценки решения сверху предлагается в формулу метода Адамса добавлять константу с избытком компенсирующую ошибку метода Адамса и ошибку округления. Величина этой постоянной подсчитывается на основе уточненного выражения локальной ошибки. То есть по формуле ут = + + + 1пт (0.8) г на каждом шаге находится ут > у{хт). Аналогично получается оценка решения снизу.

Этот метод в [13] и [14] О.Б. Ермаковым обобщается на случай системы дифференциальных уравнений, решаемых явным и неявным методом Адамса соответственно, причем на систему накладывается условие > 0 при г Ф 3. Предлагается многомерные интервалы [у¿] такие, что у» £ [г/»], находить по формуле

Ут] = [;Ут-1] + Кя + XI [Ут-<])) + г где ^ - постоянный вектор, а [г/тг]) - интервал значений функции /.

Однако, в силу накладываемых ограничеий на правую часть системы применить предложенный метод к задаче (0.1) при п — 1 в интересующем нас случае, когда < 0, нельзя.

При покоординатном оценивании численное решение заключается в движущийся параллелипипед, грани которого параллельны координатным плоскостям. В случае, когда решение задачи колеблется, возникает эффект раскрутки (эффект Мура [12]).

Этого можно избежать, если [16] заключать численное решение в произвольный параллелипипед. То есть параллельно с подсчетом численного решения на каждом шаге нужно производить подсчет параллелипипеда, содержащего точное решение. Этот метод не приводит к эффекту раскрутки, однако, для оценки ошибки удобнее использовать эллипсоиды, поскольку для задания эллипсоида нужно меньшее число параметров, чем для задания параллелипипеда.

В статье [18] предлагается сочетать метод Адамса (при уточненной по сравнению с [19] оценкой локальной ошибки и с учетом ошибки округления) с методом эллипсоидов, применявшимся в [9] для другой цели. Точное решение у(х) 6 К" заключается в движущийся эллипсоид, центром которого является численное решение.

Уравнение, определяющее движение эллипсоида с ростом ж, выводится с учетом ошибок метода и ошибок округления как при решении самой задачи (0.1), так и при решении уравнения изменения эллипсоидов. Это позволяет одновременно с получением численного решения строить его окрестность, содержащую точное решение. Недостатком метода эллипсоидов является его громоздкость, так как на каждом шаге приходится решать дополнительно 2 ' уравнение.

Ошибку численного решения, полученного специальным методом для решения дифференциальных уравнений второго порядка, можно оценить покоординатно, как предложенно в [11]. Для этого ^-шаговый метод в пространстве К1 представляется в виде одношагового метода в пространстве Як, и производится оценка каждой координаты в отдельности. Эта оценка ошибки имеет экспоненциальный рост, даже тогда, когда сама ошибка численного решения растет значительно медленнее (например, в случае уравнения у" = —ги2 у).

4. Асимптотическая оценка ошибки метода. Получить представление об асимптотическом поведении ошибки гт = у(хт) — ут можно следующим образом [33]: предположим, что при /г стремящемся к нулю ошибка начальных данных и ошибка округления на каждом шаге являются величинами 0(Д5+1), в то время как сама ошибка решения является величиной порядка к3. То есть гт = г(хт)к3 + 0(к3+1). (0.9)

Тогда при малых к справедливо ^ г{хт)к3. Множество таких к называется областью асимптотики. Однако, в реальных условиях ошибка начальных данных не зависит от шага к, а ошибка округления при фиксированной разрядной сетке ограничена по абсолютной величине снизу. Поэтому реальная область асимптотики, когда главный член (0,9) хорошо представляет полную ошибку приближенного решения, ограничена не только сверху, но и снизу.

Рассмотрим к такое, что к и к/2 принадлежат области реальной асимптотики. Пусть ^ и ^ - значения приближенного решения в точке жт, вычисленного с шагами к и к¡2 соответственно. Тогда к

У{Хт) -Уш~ Фт)^, У{хт) ~ & « г(хт)(-)'.

Исключая из этих равенств значение точного решения у{хт), получим

Откуда

Следовательно, для ошибок приближенных решений у^ и у2т справедливо следующее

Ут Ут

1 „,(Т \ ~ у™ У

Хт ~ У\Хт) - Ут--1 2~ г1 = у(хт)-у1~у1

2* - 1

Если мы хотим получить численное решение с ошибкой по абсолютной величине не превосходящей €, то предлагается искать подходящее значение шага интегрирования следующим образом: положим

Подставив значение z(xm) из (0.10), получим

О.И)

2 У \Ут ~ Ут\

Видно, что если 12> €, то новое значение шага уменьшается: he < h/2, а если < е, то увеличивается: he> h. Таким образом, (0.11) дает более подходящее значение шага интегрирования.

Недостатком предложенного метода является то, что он не позволяет получать численное решение с гарантированной точностью. Кроме того, уменьшение шага в некоторых случаях может приводить к выходу за границу области реальной асимптотики.

5. Содержание работы. Диссертация состоит из введения, трех глав, заключения и списка литературы. Предлагаемые в данной работе методы строгой оценки ошибки численного решения ут основаны на построении разностного уравнения, которому удовлетворяет ошибка.

Заключение

Предложенные в работе способы оценки ошибки метода Штер-мера позволяют одновременно с получением приближенного решения строить его окрестность, содержащую точное решение. При решении систем дифференциальных уравнений второго порядка ошибка оценивается с помощью эллипсоидов. В случае одного дифференциального уравнения предлагается применять менее трудоемкий метод, основанный на использовании оценок двух линейно независимых решений уравнения в вариациях.

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

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

1. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М., 1987.

2. Гелъфонд А.О. Исчисление конечных разностей. М., 1967.

3. Золотарева Н.Д. Строгая оценка ошибки численного решения линейного дифференциального уравнения второго порядка. Вестн.Моск.Ун-та.сер. 15, вычисл. матем. и киберн. 1999.N1.

4. Золотарева Н.Д. Строгая оценка ошибки численного решения задачи Коши для нелинейного дифференциального уравнения второго порядка. Вестн.Моск.Ун-та.сер. 15, вычисл. матем. и киберн. 2000.N1.

5. Золотарева Н.Д. Строгая оценка ошибки метода Штермера. Вестн.Моск.Ун-та.сер. 15, вычисл. матем. и киберн. 2001.N2.

6. Бахвалов Н.С. К оценке ошибки при численном интегрировании дифференциальных уравнений экстраполяционным методом Адамса. ДАН СССР, 1955, том 104, N5, 683-686.

7. Березин И. С., Жидков Н.П. Методы вычислений. Т.2 М., Физ-матгиз, 1960.

8. Черноусъко Ф.Л. Оценивание фазового состояния динамических систем. Москва "Наука", 1988.

9. Решетняк Ю.Н. Суммирование эллипсоидов в задаче гарантированного оценивания. ПММ, 1989, т.53, N2, 249-254.

10. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. М., 1990.

11. Moore R.E. Interval analysis. Prentice-Hall, Englewood Cliffs. NJ, 1966.

12. Ermakov O.B. Two-sided method for solving system of ordinary differential equations with automatic determination of guaranteed estimates. Interval Computations, 3(5),1992,pp.63-69.

13. Ermakov O.B. Solving systems of ordinary differential equations using Adams' interpolation method with guaranteed accuracy. Interval Computations, 1,1994,pp.90-95.

14. Филиппов А.Ф. Получение на ЭВМ строгих оценок решений дифференциальных уравнений. Журнал вычислительной математики и математической физики, том 31, 1991, с.994-1005.

15. Cambill T.N., Skeel R.D. Logarithmic reduction of the wrapping effect with application to ordinary differential equation. SIAM Journ. Numer. Anal., 1988, (25), N1, 153-162.

16. Franco J.M. An explicit hybrid method of Numerov type for second-order periodic initial-value problems. Journal of Computational and Applied Mathematics 59 (1995) 79-90.

17. Filippov A.F. Ellipsoidal error estimates for Adams method. Interval Computations, N 3(5), 1992.

18. Stewart N.F. Computable, guaranteed local error bounds for the Adams method. Machr.,1974, v.60, N 1-6 , 145-153.

19. Коллатц JI. Численные методы решения дифференциальных уравнений. М., И.Л., 1953.

20. Dahlquist G. Convergence and stability in the numerical integration of ordinary differential equations. Math. Scand.4 (1956), 33-53.

21. Milani A., Nobili A. Integration error over very long time spans. Celestial Mechanics,43 (1988), 1-34.

22. Tamarkine J. Math. Zeitschr., 1923, (16), 214-219.

23. Степанов B.B. Курс дифференциальных уравнений, 1968.

24. Щ Mises R. ZAMM, 1930, (10), 81-92.

25. Schulz G. ZAMM, 1932, (12), 44-59.

26. Schulz G. ZAMM, 1934, (14), 224-234.

27. Tollmien W. ZAMM, 1938, (18), 83-90.

28. Mohr E. Math. Nachr., 1951, (5), 209-218.

29. Weissinger 3. ZAMM, 1950, (30), 356-363.