Разностная схема для решения обыкновенного дифференциального уравнения

О построении разностных схем

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

Пример простейшей прямоугольной области G(z, у) с границей Г в двумерном случае показан на рис. 2.1. Стороны прямоугольника делятся на элементарные отрезки точками и . Через эти точки проводятся два семейства координатных прямых х = const и у = const, образующих сетку с прямоугольными ячейками. Любой ее узел, номер которого (i,j), определяется координатами (xi, yj). Поскольку все ячейки показанной на рис. 2.1 сетки одинаковы, такую сетку называют равномерной сеткой.

Рис. 2.1. Прямоугольная сетка

Аналогично вводятся сетки для многомерных областей, содержащих более двух измерений. На рис. 2.2 показан элемент сетки в виде прямоугольного параллелепипеда для трехмерной области.

Рис. 2.2. Элемент сетки

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

Узлы сетки, лежащие на границе Г области G, называются граничными узлами. Все остальные узлы — внутренними. Поскольку начальные и граничные условия при постановке задач формулируются на границе расчетной области, то их можно считать заданными в граничных узлах сетки. Иногда граничные точки области не являются узлами сетки, что характерно для областей сложной формы. Тогда либо вводят дополнительные узлы на пересечении координатных линий с границей, либо границу приближенно заменяют ломаной, проходящей через близкие к границе узлы. На эту ломаную переносятся граничные условия.

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

Рис. 2.3. Преобразование расчетной области

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

Вобласти G‘ можно ввести прямоугольную сетку, при этом в области Gей будет соответствовать сетка с неравномерно расположенными узлами и криволинейными ячейками.

В дальнейшем при построении разностных схем мы для простоты будем использовать прямоугольные сетки (или с ячейками в виде прямоугольных параллелепипедов в трехмерном случае), а уравнения будем записывать в декартовых координатах (х, у, z). На практике приходится решать задачи в различных криволинейных системах координат: полярной, цилиндрической, сферической и др. Например, если расчетную область удобно задать в полярных координатах (r,φ), то в ней сетка вводится с шагами и , соответственно, по радиус-вектору и полярному углу.

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

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

В качестве примера построим некоторые разностные схемы для решения уравнения теплопроводности при заданных начальных и граничных условиях. Запишем смешанную краевую задачу в виде

(2.2)

где — начальное распределение температуры U(при t= 0); , — распределение температуры на концах рассматриваемого отрезка (х = 0, 1) в любой момент времени t. Заметим, что начальные и граничные условия должны быть согласованы, т.е. .

Введем равномерную прямоугольную сетку с помощью координатных линий , ; hи τ — соответственно шаги сетки по направлениям х и t. Значения функции в узлах сетки обозначим . Эти значения заменим соответствующими значениями сеточной функции , которые удовлетворяют уравнениям, образующим разностную схему. Часто верхний индекс заключают в скобки, чтобы не путать его с показателем степени. Здесь и далее скобки для краткости опушены.

Заменяя в исходном уравнении (2.2) частные производные искомой функции с помощью отношений конечных разностей, получаем разностную схему

. (2.3)

В записи этой схемы для каждого узла использован шаблон, изображенный на рис. 2.4, а.

Для одного и того же уравнения можно построить различные разностные схемы. В частности, если воспользоваться шаблоном, изображенным на рис. 2.4, б,т.е. аппроксимировать производную при :

,

то вместо (2.3) получим разностную схему

(2.4)

И в том и другом случае получается система алгебраических уравнений для определения значений сеточной функции во внутренних узлах. Значения в граничных узлах находятся из граничных условий:

(2.5)

Совокупность узлов при t = const, т. е. при фиксированном значении j, называется слоем (или, поскольку переменная tсоответствует времени, временным слоем).Схема (2.3) позволяет последовательно находить значения на -ом слое через соответствующие значения на j-ом слое. Такие схемы называются явными.

Для начала счета по схеме (2.3) при j= 1 необходимо знать решение на начальном слое при j= 0. Оно определяется начальным условием (2.2), которое запишется в виде:

(2.6)

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

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

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

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

Разностная схема третьего порядка точности для решения жесткого обыкновенного дифференциального уравнения с линейными коэффициентами Текст научной статьи по специальности « Математика»

Аннотация научной статьи по математике, автор научной работы — Зверев В. Г.

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

Похожие темы научных работ по математике , автор научной работы — Зверев В. Г.

A third order difference scheme for solution of a first order stiff ordinary differential equation with linear coefficients

A new implicit difference scheme for numerical analysis of the stiff Cauchy problem for a first order ordinary linear differential equation is proposed. This scheme is constructed using the increased accuracy Taylor expansion of function in the vicinity of the desired solution and the integration of the differential equation. In the case of linear coefficients the obtained difference scheme has the third order of approximation. Good rate of convergence to the exact solutions is shown on test examples in a wide range of a grid parameter. Comparison with other known one-step methods is carried out.

Текст научной работы на тему «Разностная схема третьего порядка точности для решения жесткого обыкновенного дифференциального уравнения с линейными коэффициентами»

РАЗНОСТНАЯ СХЕМА ТРЕТЬЕГО ПОРЯДКА

ТОЧНОСТИ ДЛЯ РЕШЕНИЯ ЖЕСТКОГО ОБЫКНОВЕННОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ С ЛИНЕЙНЫМИ КОЭФФИЦИЕНТАМИ

В. Г. Зверев Томский государственный университет, Россия e-mail: zverev@niipmm.tsu.ru

A new implicit difference scheme for numerical analysis of the stiff Cauchy problem for a first order ordinary linear differential equation is proposed. This scheme is constructed using the increased accuracy Taylor expansion of function in the vicinity of the desired solution and the integration of the differential equation. In the case of linear coefficients the obtained difference scheme has the third order of approximation. Good rate of convergence to the exact solutions is shown on test examples in a wide range of a grid parameter. Comparison with other known one-step methods is carried out.

Задача Коши для обыкновенного дифференциального уравнения (ОДУ) первого порядка является математической моделью широкого круга прикладных проблем в различных областях науки. В качестве примера могут быть названы механика жидкости и газа, химическая кинетика, биология и т. д. Классы уравнений, для которых возможны точные решения в элементарных функциях, крайне узки и не охватывают всего многообразия возникающих на практике ситуаций. Поэтому наиболее общим подходом к получению решения ОДУ является применение численных методов [1, 2].

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

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

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

В качестве примера можно назвать известную в практике газодинамических расчетов неравновесных течений полидисперсных сред полностью неявную разностную схему второго порядка точности [14, 15]. Особенностью ее построения является использование отрезка ряда Тейлора для вычисления зависимой переменной в полуцелом узле через искомое решение с заменой первой производной соответствующим выражением из дифференциального уравнения. Использование этого приема позволило обеспечить полностью неявную форму разностного аналога, второй порядок аппроксимации и правильную асимптотику численного решения в случае малых значений параметра при производной, что очень важно для задач с пограничными слоями. Учитывая эти обстоятельства, представляют интерес дальнейшее повышение точности и расширение возможностей схем такого класса.

Цель настоящей работы — в развитие исследований [16, 17] разработка неявного экономичного А-устойчивого метода повышенного порядка точности для численного интегрирования жесткой задачи Коши для обыкновенного линейного дифференциального уравнения первого порядка.

1. Математическая постановка задачи и анализ предельных случаев

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

Ьп(х) = е-—+ а(х)и(х) = f (х), и(0) = и0, х Е (0,Х). (1)

Здесь е Е (0,1] — малый параметр; а(х), f (х) — гладкие функции, причем а(х) > ао > 0 отделена от нуля для всех х > 0. Задача (1) имеет единственное решение [8, 9], согласно принципу максимума для него имеет место априорная оценка, не зависящая от е:

||и(х)||те Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

В предельном случае при г — 0 из (19) следует ^¿+1 — (//а)^+1, что отвечает правильной асимптотике.

Если в (13) ограничиться разложением второго порядка, то данный способ построения приводит к схеме вида

иг + кг^-1 [/г+1/2 + /г+1 (г/2)] ^¿+1 =—-. (20)

¿+1 1 + £¿+1/2 + ¿¿+1^/2 У У

Сравнивая (20) с (8), видим, что при одинаковой структуре записи схема (20) имеет отличие в способе усреднения коэффициента а при старшей степени к2. Любопытно, что разностный аналог (8) также может быть получен на основе выражения (15). Для этого в интеграле при М^+1 функцию а(£) на [х^,х^+1 ] следует заменить приближенным значением

Наконец, если в (15) ограничиться двумя первыми интегралами, для вычисления которых использовать формулу правых прямоугольников, то получим неявную схему Эйлера первого порядка

иг+1 = (и + кгг-1/г+1)/[1 + ¿¿+1]. (21)

4. Анализ упрощенных вариантов схемы

Рассмотрим случай кусочно-постоянного коэффициента а(£) = а и линейного источника /(£). Точное решение задачи Коши в этом случае можно представить в виде

иг+1 = иге(г) + кгг 1[/г+1^(г) + /¿п(г)],

Рис. 4. Сеточные функции и их аппроксимации: £1^), — первый, £2^), —

второй, пэ(^) — третий порядки точности.

е(а) = exp(—а), £(а) = [а — 1 + ехр(-а)]/а2 п(а) = [1 — ехр(-а)(а + 1)]/а2

£(а) + п(а) = в (а) = [1 — ехр(-а)]/а, а = аЛ^е

где е(а), £(а), п(а), в (а) — сеточные функции [16, 17]. Их графики показаны на рис. 2-4. Функция е(а) является решением однородного уравнения (1), £(а), п(а) играют роль весовых множителей при источнике на (г + 1)-м и г-м слоях.

Нетрудно видеть, что схема (19) соответствует записи (22) при следующем представлении сеточных функций:

Отметим, что для схем второго порядка (8), (20) эти выражения имеют вид б2(а)= [1 + а + а2/2]-1, £2(а) = 0.5[1 + а] б2(а), П2(а) = 0.5в2(а).

На рис. 2-4 хорошо видно, что приближение (23) значительно точнее, чем (24), воспроизводит сеточные функции точного решения (22) во всем диапазоне значений аргумента, особенно при малых и больших а. Асимптотический анализ (22) показывает, что аппроксимации (23) следуют из точного решения (22) при использовании разложения экспоненты в ряд Тейлора до третьего порядка точности.

Для кусочно-постоянного источника точное решение (22) переходит в (11) и разностная схема (19) принимает вид

1 + а + а2/2 + а з/6

что соответствует следующей аппроксимации сеточных функций е(а), в (а): ез(а) = [1 + а + а 2/2 + аз/6]-1, вз(а) = [1 + а/2 + а2/6]вз (а).

5. Результаты численных расчетов и их анализ

Результаты численного решения рассмотренной задачи 3 с помощью предложенных схем (20) и (19) при е = 0.1 и грубом шаге Л = 0.25 показаны на рис. 1. Для сравнения приведены данные по неявной схеме Эйлера (21) первого порядка. С теоретической точки зрения схемы (8) и (20) дают одинаковые результаты, качественно и количественно превышающие возможности неявной схемы Эйлера. Однако они уступают результатам, полученным по схеме (19) третьего порядка аппроксимации. Ее применение позволило сократить относительную ошибку численного решения на первом грубом расчетном шаге до

12 % и с приемлемой для практики точностью воспроизвести особенности формирования (кривая 2) пограничного слоя. В числе других преимуществ схемы (19) — высокая экономичность, так как используются только рациональные выражения. Необходимо заметить, что более трудоемкая разностная схема специального вида, предложенная в [16, 17], решает эту задачу точно для любых значений Л и е.

Рис. 5. Решение задачи (27) при различных е: кривые 1-3 соответствуют е = 1; 0.1; 0.01; расчеты: Д — метод Эйлера [12], V — схема [12], * — схема (8) [14, 15], □ — схема (20), о — схема (19); е = 1, Н = 1; е = 0.1, Н = 0.25.

Для более детальной оценки возможностей предложенных разностных схем (20) и (19) рассмотрим тестовую задачу Коши с переменным коэффициентом а(х) [12]:

Аналитическое решение п(х) при различных значениях е показано на рис. 5. Для сравнения здесь в соответствии с работой [12] приведены результаты по неявной схеме Эйлера и равномерной по е схеме (11) первого порядка. Видно, что полученные при грубом шаге по схеме (19) данные графически близки к кривым 1, 2 точного решения. Более подробно ошибка Ди численного решения пи при различных значениях к и е указана в таблице:

Ее анализ показывает, что, несмотря на иной способ усреднения коэффициента а(х) в (20), ошибки по схемам (8) и (20) второго порядка незначительно отличаются друг от друга. Схема (19) третьего порядка по сравнению со схемами (8) и (20) имеет гораздо более низкий уровень ошибки во всем рассматриваемом диапазоне значений к и е.

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

В основе построения схемы лежат разложение функции в ряд Тейлора в окрестности искомого решения и прямое интегрирование дифференциального уравнения.

ей'(х) + (1 + х)п(х) = (1 + х), п(0) = 0, х € [0, 2], п(х) = 1 — ехр(—(2х + х2)/(2е)).

Ди = тах |п(х) — пи(х)|, х € [0, 2].

Погрешность численного решения задачи (27) при различных значениях Н и е

Схема (8) [14, 15] Схема (20) Схема (19)

h £ = 1 £ = 0.1 £ = 0.01 £ = 1 £ = 0.1 £ = 0.01 £ = 1 £ = 0.1 £ = 0.01

1 2.7e-2 6.0e-3 6.6e-5 3.8e-2 6.7e-3 7.4e-5 4.1e-3 1.0e-3 1.2e-6

0.1 6.2e-4 3.1e-2 1.4e-2 8.1e-4 3.2e-2 1.5e-2 2.0e-5 6.2e-3 3.6e-3

0.01 6.8e-6 5.4e-4 3.2e-2 8.9e-6 5.7e-4 3.2e-2 2.3e-8 1.2e-5 7.0e-3

0.001 6.9e-8 5.8e-6 5.7e-4 9.0e-8 6.1e-6 5.7e-4 2.4e-11 1.3e-8 1.4e-5

0.0001 6.9e-10 5.9e-8 6.1e-6 9.0e-10 6.2e-8 6.1e-6 2.5e-14 1.3e-11 1.5e-8

Тестовыми расчетами подтверждена хорошая для практики сходимость численных результатов к точным решениям, в том числе и при малых значениях параметра при производной.

[1] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.

[2] КАлиткин Н.Н. Численные методы. М.: Наука, 1978.

[3] РАкитский Ю.В., Устинов С.М., Черноруцкий И.Т. Численные методы решения жестких систем. М.: Наука, 1979.

[4] Евсеев Г.А., КАлюжный В.В. Экономичный метод численного интегрирования уравнений химической кинетики // Численные методы механики сплошной среды: Сб. науч. тр. АН СССР. ВЦ; ИТПМ. 1974. Т. 5, № 3. С. 21-28.

[5] Пирумов У.Г., Росляков Г.С. Газовая динамика сопел. М.: Наука, 1990.

[6] Васильева А.Б., Бутузов В.Ф. Асимптотические разложения сингулярно возмущенных уравнений. М.: Наука, 1973.

[7] Ломов С.А. Введение в общую теорию сингулярных возмущений. М.: Наука, 1981.

[8] Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.

[9] Багаев Б.М., ШАйдуров В.В. Сеточные методы решения задач с пограничным слоем. Ч. 1. Новосибирск: Наука, 1998.

[10] Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems. Berlin: Springer-Verlag, 1991.

[11] Roos H.G., Stynes M., Tobiska L. Numerical Methods for Singularly Perturbed Differential Equations. Berlin: Springer-Verlag, 1996.

[12] Титов В.А., Шишкин Г.И. Численное решение задачи Коши для обыкновенного дифференциального уравнения с малым параметром при производной // Численные методы механики сплошной среды: Сб. науч. тр. АН СССР. ВЦ; ИТПМ. 1978. Т. 9, № 7. С. 112-121.

[13] Боглаев И.П. О численном интегрировании сингулярно возмущенной задачи Коши для обыкновенного дифференциального уравнения // Журн. вычисл. математики и мат. физики. 1985. Т. 25, № 7. С. 1009-1022.

[14] ВАСЕНИН И.М., Архипов В.А., Бутов В.Г. и др. Газовая динамика двухфазных течений в соплах. Томск: Изд-во Том. ун-та, 1986.

[15] Рычков А.Д. Математическое моделирование газодинамических процессов в каналах и соплах. Новосибирск: Наука, 1988.

[16] Зверев В.Г. О численном решении сингулярно возмущенной задачи Коши для дифференциального уравнения первого порядка // Вычисл. гидродинамика. Томск: Изд-во Том. ун-та, 1999. С. 73-82.

[17] Зверев В.Г., Гольдин В.Д. Об одной разностной схеме для решения обыкновенного дифференциального уравнения первого порядка с малым параметром при производной // Сопряженные задачи механики и экологии. Томск: Изд-во Том. ун-та, 2000. С. 166-174.

Поступила в редакцию 10 августа 2005 г., в переработанном виде —13 января 2006 г.

Разностные схемы для решения задачи Коши

Рассмотрим основные принципы конструирования разностных схем.

1. Производная [math]\left.<\left(\frac\right)\!>\right|_>[/math] аппроксимируется на заранее выбранном шаблоне [math](x_, x_,\ldots,x_)[/math] , а правая часть [math]f(x,y)[/math] обыкновенного дифференциального уравнения заменяется сеточным представлением [math]f(x_,\widehat_)[/math] в точке [math](x_,\widehat_)[/math] . Условно этот принцип называется аппроксимационным. С его помощью строятся одношаговые и многошаговые схемы в зависимости от шаблона, по которому аппроксимируется производная.

2. Дифференциальное уравнение заменяется разностным аналогом на одношаговом (двухточечном) шаблоне [math](x_,x_)[/math] (возможно использование дополнительных точек внутри шаблона) путем согласования с разложением функции [math]y=y(x)[/math] при [math]x=x_[/math] no формуле Тейлора относительно точки [math]x=x_[/math] . По этому принципу строятся одношаговые схемы, в частности семейство схем Рунге-Кутты.

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

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

Принцип аппроксимаций

Рассмотрим проблему нахождения численного решения задачи Коши (6.9): [math]\frac= f(x,y),

y(x_0)=y_0[/math] . Для аппроксимации производной [math]\left.<\left(\frac\right)\!>\right|_>[/math] используем формулу (5.4), записанную на двухточечном шаблоне [math]H_<2,i>= (x_,x_)\colon[/math]

и формулу (5.10), записанную на трехточечном нерегулярном шаблоне [math]H_<3,i>= (x_, x_, x_)\colon[/math]

\delta_= \frac>>[/math] — параметр нерегулярности.

При [math]\delta_=1[/math] шаблон становится регулярным, а среднее слагаемое с входящим в него [math]y_[/math] — равным нулю. При этом формула (6.17) сводится к традиционной:

Аппроксимация формульных функций, входящих в правую часть (6.9), осуществляется путем перехода к их сеточному представлению.

Методика построения разностных схем с помощью аппроксимаций производной

1. Вводится в общем случае неравномерная сетка [math]\Omega_n= (x_0, x_1,\ldots, x_, x_,\ldots, x_)[/math] . Величина шага [math]h_=x_-x_[/math] выражается через узловые точки.

2. Для получения рассматриваемых здесь двух схем (одношаговой и двух-шаговой) выбирается соответственно двухточечный (одношаговый) и трехточечный (двухшаговый) шаблоны [math]H_<2,i>= (x_,x_)[/math] и [math]H_<3,i>= (x_,x_,x_)[/math] , на которых в точке [math]x_[/math] производится аппроксимация производной [math]\left.<\left(\frac \right)\!>\right|_>[/math] по формулам. (6.16), (6.17) соответственно в левой и центральной точках [math]x_[/math] . Далее заменяется правая часть уравнения (6.9) ее сеточным представлением, т.е. [math]f(x,y)\to f(x_,y_)[/math] , а вместо [math]y(x)[/math] рассматривается сеточная функция [math]\widehat_\approx y(x_)[/math] , которая определяется только в точках сетки.

3. Выполняется подстановка аппроксимаций (6.16), (6.17) и [math]f(x_,y_)[/math] в дифференциальное уравнение:

После отбрасывания остаточных слагаемых получаются разностные схемы: — явная схема Эйлера первого порядка (явный метод Эйлера):

– обобщенная на нерегулярный шаблон двухшаговая явная схема Эйлера второго порядка (где [math]\Delta \widehat_=\widehat_-\widehat_[/math] ):

Обозначим схему (6.20) символами 2Y2B (суть этого обозначения поясняется далее). Анализ остаточного слагаемого [math]\left(\frac^2><6>\delta_M_<3,i>\right)[/math] обобщенной схемы Эйлера указывает на то, что при условном выборе шаблона, когда достигается соотношение [math]\delta_\leqslant h_[/math] (то есть [math]h_\leqslant h_^2[/math] ) , порядок аппроксимации этой схемы повышается на единицу без увеличения количества точек шаблона. Этот прием повышения точности может быть использован при конструировании вложенных алгоритмов путем измельчения шага.

При [math]\delta_=1[/math] (регулярный шаблон) схема (6.20) соответствует методу Эйлера-Коши:

Для начала расчетов по формулам (6.20),(6.21) требуется иметь две «разгонные» точки [math]\widehat_<0>,\,\widehat_<1>[/math] . Первая определяется известным начальным условием [math]\widehat_0=y_0[/math] , а вторая может быть найдена с помощью другого метода, например, по формуле (6.19): [math]\widehat_1=y_0+h_1f(x_0,y_0)[/math] .

Приведем геометрическую интерпретацию явного метода Эйлера (6.19), которая показана на рис. 6.5.

Пусть известна точка [math](x_,\widehat_)[/math] на искомой интегральной кривой [math]y=y(x)[/math] . Через эту точку проведем к кривой [math]y=y(x)[/math] касательную [math]L[/math] , тангенс угла наклона которой равен [math]y'(x_)= f(x_, \widehat_)[/math] . Ее уравнение имеет вид

Следующая точка [math]\widehat_[/math] приближенного решения получается при [math]x=x_[/math] (в результате нахождения точки пересечения прямых [math]L[/math] и [math]x_=\text[/math] ): [math]\widehat_= \widehat_+ h_f(x_,\widehat_)[/math] . Это соотношение совпадает с формулой (6.19).

Схема (6.19) иногда называется схемой ломаных (рис. 6.6), так как очередное значение [math]\widehat_[/math] получается в результате проведения касательной к интегральной кривой в точке [math](x_,\widehat_)[/math] . После выполнения каждого шага метода Эйлера приближенное решение переходит с одной интегральной кривой на другую, что обусловлено погрешностью численного решения. Большая погрешность схемы связана с тем, что порядок этой схемы равен единице, и поэтому ломаная, определяемая по (6.19), довольно значительно отклоняется от точного решения [math]y=y(x)[/math] . В связи с этим данная схема используется, как правило, в качестве вспомогательной при конструировании схем повышенного порядка.

Для анализа устойчивости явного метода Эйлера применим формулу (6.19) при [math]h_=h[/math] к решению задачи (6.12). В результате разностное уравнение будет иметь вид [math]\widehat_= \widehat_+ h\,\mu\,\widehat_[/math] или [math]\widehat_-(1+h\,\mu)\widehat_=0[/math] . Корень соответствующего характеристического уравнения [math]\lambda_1=1+h\,\mu[/math] . Условие (6.13) принимает форму [math]|1+h\,\mu| . Отсюда следует, что область устойчивости представляется внутренней частью круга единичного радиуса с центром в точке [math](-1;0)[/math] (рис. 6.7). Если [math]\mu[/math] действительная константа, то можно указать интервал устойчивости [math]-2 , то есть [math]h\,\mu\in (-2;0)[/math] , а при [math]\mu имеем [math]0 . Поэтому явный метод Эйлера является ограниченно устойчивым с критическим шагом [math]h_<\text>=-\frac<2><\mu>[/math] .

Дадим геометрическую интерпретацию метода Эйлера-Коши (6.21) применительно к решению задачи (6.9) (рис. 6.8). Пусть известны две точки [math](x_,\widehat_),\, (x_, \widehat_)[/math] на искомой интегральной кривой [math]y=y(x)[/math] . Проведем через точку [math](x_,\widehat_)[/math] касательную [math]L[/math] , тангенс угла наклона которой равен [math]y'(x_)= f(x_,\widehat_)[/math] . После этого через точку [math](x_, \widehat_)[/math] проведем прямую [math]L_1[/math] , параллельную [math]L[/math] . Ее уравнение имеет вид [math]y=\widehat_+ f(x_,\widehat_)(x-x_)[/math] . Следующая точка [math]\widehat_[/math] , приближенного решения получается при [math]x=x_\colon[/math] [math]\widehat_= \widehat_+2h f(x_, \widehat_)[/math] . Это соотношение совпадает с формулой (6.21).

Замечание. Все описанные здесь и далее методы легко переносятся на системы обыкновенных дифференциальных уравнений (6.7). Под [math]y[/math] и [math]f(x,y)[/math] следует подразумевать векторы [math]Y=(y_1,\ldots,y_n)^T[/math] и [math]F(x,Y)= \bigl(f_1(x,Y), \ldots, f_n(x,Y)\bigr)^T[/math] . Например, для системы двух дифференциальных уравнений

явный метод Эйлера (6.19) записывается в форме

Пример 6.3. Исследовать устойчивость явного метода Эйлера на примере решения уравнения [math]Ty’+y=1,[/math] [math]y(0)=y_0[/math] , где [math]T[/math] — действительное положительное число. Найти приближенное решение задачи Коши [math]0,\!1y’+y=1,

y(0)=0[/math] на отрезке [math][0;1][/math] явным методом Эйлера.

Точное решение задачи получено в примере 6.1: [math]y(x)=1-e^<-10x>[/math] , поскольку [math]T=0,\!1[/math] . Уравнение [math]Ty’+y=1[/math] перепишем в виде, аналогичном (6.12): [math]y’=-\frac<1>y+\frac<1>[/math] , поэтому [math]\mu=-\frac<1>[/math] , а критический шаг равен [math]h_<\text>=-\frac<2><\mu>=2T[/math] . Следовательно, при [math]h метод устойчив, а при 2T»>[math]h>2T[/math] неустойчив. На рис. 6.9 схематически изображены приближенные решения уравнения [math]Ty’+y=1,[/math] [math]y(0)=0[/math] , соответствующие устойчивому процессу и хорошей точности численного решения (рис.6.9,а), устойчивому процессу и плохой точности (рис.6.9,б), а также неустойчивому процессу и плохой точности (рис. 6.9,в).

Для уравнения [math]0,\!1y’+y=1[/math] постоянная времени [math]T=0,\!1[/math] , поэтому [math]h_<\text>==0,\!2[/math] . Решим задачу для значений шага [math]h=0,\!05; 0,\!2; 0,\!5[/math] . При [math]h=0,\!05[/math] метод является устойчивым, при [math]h=0,\!2[/math] находится на фанице устойчивости, а при [math]h=0,\!5[/math] является неустойчивым.

Переписав уравнение в виде [math]y’=10-10y[/math] , получим из (6.19) формулу для нахождения приближенного решения:

6.1>>\\\hline x_& \widehat_

(h=0,\!5)& y(x_)\\\hline 0,\!00& 0,\!000000& 0,\!000& 0,\!000& 0,\!000000 \\\hline 0,\!05& 0,\!500000& & & 0,\!393469 \\\hline 0,\!01& 0,\!750000& & & 0,\!632121 \\\hline 0,\!15& 0,\!875000& & & 0,\!776870 \\\hline 0,\!20& 0,\!937500& 2,\!000& & 0,\!864665 \\\hline 0,\!25& 0,\!968750& & & 0,\!917915 \\\hline 0,\!30& 0,\!984375& & & 0,\!950213 \\\hline 0,\!35& 0,\!992187& & & 0,\!969803 \\\hline 0,\!40& 0,\!996094& 0,\!000& & 0,\!981684 \\\hline 0,\!45& 0,\!998047& & & 0,\!988891 \\\hline 0,\!50& 0,\!999023& & 5,\!000& 0,\!993262 \\\hline 0,\!55& 0,\!999511& & & 0,\!995913 \\\hline 0,\!60& 0,\!999755& 2,\!000& & 0,\!997521 \\\hline 0,\!65& 0,\!999877& & & 0,\!998497 \\\hline 0,\!70& 0,\!999938& & & 0,\!999088 \\\hline 0,\!75& 0,\!999969& & & 0,\!999446 \\\hline 0,\!80& 0,\!999984& 0,\!000& & 0,\!999665 \\\hline 0,\!85& 0,\!999992& & & 0,\!999797 \\\hline 0,\!90& 0,\!999996& & & 0,\!999877 \\\hline 0,\!95& 0,\!999998& & & 0,\!999925 \\\hline 1,\!00& 0,\!999999& 2,\!000& -15,\!000& 0,\!999955 \\\hline \max\varepsilon_& 0,\!117879& 1,\!135335& 15,\!99995& \\\hline \end[/math]

Результаты расчетов, выполненных для различных значений [math]h[/math] , приведены в табл. 6.1. При [math]h=0,\!05[/math] приближенное решение достаточно близко к точному, при [math]h=0,\!2[/math] «пила» не сходится и не расходится, а при [math]h=0,\!5[/math] «пила», аппроксимирующая точное решение, расходится (рис. 6.9,в).

Замечание. Порядок точности метода, как правило, определяется порядком аппроксимации схем. Например, для схемы (6.19) первого порядка, записанной в виде, соответствующем исходному уравнению (6.9), получаем (первое слагаемое аппроксимирует производную):

то есть [math]y’_=f_+O(h_)[/math] . Здесь были использованы формула Тейлора и свойства, приведенные ранее.

Для двухшаговой схемы (6.21) аналогично получается второй порядок аппроксимации. Действительно, воспользуемся частным случаем формулы Тейлора (см. (В. 18)) при

Подставляя полученные соотношения в схему (6.21), записанную в виде [math]\frac-y_><2h>-f(x_,y_)=0[/math] , соответствующем дифференциальному уравнению (6.9), с учетом свойств, приведенных ранее, имеем

Также легко проверить, что и схема (6.20) на нерегулярном шаблоне имеет второй порядок при безусловной аппроксимации. Если же [math]h_\leqslant h_^2[/math] (условная аппроксимация), то этот порядок равен трем.

Интегрально-интерполяционный принцип

В данном подходе явные или неявные методы соответствующего порядка получаются путем интегрирования уравнения (6.9) вдоль решения [math]y(x)[/math] . Так, на двухточечном шаблоне

Функция [math]f(x,y(x))[/math] в (6.22) заменяется интерполяционным многочленом соответствующей степени [math]p[/math] ( [math]p=0[/math] для метода первого порядка, [math]p=1[/math] для метода второго порядка, [math]p=2[/math] для метода третьего порядка и т.д.). Подчеркнем, что замена [math]f(x,y)[/math] интерполяционным многочленом в (6.22) становится возможной благодаря тому, что [math]f(x,y)[/math] рассматривается не на всей плоскости [math](x,y)[/math] , а только вдоль интегральной кривой [math]y(x)[/math] . Для вычисления интеграла в правой части могут быть использованы известные квадратурные формулы прямоугольников, трапеций, парабол (Симпсона) и др.

Примем [math]x_[/math] в качестве узловой точки многочлена [math]L_<0>(x)[/math] , то есть [math]L_0(x)= f(x_,\widehat_)[/math] . Тогда из (6.22) следует формула неявного метода Эйлера:

Принимая для интерполяционного многочлена [math]L_1(x)[/math] в качестве узловых точек [math]x_[/math] и [math]x_[/math] и используя квадратурную формулу трапеций (5.46), получаем неявную одношаговую схему второго порядка (метод трапеций):

где [math]f_= f(x_,\widehat_)[/math] . Формулы (6.23),(6.24) имеют вид (6.11) при [math]k=1[/math] . Подчеркнем, что свойство неявности схем (6.23),(6.24) обусловлено наличием искомой величины [math]\widehat_[/math] в левой и правой частях в общем случае нелинейных уравнений.

При реализации алгоритма решения задачи Коши неизвестное значение [math]\widehat_[/math] вычисляется одним из методов решения нелинейных уравнений. Применение метода Ньютона связано с записью уравнения в форме

и с дифференцированием функции [math]F(\widehat_)[/math] , что увеличивает время расчетов из-за возможной сложности производных или увеличивает погрешность в силу неточного задания правых частей.

Как правило, используется метод простых итераций:

Таким образом, уравнение (6.24) решается с помощью формулы (рекуррентного соотношения):

При применении методов Ньютона и простых итераций вначале задается или находится нулевое приближение решения по формуле [math]y_^<\,(0)>= \widehat_[/math] , (так называемый «постоянный» прогноз) или явным методом Эйлера (6.19): [math]\widehat_^<\,(0)>= \widehat_+ h_ f(x_, \widehat_)[/math] .

Итерации завершаются при выполнении условия окончания [math]\bigl|\widehat_^<\,(k+1)>-\widehat_^<\,(k)>\bigr| \leqslant \varepsilon[/math] , где [math]\varepsilon[/math] — малое положительное число.

Для сокращения описания и удобства ссылок получаемые здесь схемы обозначаются буквенно-цифровыми символами. Первая цифра в нем обозначает «шаговость» схемы, следующая за ней буква «Я» или буквы «НЯ» — явный или неявный характер, затем цифра указывает порядок точности схемы и за ней буква — возможную модификацию. Таким образом, схема (6.24) обозначается 1НЯ2, т.е. она является одношаговой, неявной и имеет второй порядок точности. Последняя буква указывается только для тех схем, которые имеют модификации.

Приведем геометрическую интерпретацию неявного метода Эйлера (6.23) (рис.6.10). Пусть известна точка [math](x_, \widehat_)[/math] на искомой интегральной кривой [math]y=y(x)[/math] . Через эту точку проведем прямую, тангенс угла наклона которой равен [math]y'(x_)= f(x_, \widehat_)[/math] , т.е. следующая точка [math]\widehat_[/math] приближенного решения берется на той кривой из семейства интегральных кривых, касательная к которой в точке [math](x_,\widehat_)[/math] проходит через точку [math](x_,\widehat_)[/math] . Поэтому [math]\frac<\widehat_-\widehat_>>= f(x_, \widehat_)[/math] . Полученное соотношение совпадает с формулой (6.23).

Исследуем устойчивость неявного метода Эйлера и метода трапеций. Применим формулу (6.23) при [math]h_=h[/math] к решению задачи (6.12). В результате получим разностное уравнение [math]\widehat_= \widehat_+h\cdot\mu\cdot \widehat_[/math] или [math]\widehat_(1-h\cdot\mu)-\widehat_=0[/math] . Корень характеристического уравнения [math]\lambda_1= \frac<1><1-h\cdot\mu>[/math] . Условие устойчивости (6.13) имеет вид [math]\left|\frac<1><1-h\cdot\mu>\right| или 1″>[math]|1-h\cdot\mu|>1[/math] . Область устойчивости представляется внешней частью круга единичного радиуса с центром в точке [math](1;0)[/math] . Из геометрической интерпретации (рис. 6.11) следует, что неявный метод Эйлера обладает свойством A-устойчивости (см. рис.6.4,а), поскольку выделенная на рис. 6.11 область включает левую полуплоскость.

Применим формулу (6.24) метода трапеций при [math]h_=h[/math] к решению задачи (6.12). Получим разностное уравнение

Корень характеристического уравнения [math]\lambda_1= \left(1+ \frac<2>\right)\,\colon \left(1-\frac<2>\right)[/math] . Условие устойчивости (6.13) принимает вид [math]\left|\left(1+ \frac<2>\right)\,\colon \left(1-\frac<2>\right)\right| или [math]|2+m\mu| . Полученное неравенство выполняется в области, изображенной на рис. 6.4,а, т.е. метод трапеций является A-устойчивым.

Продолжим рассмотрение интегрально-интерполяционного метода. Использование квадратурной формулы прямоугольников дает модифицированный метод Эйлера второго порядка:

Дадим геометрическую интерпретацию модифицированного метода Эйлера второго порядка (рис. 6.12). Пусть известна точка [math](x_,\widehat_)[/math] на искомой интегральной кривой [math]y=y(x)[/math] . Через эту точку проведем касательную [math]L[/math] , тангенс угла наклона которой равен [math]y'(x_)= f(x_,\widehat_)[/math] . При [math]x=x_+ \frac><2>[/math] получаем точку [math]\widehat_>\,\,2>[/math] . Далее проведем прямую [math]L_1[/math] с тангенсом угла наклона [math]f\! \left(x_+ \frac><2>, \widehat_>\,\,2>\right)[/math] , а затем прямую [math]L_2[/math] , параллельную [math]L_1[/math] и проходящую через точку [math](x_, \widehat_)[/math] . Ее уравнение имеет вид [math]y=\widehat_+ f\! \left(x_+ \frac><2>, \widehat_>\,\,2>\right)\!\cdot (x-x_)[/math] . Следующая точка [math]\widehat_[/math] приближенного решения находится при [math]x=x_[/math] (это соотношение совпадает с формулой (6.28)):

Исследуем устойчивость модифицированного метода Эйлера второго порядка. Применим формулы (6.28) при [math]h_=h[/math] к решению задачи (6.12). В результате получим разностное уравнение:

Корень соответствующего характеристического уравнения [math]\lambda_1= 1+h\mu+ \frac<(h\mu)^2><2>[/math] . Условие устойчивости (6.13) имеет вид [math]\left|1+h\mu+ \frac<(h\mu)^2><2>\right| . Если [math]\mu[/math] — действительное число, то, раскрывая модуль, получаем [math]-1 . Неравенство 0″>[math]2+h\mu+ \frac<(h\mu)^2><2>>0[/math] выполняется для любых [math]h\mu[/math] , а неравенство [math]h\mu+ \frac<(h\mu)^2> <2>, то есть [math]h\mu\cdot (2+h\mu) , справедливо при [math]-2 . Таким образом, интервал устойчивости [math]h\mu\in (-2;0)[/math] модифицированного метода Эйлера совпадает с интервалом устойчивости явного метода Эйлера (6.19).

Для использования регулярного трехточечного шаблона [math](x_,x_,x_)[/math] заменим (6.22) аналогичным выражением: [math]\textstyle<\int\limits_>^> dy= \int\limits_>^> f(x,y(x))dx>[/math] . Если функцию [math]f(x,y(x))[/math] , входящую в правую часть, заменить многочленом Лагранжа второй степени и применить формулу (5.47), то получается неявная схема парабол (Симпсона) четвертого порядка (2НЯ4):

Отметим, что схема парабол имеет четвертый, а не третий порядок точности из-за регулярности шаблона.

В соответствии с классификацией, схема (6.29) относится к классу неявных схем Адамса. Для получения других схем Адамса следует переписать (6.22) в форме

и заменить [math]y'(x)[/math] соответствующим интерполяционным полиномом, построенным на регулярном шаблоне с [math]h=\text[/math] .

В силу формулы (4.36) второго интерполяционного многочлена Ньютона получаем (где [math]\widehat= \frac>[/math] — фаза интерполяции)

Подставляя полученное выражение в правую часть (6.30) и учитывая, что [math]dx= h\cdot d\widehat[/math] , вычисляем [math]\textstyle<\int\limits_<0>^ <1>y'(\widehat)\cdot h\, d\widehat>[/math] , поскольку при [math]x=x_[/math] и [math]x=x_[/math] справедливо [math]\widehat=0[/math] и [math]\widehat=1[/math] соответственно. В результате имеем

или с учетом (6.9) общую формулу семейства явных схем Адамса:

где [math]\eta_= h\cdot f(x_, \widehat_)= h\cdot f_[/math] , а конечные разности вычисляются с помощью табл. 6.2. Подчеркнем, что здесь интерполяционный многочлен записан на шаблоне, не включающем точку [math]x_[/math] .

Если взять в правой части (6.31) два слагаемых, то получится явный метод Эйлера (6.19), если три слагаемых, то явный метод, описываемый формулой

а если четыре слагаемых, то явный метод вида

В результате имеем схемы Адамса-Бэшфорта

– второго порядка (2Я2В):

– третьего порядка (ЗЯЗА):

– четвертого порядка (4Я4):

Для начала расчетов по формуле (6.32) требуются две «разгонные» точки: [math]\widehat_<0>,\, \widehat_<1>[/math] , по формуле (6.33) — три «разгонные» точки: [math]\widehat_<0>,\, \widehat_<1>,\, \widehat_<2>[/math] , по формуле (6.34) — четыре «разгонные» точки: [math]\widehat_<0>,\, \widehat_<1>,\, \widehat_<2>,\, \widehat_<3>[/math] . Их необходимо вычислить с порядком точности не меньше порядка точности схемы. Далее схема (6.32) обобщается на нерегулярный шаблон.

Если для построения интерполяционного многочлена использовать шаблоны, включающие точку [math]x_[/math] , то получится следующая формула семейства неявных схем Адамса :

где [math]\eta_= h\cdot f(x_, \widehat_)= h\cdot f_[/math] , а конечные разности вычисляются с помощью табл. 6.3.

Если взять в правой части (6.35) два слагаемых, то получится неявный метод Эйлера (6.23), если три слагаемых, — метод трапеций (6.24):

если четыре слагаемых, — неявный метод, описываемый соотношением (и т.д.)

– первого порядка (6.23);

– второго порядка (6.24);

– третьего порядка (2НЯЗ):

– четвертого порядка (ЗНЯ4):

Для расчетов по формуле (6.36) требуются две «разгонные» точки: [math]\widehat_<0>,\, \widehat_<1>[/math] , по формуле (6.37) — три «разгонные» точки: [math]\widehat_<0>,\, \widehat_<1>,\, \widehat_<2>[/math] . Их необходимо вычислить с порядком точности не меньше порядка точности схемы. Чтобы найти искомое значение [math]\widehat_[/math] с помощью формул (6.36),(6.37), так же как в неявном методе Эйлера и методе трапеций, требуется решить в общем случае нелинейное уравнение.

Замечание. Аналогичные методы могут быть получены из интегрального уравнения вида [math]\textstyle<\int\limits_>^> dy= \int\limits_>^> f(x,y(x))\,dx,

k \geqslant 2>[/math] .

Пример 6.4. Найти приближенное решение задачи Коши [math]\beginy’=z-1,&y(0)=1,\\ z’=-y-2z,&z(0)=-1\end[/math] на отрезке [math][0;1][/math] с шагом [math]h=0,\!1[/math] методами Эйлера (неявным и модифицированным), Адамса-Бэшфорта третьего порядка и трапеций.

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

Используя соотношения (6.23),(6.24),(6.28),(6.33), выписываем формулы для нахождения приближенного решения указанными методами (при этом применяется векторная форма записи).

Для неявного метода Эйлера из (6.23) имеем

Разрешая эту систему относительно [math]\widehat_,\, \widehat_[/math] , окончательно получаем

Соотношения (6.28) для модифицированного метода Эйлера принимают вид

Для метода Адамса-Бэшфорта третьего порядка из (6.33) находим

Для определения «разгонных» точек [math](x_0,\widehat_0,\widehat_0),\, (x_1,\widehat_1, \widehat_1),\, (x_2,\widehat_2,\widehat_2)[/math] воспользуемся модифицированным методом Эйлера. Точка [math](x_0,\widehat_0,\widehat_0)[/math] определяется начальными условиями [math](0;1;-1)[/math] .

Выпишем соотношения для метода трапеций (6.24) и разрешим их относительно неизвестных:

откуда путем разрешения системы относительно [math]\widehat_,\, \widehat_[/math] получаем

6.4>>\\\hline x_& \begin\text\\ \text\end & \begin\text\\ \text\end & \begin\text\\ \text \end & \begin\text\\ \text\end & \begin\text\\ \text\end\\\hline 0,\!0&1,\!00000&1,\!00000&1,\!00000&1,\!00000&1,\!00000\\\hline 0,\!1&0,\!80992&0,\!80500&0,\!80499&0,\!80499&0,\!80500\\\hline 0,\!2&0,\!62960&0,\!61997&0,\!61991&0,\!61991&0,\!61994\\\hline 0,\!3&0,\!45885&0,\!44479&0,\!44469&0,\!44464&0,\!44470\\\hline 0,\!4&0,\!29740&0,\!27924&0,\!27908&0,\!279000&0,\!27909\\\hline 0,\!5&0,\!14500&0,\!12309&0,\!12285&0,\!12273&0,\!12286\\\hline 0,\!6&0,\!00131&-0,\!02397&-0,\!02428&-0,\!02444&-0,\!02428\\\hline 0,\!7&-0,\!13397&-0,\!16224&-0,\!16265&-0,\!16284&-0,\!16263\\\hline 0,\!8&-0,\!26119&-0,\!29208&-0,\!29257&-0,\!29279&-0,\!29255\\\hline 0,\!9&-0,\!38072&-0,\!41384&-0,\!41441&-0,\!41465&-0,\!41438\\\hline 1,\!0&-0,\!49287&-0,\!52787&-0,\!52853&-0,\!52879&-0,\!52848\\\hline \max \varepsilon_&0,\!03720415&0,\!00067291&0,\!00005942&0,\!00033550&-\\\hline \end[/math]

6.5>>\\\hline x_& \begin\text\\ \text\end & \begin\text\\ \text\end & \begin\text\\ \text \end & \begin\text\\ \text\end & \begin\text\\ \text\end\\\hline 0,\!0&-1,\!00000&-1,\!00000&-1,\!00000&-1,\!00000&-1,\!00000\\\hline 0,\!1&-0,\!90083&-0,\!90000&-0,\!90023&-0,\!90023&-0,\!90016\\\hline 0,\!2&-0,\!80316&-0,\!80095&-0,\!80132&-0,\!80132&-0,\!80121\\\hline 0,\!3&-0,\!70753&-0,\!70357&-0,\!70402&-0,\!70401&-0,\!70388\\\hline 0,\!4&-0,\!61440&-0,\!60844&-0,\!60893&-0,\!60890&-0,\!60877\\\hline 0,\!5&-0,\!52408&-0,\!51601&-0,\!51650&-0,\!51645&-0,\!51632\\\hline 0,\!6&-0,\!43684&-0,\!42663&-0,\!42709&-0,\!42702&-0,\!42691\\\hline 0,\!7&-0,\!35287&-0,\!34054&-0,\!34095&-0,\!34087&-0,\!34078\\\hline 0,\!8&-0,\!27229&-0,\!25794&-0,\!25828&-0,\!25818&-0,\!25812 1\\\hline 0,\!9&-0,\!19518&-0,\!17894&-0,\!17920&-0,\!17908&-017905\\\hline 1,\!0&-0,\!12158&-0,\!10357&-0,\!10378&-0,\!10364&-0,\!10364\\\hline \max \varepsilon_& 0,\!01958133& 0,\!00032586& 0,\!00012145&0,\!00003005&-\\\hline \end[/math]

Результаты проведенных расчетов даны в табл. 6.4 для [math]\widehat_,\, y(x)[/math] и табл. 6.5 для [math]\widehat_,\, z(x)[/math] соответственно. Из их анализа вытекает, что наиболее точно величину [math]\widehat_[/math] можно рассчитать методом Адамса-Бэшфорта, а величину [math]\widehat_[/math] — методом трапеций.

Пример 6.5. Исследовать устойчивость метода Адамса-Бэшфорта третьего порядка на примере решения задачи Коши

\mu=-1[/math] и [math]\mu=-100[/math] .

Точное решение задачи: [math]y(x)=x^3[/math] . При [math]\mu=-1

\left(h\mu=-\frac<1> <8>\right)[/math] дифференциальное уравнение имеет вид [math]y’=x^3-y+3x^2[/math] , а при [math]\mu=-100

Формула (6.33) при [math]h=1\!\!\not<\phantom<|>>\,8[/math] принимает форму

Положим [math]f_0=y’\! \left(-\frac<1><4>\right)=0,\!1875;

\widehat_0=y_0=-0,\!15625; f_1= y’\! \left(-\frac<1><8>\right)= 0,\!046875;[/math] [math]\widehat_1= y_0=-0,\!001953121;[/math] [math]f_2=y'(0)=0;[/math] [math]\widehat_2=0[/math] , т.е. используем в качестве «разгонных» точек значения, соответствующие точному решению. Результаты нахождения приближенного решения задачи Коши приведены в табл. 6.6.

6.6>>\\\hline i& \widehat_

(h\mu=-100\!\!\not<\phantom<|>>\,8)& \text \\\hline 0& -0,\!01562500& 0,\!00000000& -0,\!01562500& 0,\!00000000 \\\hline 1& -0,\!001953125& 0,\!00000000& -0,\!001953125& 0,\!00000000 \\\hline 2& 0,\!0000000& 0,\!00000000& 0,\!00000000& 0,\!00000000 \\\hline 3& 0,\!00195317& 0,\!00000004& 0,\!00195317& 0,\!00000004 \\\hline 4& 0,\!01562501& 0,\!00000001& 0,\!01562409& 0,\!00000091 \\\hline 5& 0,\!05273429& 0,\!00000009& 0,\!05275568& 0,\!00002130 \\\hline 6& 0,\!12499996& 0,\!00000004& 0,\!12449656& 0,\!00050437 \\\hline 7& 0,\!24414058& 0,\!00000005& 0,\!25408001& 0,\!00993938 \\\hline 8& 0,\!42187492& 0,\!00000008& 0,\!18516620& 0,\!23670880 \\\hline 9& 0,\!66992193& 0,\!00000005& 6,\!27264466& 5,\!60272278 \\\hline \end[/math]

Из анализа решения, приведенного в табл. 6.6, видно, что в случае [math]h\mu=-100\!\!\not<\phantom<|>>\,8[/math] процесс не является устойчивым. Дадим обоснование этому факту.

Запишем формулу (6.33) в виде (6.14):

то есть [math]\alpha_0=1,

Составим уравнение (6.15) при каждом значении [math]h\mu[/math] .

При [math]h\mu=-\frac<1><8>[/math] все корни уравнения [math]\xi^3-\frac<73><96>\xi^2-\frac<16><96>\xi+ \frac<5><96>=0[/math] лежат внутри круга единичного радиуса с центром в начале координат, поэтому процесс устойчив.

При [math]h\mu=-\frac<100><8>[/math] один корень уравнения [math]\xi^3+ \frac<2204><96>\xi^2-\frac<1600><96>\xi+ \frac<500><96>=0[/math] , равный [math]\frac<1><3>\cdot \frac<2204><96>[/math] , очевидно не лежит внутри указанного круга. Поэтому процесс неустойчив. Эти выводы подтверждают результаты численного эксперимента. Ошибка вычисления значении [math]\widehat_[/math] при [math]h\mu=-\frac<100><8>[/math] возрастает (см. правую колонку табл. 6.6).

Принцип согласования с разложением по формуле Тейлора

Рассмотрим задачу Коши (6.9) в общем случае на неравномерной сетке [math]\Omega_[/math] . Для получения формулы нахождения последующего значения [math]\widehat_[/math] на шаблоне [math](x_, x_)[/math] разложим решение [math]y(x)[/math] при [math]x=x_[/math] относительно точки [math]x_[/math] по формуле Тейлора на промежутке [math]x_ \leqslant x \leqslant x_\colon[/math]

Входящие в правую часть производные можно определить, дифференцируя уравнение (6.9) требуемое число раз:

Если функция [math]f(x,y)[/math] имеет p-е непрерывные производные по обоим аргументам, то в разложении (6.38) можно учесть слагаемые до [math]O(h_^)[/math] . Чем больше членов разложения используется для вычисления, тем точнее будет приближение и тем выше порядок схемы. Однако рассчитывать производные путем дифференцирования правой части уравнения (6.9) невыгодно. Поэтому Рунге и Кутта независимо друг от друга предложили более удобную и рациональную форму представления численного решения, при использовании которой дифференцировать функцию [math]f(x,y)[/math] не требуется.

В общем случае методы Рунге-Кутты имеют следующую структуру (где [math][/math] ):

где [math]s[/math] — число стадий (этапов), [math]K_[/math] — значения коэффициентов схемы Рунге-Кутты, вычисленные на основе правой части уравнения (6.9), [math]c_,

\ell= \overline<2,s>;[/math] [math]m=\overline<1,s-1>;[/math] [math]b_k,

k=\overline<1,s>[/math] . Первый индекс в обозначениях коэффициентов является порядковым номером, а второй соответствует индексу точки [math]x_[/math] — началу отрезка [math][x_,x_][/math] , на котором производится расчет.

В некоторых методах кроме вычисления приближенного решения [math]\widehat_[/math] определяется еще дополнительное значение [math]\widetilde_[/math] по формуле

порядок которого, как правило, на единицу больше или меньше обеспечиваемого выражением для [math]\widehat_[/math] . Величина [math]|\widehat_-\widetilde_|[/math] служит для учета погрешности и управления величиной шага.

Коэффициенты схемы подбираются так, чтобы выражение для [math]\widehat_[/math] согласовывалось с формулой Тейлора (6.38) вплоть до членов степени [math]h_^

[/math] , где степень [math]p[/math] меняется в зависимости от метода и соответствует порядку схемы.

Наибольшее распространение в вычислительной практике нашел метод Рунге-Кутты четвертого порядка :

Схема (6.41) является четырехчленной, первый коэффициент [math]K_<1,i>[/math] относится к точке [math]x_[/math] , второй и третий — к средней точке [math]x_+ \frac><2>[/math] четвертый — к точке [math]x_[/math] . Для этой схемы выбираются следующие параметры, входящие в (6.39):

Приведем геометрическую интерпретацию метода Рунге-Кутты четвертого порядка (рис. 6.13). Пусть известна точка [math](x_,\widehat_)[/math] на искомой интегральной кривой [math]y=y(x)[/math] . Прямая [math]AE[/math] , проходящая через точку [math](x_, \widehat_)[/math] и определяющая величину [math]\widehat_[/math] , имеет угловой коэффициент. Этот коэффициент рассчитывается по угловым коэффициентам касательных в точках [math]A,B,C,D[/math] , проведенным к проходящим через эти точки интегральным кривым (штриховые линии). При этом коэффициент [math]K_<1,i>[/math] рассчитывается в точке [math]A[/math] , коэффициент [math]K_<2,i>[/math] — в точке [math]C[/math] , коэффициент [math]K_<3,i>[/math] — в точке [math]B[/math] , а коэффициент [math]K_<4,i>[/math] — в точке [math]D[/math] .

1. Схемы Рунге-Кутты легко обобщаются на системы дифференциальных уравнений. Так, применительно к системе двух уравнений

схема Рунге-Кутты четвертого порядка имеет вид [math](i=\overline<0,n-1>)[/math]

2. Соотношения модифицированного метода Эйлера (6.28) можно записать в форме схемы Рунге-Кутты:

Покажем, что эта схема второго порядка. Разложение (6.38) с учетом равенства [math]y»= f_x+ f_\cdot f_[/math] , где [math]f_= f(x_,y_)[/math] , частные производные [math]f_x, f_y[/math] и вторая производная [math]y»[/math] вычисляются в точке [math](x_, y_)[/math] , принимает вид

Разложение функции [math]f(x,y)[/math] относительно точки [math](x_,y_)[/math] можно записать следующим образом:

где частные производные [math]f_x,f_y[/math] вычисляются в точке [math](x_,y_)[/math] . Тогда при [math]x=x_+\frac><2>,

Следовательно, формулу (6.43) можно переписать в виде

Сравнивая последнее соотношение с (6.44), можно сделать вывод о том, что схема модифицированного метода Эйлера согласуется с разложением по формуле Тейлора вплоть до членов степени [math]h_^2[/math] и поэтому является методом Рунге-Кутты второго порядка.

3. Коэффициенты семейства методов Рунге-Кутты (6.39),(6.40) удобно записывать в виде табл. 6.7. Табл. 6.8 соответствует классическому методу Рунге-Кутты четвертого порядка (6.41), а в табл. 6.9 приведены коэффициенты другого метода четвертого порядка (правило [math]3\!\!\not<\phantom<|>>\,8[/math] ). Табл. 6.10 соответствует методу Рунге–Кутты (6.43) второго порядка, табл. 6.11 и 6.12 — методам третьего порядка. В табл. 6.13 даны коэффициенты метода Фельберга четвертого порядка, где формула для [math]\widetilde_[/math] имеет пятый порядок, а в табл. 6.14 — семистадийного метода Батчера шестого порядка. В табл. 6.15 приведены коэффициенты метода Рунге-Кутты-Мерсона четвертого порядка, где формула для [math]\widetilde_[/math] имеет пятый порядок для систем линейных дифференциальных уравнений с постоянными коэффициентами и третий порядок для нелинейных систем. Параметры метода Рунге-Кутты-Дормана-Принса пятого порядка приведены в табл. 6.16. Метод Фвльберга седьмого порядка и метод Дормана-Принса восьмого порядка изложены.

4. Для s-стадийных методов Рунге-Кутты справедливо следующее условие устойчивости, определяемое для задачи (6.12):

Если [math]\mu[/math] — действительное число, можно указать интервалы устойчивости: [math]h\mu\in (-2;0)[/math] при [math]s=1,\!2;

h\mu\in(-2,\!51;0)[/math] при [math]s=3;

h\mu\in (-2,\!78;0)[/math] при [math]s=4[/math] и т.д.

Пример 6.6. Для задачи Коши [math]y’=x+y,

x\in[0;0,\!2][/math] выполнить два шага методом Рунге-Кутты второго порядка (6.43) и один (первый шаг) методом четвертого порядка (6.41). Шаг интефирования принять постоянным и равным 0,1. Сравнить оба численных решения с точным [math]y(x)= 2e^x-x-1[/math] .

В данной задаче [math]x_0=0,

y_0=1[/math] . Применим схему (6.43).

Первый шаг [math](i=0)[/math] . В соответствии с постановкой задачи принимается [math]h_1=0,\!1[/math] . Вычислим

Тогда [math]\widehat_1= \widehat_0+ h_1K_<2,0>= 1+0,!1\cdot 1,!1= 1,\!11[/math] .

Второй шаг [math](i=1)[/math] . Так как шаг постоянный, то [math]h_2=0,!1[/math] , a [math]x_1=0,\!1[/math] и, следовательно,

Тогда [math]\widehat_2= \widehat_1+ h_2K_<2,1>= 1,\!11+ 0,\!1\cdot 1,\!316= 1,\!2416[/math] . Точное решение: [math]y(x_2)= 1,\!2428055[/math] . Отличие приближенного решения от точного составляет [math]0,\!096\%[/math] .

Выполним теперь первый шаг методом четвертого порядка по формуле (6.41), записанной при [math]i=0\colon[/math]

При этом вычислим [math]K_<1,0>= x_0+ y_0=0+1=1[/math] ;

Точное решение: [math]y(0,\!1)= 1,\!1103418[/math] .

Сопоставляя точное решения [math]y(x_1)[/math] в точке [math]x_1[/math] с результатами [math]\widehat_1[/math] , полученными по схемам второго и четвертого порядка, можно увидеть, что для схемы четвертого порядка результат [math]\widehat_1[/math] совпадает с точным решением в семи цифрах, а для схемы второго порядка — в трех знаках.

В заключение данного раздела отметим преимущества и недостатки одношаговых и многошаговых схем:

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

2. В многошаговых схемах, напротив, учитываются значения искомой функции, рассчитанные на предыдущих шагах интегрирования, что повышает порядок точности этих схем без использования дополнительных точек на отрезке [math][x_,x_][/math] . Однако эти схемы требуют усложнения расчетных алгоритмов при реализации переменного шага.

Замечание. В явных одношаговых методах обеспечить заданную точность на каждом очередном шаге из текущей точки [math](x_, \widehat_)[/math] можно по следующей методике.

1. Задать начальную величину шага [math]h^<(0)>= x_^<(0)>-x_[/math] и рассчитать значение [math]\widehat_^<(0)>[/math] в точке [math]x_^<(0)>[/math] .

2. Вычислить величину [math]\widehat_^<(1)>[/math] в точке [math]x_^<(0)>[/math] с половинным значением шага [math]h^<(1)>= \frac<1><2>h^<(0)>[/math] (до достижения точки [math]x_^<(0)>[/math] выполняются два шага).

3. Если [math]\bigl|\widehat_^<(1)>-\widehat_^<(0)>\bigr| , то перейти к нахождению приближенного решения в точке [math]x_^<(0)>= x_^<(0)>+ h^<(0)>[/math] . Иначе осуществить еще одно дробление шага, т.е. вычислить [math]h^<(2)>= \frac<1><2>h^<(1)>[/math] , и найти [math]\widehat_^<(2)>[/math] в точке [math]x_^<(1)>= x_+ h^<(1)>[/math] .

4. Сравнить решения, полученные с шагом [math]h^<(1)>[/math] и [math]h^<(2)>[/math] в точке [math]x_^<(1)>[/math] . Если желаемая точность достигнута, расчеты продолжить с шагом [math]h^<(1)>[/math] , а если нет, то процесс дробления шага продолжить до достижения заданной точности.

Принцип аналогий

Данный метод следует из обобщения интегрально-интерполяционного принципа. Он основан на применении интерполяционных интегрально-дифференциальных сплайн-функций, преобразования которых — параметрические соотношения, аппроксимационные формулы для производных и квадратурные формулы — используются для получения разностных схем решения задачи Коши. Приведем некоторые из этих соотношений, соответствующих задаче аппроксимации некоторой произвольной функции [math]\Phi(x)[/math] многочленами второй степени:

а) параметрическое соотношение, связывающее интерполируемую функцию [math]\Phi(x_)[/math] и ее производную [math]\overline_= \Pho'(x_)\colon[/math]

б) формула аппроксимации первой производной в точке [math]x_[/math] на трехточечном нерегулярном шаблоне [math](x_, x_,x_)\colon[/math]

s,\,p[/math] — целые числа [math]\left(H_^= h_(s+ p\cdot \delta_),

в) квадратурная формула экстраполяционного типа, имеющая четвертый порядок:

Эта формула при [math]h_= h=\text[/math] упрощается:

г) двухинтервальная (трехточечная) обобщенная на нерегулярный шаблон квадратурная формула парабол четвертого порядка на отрезке [math][x_,x_][/math] , выражающая сумму «удельных» значений интегралов через линейную комбинацию значений [math]\Phi_, \Phi_, \Phi_\colon[/math]

д) одноинтервальная (трехточечная) функциональная квадратурная формула четвертого порядка на отрезке [math][x_,x_][/math]

которая при [math]h_= h=\text[/math] преобразуется к более простому виду:

Анализ остаточного слагаемого, получающегося при оценке погрешности квадратурной формулы (6.50) (см. мажоранту в скобках справа от данной формулы), указывает на то, что при сгущающейся вправо сетке (условная аппроксимация интеграла) порядок точности формулы может быть повышен на единицу без увеличения количества точек шаблона. Действительно, при [math]\delta_ легко получить приближенное условие повышения порядка: [math]h_ \leqslant h_^2[/math] . Таким образом, в данном случае порядок повышается за счет значительного уменьшения шага, такого, что [math]h_[/math] не превышает квадрата шага [math]h_[/math] . Этот факт может быть использован при решении задачи Коши на нерегулярном шаблоне;

е) две одноинтервальные (двухточечные) квадратурные формулы функционально-дифференциального типа:

ж) интегрально-функциональное параметрическое соотношение, осуществляющее связь интегралов и производных:

Дискретные явные или неявные схемы для решения задачи (6.9) получаются из приведенных соотношений либо посредством аппроксимации производной, либо путем соответствующей замены одних параметров, относящихся к параметрическому соотношению или к квадратурной формуле, на другие параметры. Этот способ построения схем называется способом аналогий .

Так, заменяя в (6.46) [math]\Delta\Phi_[/math] на [math]\Delta y_[/math] , а [math]\overline_[/math] на [math]f_

(k=i,i-1)[/math] и разрешая относительно [math]\widehat_[/math] , получаем двухшаговую явную схему второго порядка (2Я2А):

которая при [math]h_= h= \text[/math] преобразуется к виду [math]\widehat_= 2\widehat_-\widehat_+ h(f_-f_)[/math] .

Комбинация схем (6.54) и (6.24) дает обобщенную на нерегулярный шаблон схему Адамса—Бэшфорта (2Я2В):

которая при [math]h_= h= \text[/math] преобразуется к виду (6.32): [math]\widehat_= \widehat_+ \frac<2>(3f_-f_)[/math] .

Аппроксимационная формула (6.47) определяет неявную двухшаговую схему второго порядка (2НЯ2):

которая при [math]h_= h= \text[/math] принимает вид

Замечание. Схема (6.57) является одним из методов Гира , получаемых на основе формул дифференцирования назад, связывающих значение производной в точке [math]x_[/math] со значениями функции в предыдущих точках [math]x_,\ldots, x_[/math] (см. рис.6.2). Схема (6.57) является схемой второго порядка, а неявный метод Эйлера — схемой первого порядка. Приведем более точные многошаговые схемы, записанные на регулярных шаблонах:

– схема третьего порядка:

– схема четвертого порядка:

– схема пятого порядка:

– схема шестого порядка:

Известно, что схема (6.58) является A(α)-устойчивой с [math]\alpha\cong 68^<\circ>[/math] .

Все рассматриваемые ниже схемы получаются из квадратурных формул (6.48)-(6.52), в которых вместо разности первообразных [math]F_-F_= I_^[/math] подставляется [math]\Delta \widehat_[/math] , а вместо [math]\Phi_(k=i-2,i-1,i)[/math] и [math]\Phi’_[/math] подставляются [math]f_[/math] и [math]f’_[/math] . При этом порядки точности схем на отрезке [math][a,b][/math] соответствуют порядкам точности вычисления интегралов на всем отрезке интегрирования.

Преобразуя таким образом квадратурную формулу (6.48), получаем явную трехшаговую схему третьего порядка (ЗЯЗ):

При [math]h_= h= \text[/math] она является традиционной, получающейся интегрально-интерполяционным методом (6.33):

Из обобщенной формулы (6.49) следует первая двухшаговая неявная схема третьего порядка (2НЯЗА), которая при [math]h_= h= \text[/math] преобразуется к виду (6.29)

Одноинтервальная (трехточечная) квадратурная формула (6.50) определяет вторую двухшаговую неявную схему третьего порядка (2НЯЗБ)

Эта схема в соответствии с вышеприведенной оценкой точности интеграла при [math]\delta_ \leqslant h_[/math] имеет не третий, а четвертый порядок точности. При [math]h=\text[/math] эта схема преобразуется к виду (6.36):

И наконец, последние приводимые здесь схемы следуют из квадратурных формул (6.51), (6.52) (неявные одношаговые схемы) и параметрического соотношения (6.53) (неявная двухшаговая схема):

где [math]f’=f_+f_\cdot f[/math] . Эти схемы обозначаются 1НЯЗА, 1НЯЗБ, 2НЯЗД соответственно.

При [math]h_= h= \text[/math] последняя схема преобразуется к виду

Для применения разностных схем, полученных с помощью четырех различных подходов, можно использовать общую методику.

Методика решения задачи Коши методами сеток

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

2. Задать сетку [math]\Omega_[/math] или выбрать ее параметр ( [math]h[/math] , если сетка равномерная; [math]h_,

i=\overline<0,n-1>[/math] , если сетка неравномерная). Если выбранный метод является ограниченно устойчивым, то для задания параметра сетки оценить величину критического шага.

3. Если выбранный метод многошаговый, найти требуемое число «разгонных» точек.

4. Найти приближенное решение задачи Коши во всех точках сетки [math]\Omega_[/math] . Если выбран явный метод, каждая следующая точка [math]\widehat_[/math] рассчитывается по явной формуле, а если неявный, то с помощью решения в общем случае нелинейного алгебраического уравнения (6.11) одним из численных методов.

5. Оценить достигнутую точность. В случае невыполнения требований повторить расчеты, меняя значение шага или используя более точный метод.