Методы численного интегрирования уравнений второго порядка

Моделирование динамических систем: численные методы решения ОДУ

Очень кратко рассмотрев основы механики в предыдущей статье, перейдем к практике, ибо даже той краткой теории что была рассмотрена хватит с головой.

Камень бросают вертикально, без начальной скорости с высоты h = 100 м. Пренебрегая сопротивлением воздуха, определить закон движения камня, как функцию высоты камня над поверхностью Земли от времени. Ускорение свободного падения принять равным 10 м/с 2

1. Формализация задачи

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

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

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

где g — ускорение свободного падения у поверхности Земли. Теперь пришло время составить уравнения движения камня. Помните эти уравнения

Левая их часть нас пока не интересует, а вот в правой стоят суммы проекций сил, приложенных к точке на оси координат. Пусть оси x и y располагаются на поверхности, а ось z направлена вверх перпендикулярно ей. Сила одна единственная, её проекции на оси x и y равны нулю, а на ось z проекция отрицательна, так как сила направлена против направления оси, то есть

Масса камня, очевидно не равна нулю, значит можно спокойно поделить обе части получившихся уравнений на неё

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

То что у нас получилось не много не мало — математическая модель процесса происходящего в задаче. Пафосно, да?

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

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

А вот как решить это численно? И что это вообще такое — «численно»?

2. Численное интегрирование дифференциального уравнения первого порядка

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

Во-первых, это уравнение относится к такому типу, что допускает понижение порядка. Правая часть не зависит от неизвестной функции (там нет z), поэтому вспоминаем, что

проекция ускорения на ось z равна первой производной проекции скорости на ту же оcь z. Ну классно, тогда

вот вам и уравнение первого порядка. Не всегда этот номер проходит (не буду я сейчас про форму Коши!), но в данном случае всё в порядке. Будем искать не координату а скорость точки. Что дальше-то? А дальше

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

Что получается? А вот что

Мы получили приращение скорости. Отрицательное приращение. Как это так, камень, падая вниз будет разгонятся же! Да, будет. Его скорость, вектор его скорости, будет направлен вниз. А значит проекция этого вектора на ось z будет отрицательной. Всё правильно, мы получаем растущую по абсолютному значению проекцию вектора, направленного вниз. Мы знаем, начальное значение скорости — ноль, а значит

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

а ещё через 0.1 секунды

и ещё через 0.1 секунды

Хм, так мы можем продолжать довольно долго, но ограничимся промежутком времени в одну секунду

То есть, воспользовавшись формулой

мы получили зависимость скорости точки от времени. А всего-то нужно взять значение скорости в текущий момент времени, и добавить к нему приращение, которое скорость получит в новый момент времени, отстоящий от текущего на секунд. Приращение времени называется здесь шагом интегрирования. А приращение мы вычисляем как значение производной искомой функции в текущий момент времени умноженное на шаг. Просто? Да просто конечно. И та формула, которую я написал, имеет название название — явный метод Эйлера для численного решения дифференциальных уравнений. Это так называемая рекуррентная формула, когда новое значение вычисляемой величины зависит от её предыдущего значения.

А что же с высотой точки над Землей? Да всё аналогично, смотрите

ведь проекция скорости есть производная от соответствующей координаты. Применим формулу метода Эйлера для этого уравнения, ведь скорость мы уже знаем

и по этой формуле добавим в нашу таблицу ещё одну колонку

Хм, ну, во-первых, заметно, что высота меняется у нас уже неравномерно, так как скорость со временем меняется. Теперь наша производная сама зависит от времени. Но уже на первом шаге, мы замечаем неладное — скорость уже есть, а вот высота по прежнему 100 метров. Как так?

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

Да, наш камень как будто зависает в воздухе. Численное решение отстает от аналитического, и чем дальше мы считаем, тем выше погрешность счета. Погрешность накапливается, так как на каждом шаге мы берем всё более и более грубое приближение. Что делать?

Во первых, можно уменьшить шаг. Скажем в 10 раз, пусть секунды

Уже лучше, погрешность в конце счета не превышает 0,05 метров, и это в 10 раз меньше предыдущего значения. Можно предположить, что уменьшив шаг ещё в 10 раз мы получим ещё более точное решение. Я схитрил, выводя значения только для 10 точек с шагом 0.1, на самом деле, чтобы получить такую таблицу нужны уже 100 итерации а не 10. При шаге 0.001 потребуется уже тысяча итераций, а результат будет таким

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

Метод Эйлера самый простой из известных методов интегрирования дифференциальных уравнений. Из нашего простого примера видно, что погрешность метода прямо пропорциональна шагу интегрирования, и это действительно так. Такие методы называются методами 1-го порядка точности.

Точность расчетов даже на шаге 0.1 можно улучшить, если мы применим, скажем метод Рунге-Кутты 4-го порядка точности. Но это отдельная история.

Заключение

Задумайтесь… Мы рассмотрели очень простой пример. Мы даже не применяли компьютер, но уже понимаем принцип, по которому работают те самые мощные в мире суперкомпьютеры, что моделируют ранние этапы жизни Вселенной. Конечно, там всё устроено гораздо сложнее, но принцип лежит этот же самый.

Представьте себе, какой мощный инструмент вы получаете в свои руки. Эта последняя статья, где мы не будем применять компьютер. Я обещал Octave. В следующий раз будет именно он.

Методы численного интегрирования уравнений второго порядка

Системой дифференциальных уравнений называется система вида

где x — независимый аргумент,

y_i — зависимая функция, ,

Функции y_i(x), при подстановке которой система уравнений обращается в тождество, называется решением системой дифференциальных уравнений.

Численные методы решения систем дифференциальных уравнений.

Модифицированный метод Эйлера.

Метод Рунге-Кутта четвертого порядка.

Дифференциальным уравнением второго порядка называется уравнение вида

Функция y(x), при подстановке которой уравнение обращается в тождество, называется решением дифференциального уравнения.

Численно ищется частное решение уравнения (2), которое удовлетворяет заданным начальным условиям, то есть решается задача Коши.

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

.

(3)

Функция f₂(x, y₁, y) в систему (3) введена формально для того, чтобы методы, которые будут показаны ниже, могли быть использованы для решения произвольной системы дифференциальных уравнений первого порядка. Рассмотрим несколько численных методов решения системы (3). Расчетные зависимости для i+1 шага интегрирования имеют следующий вид. Для решения системы из n уравнений расчетные формулы приведены выше. Для решения системы из двух уравнений расчетные формулы удобно записать без двойных индексов в следующем виде:

Метод Рунге-Кутта четвертого порядка.

где h — шаг интегрирования. Начальные условия при численном интегрировании учитываются на нулевом шаге: i=0, x=x₀, y₁=y₁₀, y=y₀.

Контрольное задание по зачетной работе.

Колебания с одной степенью свободы

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

Задание. Численно и аналитически найти:

1. закон движения материальной точки на пружинке х(t),
2. закон изменения силы тока I(t) в колебательном контуре (RLC — цепи) для заданных в табл.1,2 режимов. Построить графики искомых функций.

Свободные незатухающие колебания

Затухающее колебательное движение

Предельное апериодическое движение

Вынужденное колебание без сопротивления

Вынужденное колебание без сопротивления, явление резонанса

Вынужденное колебание с линейным сопротивлением

Вынужденное колебание с линейным сопротивлением, явление резонанса

Методы численного интегрирования

Формулы на основе интерполяционных многочленов

Одним из классических методов вычисления определенных интегралов является применение функциональных квадратурных формул

где [math]q_[/math] — весовые коэффициенты] [math]x_,

j=\overline<1,N>[/math] — некоторые точки отрезка [math][a,b][/math] ; [math]N[/math] — число точек <узлов квадратурной формулы).

Квадратурная формула называется тонной для многочленов степени [math]m[/math] , если при замене функции [math]f(x)[/math] на произвольный алгебраический многочлен степени не выше [math]m[/math] приближенное равенство (5.42) становится точным. В этом случае говорят, что квадратурная формула обладает m-свойством.

1. При приближенном вычислении интеграла, как правило, отрезок [math][a,b][/math] представляется в виде объединения [math]I[/math] непересекающихся частичных отрезков вида [math][x_,x_][/math] , которым соответствует шаблон [math]H_= (x_,\ldots, x_,\ldots, x_)[/math] , где [math]i[/math] — номер базового узла сетки; [math]r[/math] и [math]s[/math] — количество узлов левее и правее узла с номером [math]i[/math] ; [math]k=r+s+1[/math] — общее число узлов (точек) в шаблоне (рис. 5.1). На каждом частичном отрезке с номером [math]j=1,2,\ldots,l[/math] вычисляется интеграл по соответствующей квадратурной формуле

а затем полученные значения суммируются по всем частичным отрезкам, т.е.

2. Далее в силу использования представления (5.44) проблеме вычисления интеграла (5.43) на частичном отрезке уделяется основное внимание. По заданной сеточной функции или сеточному представлению формульной функции на частичном отрезке строится интерполяционный многочлен некоторой степени. Значение [math]\widehat_^<\,i+s>[/math] определяется величиной интеграла от этого многочлена.

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

А. Двухточечный шаблон. На частичном отрезке [math][x_, x_][/math] , которому соответствует двухточечный шаблон [math]H_<2,i>= (x_,x_)[/math] , где [math]r=0,

s=1[/math] , заменим функцию [math]f(x)[/math] двумя способами:

а) интерполяционным многочленом нулевой степени: [math]L_0>(x)= f\! \left(\frac+ x_><2>\right)[/math] , построенным по значению функции [math]f_<2>>= f(x_<2>>)[/math] в середине частичного отрезка [math]x_<2>>= \frac+ x_><2>[/math] (рис. 5.2,а);

б) интерполяционным многочленом первой степени [math]L_<1>(x)[/math] с узловыми значениями [math]x_,\,x_\colon[/math]

где [math]q[/math] — фаза интерполяции, [math]q=\frac>>,

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

[math]\widehat_>^<\,i+1>= \int\limits_>^> L_<1>(x)\,dx= \int\limits_<0>^\bigl[(1-q)f_+ qf_\bigr]h_\,dq= h_\frac+f_><2>\,,[/math]

где при интегрировании учитывалось, что [math]dx=h_\cdot dq;

q=0[/math] при [math]x=x_[/math] и [math]q=1[/math] при [math]x=x_[/math] .

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

Подчеркнем, что данные формулы справедливы и для нерегулярного шаблона, хотя последующее их суммирование по всем частичным отрезкам [math][x_, x_][/math] традиционно выполняется при [math]h=\text[/math] .

Б. Трехточечный шаблон. Пусть отрезок [math][a,b][/math] разбит на четное количество одинаковых частичных отрезков, т.е. [math]n=2k[/math] , где [math]k[/math] — число пар. Проделав аналогичные выкладки для многочлена [math]L_<2>(x)[/math] второй степени, записываемого при [math]h=\text[/math] на шаблоне [math]H_<3,i>= (x_,x_, x_)[/math] (по одной паре отрезков при [math]r=1,

s=1[/math] ), получим двухинтервальную квадратурную формулу парабол, или формулу Симпсона:

Эта же формула может быть получена из (5.23) путем разрешения ее относительно [math]I_^[/math] и подстановки в нее конечно-разностной аппроксимации второй производной [math]\widehat\,»_= \frac<1>(f_-2f_+ f_)[/math] имеющей второй порядок, или из формулы (4.92) с учетом [math]I_^= I_^+ I_^[/math] и [math]h=\text[/math] .

Дадим геометрические интерпретации квадратурных формул прямоугольников (рис. 5.2,а), трапеций (рис. 5.2,б) и парабол (рис. 5.2,в). Интегралы обычно определяются не на частичных отрезках, а на всем отрезке [math][a,b][/math] , и поэтому путем суммирования левых и правых частей (5.45)–(5.47) получаются так называемые составные квадратурные формулы. Для сеточной функции [math]f(x_),

i=\overline<0,n>[/math] , заданной на регулярном шаблоне при [math]h=\text[/math] , эти формулы имеют вид

[math]\begin\widehat_>^<\,b>&= h\! \left(\frac+f_<1>><2>+ \frac+f_<2>><2>+ \ldots+ \ldots+ \frac+f_><2>\right)=\\ &=\frac <2>\bigl[f_<0>+ 2(f_1+f_2+\ldots+f_)+ f_\bigr]=\\ &= \frac<2>\! \left(f_0+ 2 \sum\limits_^ f_+ f_\right)\!, \end[/math]

Подчеркнем, что в составной квадратурной формуле парабол индекс «k» указывает на число пар отрезков разбиения, которое предполагается четным [math](n=2k)[/math] . Если это условие не выполняется, то интеграл вычисляется для четного количества отрезков и к полученному значению добавляется величина [math]I_^[/math] , рассчитанная с порядком [math]O(h^5)[/math] по формулам, приведенным далее.

Оценка погрешностей квадратурных формул

Проведем оценку погрешностей одноинтервальных квадратурных формул прямоугольников и трапеций. С этой целью составим разность [math]I_^-\widehat_^<\,i+1>[/math] , которую затем преобразуем путем разложения первообразных и функций по формуле Тейлора. Для этого [math]I_^[/math] заменим разностью [math]F_-F_[/math] , а [math]F(x)[/math] и функции, входящие в правые части квадратурной формулы, разложим по формуле Тейлора относительно точки [math]x_[/math] ( [math]F(x)[/math] — первообразная, [math]F_= F(x_)[/math] ). Так, для формулы (5.46) при условии, что [math]f(x)= C_2[a,b][/math] , получим

[math]\beginI_^-\widehat_>^<\,i+1>&= F_-F_-\frac><3>(f_+f_)=\\ &= F_+ h_f_+ \frac^2><2>f’_+ \frac^3><6>f»(\xi_)-F_-\frac><2>\! \left(f_+f_+ h_f’_+ \frac^2><2>f»(\xi_)\right)=\\ &= \frac^3><6>f»(\xi_)-\frac^3><4>f»(\xi_)=-\frac^3><12>f»(\xi_), \end[/math]

где [math]\xi_\in (x_, x_)[/math] — одна и та же точка для функций [math]F(x)[/math] и [math]f(x)[/math] . Знак полученной разности указывает на то, что если вторая производная [math]f»(x)[/math] на частичном отрезке [math][x_,x_][/math] положительна, то формула (5.46) аппроксимирует [math]I_^[/math] с избытком. Из полученного соотношения для модуля разности [math]I_^-\widehat_>^<\,i+1>[/math] вытекает оценка [math]\bigl|I_^-\widehat_>^<\,i+1>\bigr| \leqslant \frac> <12>h_^3[/math] , где [math]M_<2,i>= \max_<[x_,x_]> |f»'(x)|[/math] , которая означает, что одноинтервальная формула трапеций аппроксимирует [math]I_^[/math] с третьим порядком по [math]h[/math] .

Чтобы перейти к оценке погрешности аппроксимации интеграла [math]I_^[/math] по составной формуле трапеций на всем отрезке, соотношения для разности [math]I_^-\widehat_>^<\,i+1>=-\frac^3><12>f»(\xi_)[/math] необходимо просуммировать по всем частичным отрезкам (при этом полагается [math]h=\text[/math] ):

[math]I_^-\widehat_>^<\,b>=-\frac <12>\bigl(f»(\xi_<0>)+ f»(\xi_<1>)+ f»(\xi_<2>)+ \ldots+f»(\xi_)\bigr).[/math]

Применяя к сумме в правой части утверждение В.1, получаем

[math]I_^-\widehat_>^<\,b>=-\frac<12>nf»(\xi),\quad \xi\in[a,b].[/math]

Учитывая, что [math]n=\frac,

M_2= \max_<[a,b]>|f»(x)|[/math] , и переходя к неравенству, находим оценку

[math]\bigl|I_^-\widehat_>^<\,b>\bigr| \leqslant \frac <12>(b-a)h^2.[/math]

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

1. Порядок аппроксимации составной квадратурной формулы трапеций на единицу меньше порядка аппроксимации одноинтервальной формулы и равен двум.

2. Величина остаточного слагаемого зависит от длины отрезка интегрирования [math][a,b][/math] и [math]M_<2>[/math] . Константа в правой части равна [math]\frac><12>(b-a)[/math] .

Проделав аналогичные выкладки для одноинтервальной и составной формул прямоугольников и парабол, получим оценки:

где [math]M_<4>= \max_ <[a,b]>\bigl|f^<(4)>(x)\bigr|[/math] при [math]f(x)\in C_<4>[a,b][/math] . Отсюда вытекает, что составная формула прямоугольников имеет на отрезке [math][a,b][/math] второй порядок аппроксимации и константу [math]\frac <24>(b-a)[/math] , которая в два раза меньше соответствующей константы в оценке погрешности квадратурной формулы трапеций. При 0″>[math]f»(x)>0[/math] на частичном отрезке [math][x_, x_][/math] одноинтервальная формула (5.45) аппроксимирует [math]I_^[/math] с недостатком.

Замечание. Между максимальной степенью многочленов, для которых квадратурная формула является точной, и порядком аппроксимации по отношению к шагу [math]h[/math] имеется связь. Формулы прямоугольников и трапеций точны для многочленов первой степени и обладают вторым порядком аппроксимации. Формула Симпсона является точной для многочленов третьей степени и имеет четвертый порядок аппроксимации.

Вышеприведенные оценки погрешностей аппроксимационных квадратурных формул позволяют для непрерывных функций априорным путем вычислять шаг интегрирования, который обеспечивает заданную точность. Для сеточных функций априорные оценки погрешностей выполняются путем предварительной аппроксимации функции [math]y_= f(x_)[/math] и определения соответствующей величины [math]M_

(k=2,3,4)[/math] по формулам численного дифференцирования. Далее опишем методику вычисления интеграла с априорным нахождением шага [math]h[/math] .

Методика вычисления определенного интеграла с заданной точностью

1. Для правой части формулы оценки погрешностей вычислить константу [math]M_

= \max_ <[a,b]>\bigl|f^<(p)>(x)\bigr|[/math] с этой целью необходимо продифференцировать функцию [math]p[/math] раз и вычислить ее максимальное значение на отрезке [math]a,b[/math] , где [math]p[/math] — порядок аппроксимации квадратурной формулы.

3. По значению [math]h[/math] вычислить [math]n[/math] — количество разбиений отрезка [math][a,b][/math] и сформировать сеточное представление функции [math]y=f(x)[/math] , то есть

4. Полученную сеточную функцию подставить в правую часть соответствующей квадратурной формулы и вычислить искомое значение [math]\widehat_^<\,b>[/math] . При этом значение интеграла в силу справедливости оценки удовлетворяет заданной точности [math]\varepsilon[/math] .

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

Пример 5.8. Вычислить на основе квадратурных формул трапеций интеграл [math]\textstyle^ <2>\dfrac>[/math] точностью [math]\varepsilon= 0,\!01[/math] . Оценить фактическую погрешность вычисленного значения.

1. Поскольку [math]p=2[/math] , вычислим [math]M_2= \max_ <[1;2]>\bigl|f»(x)\bigr|[/math] . Дифференцируя [math]f(x)=\frac<1>[/math] дважды, получаем [math]f»(x)= \frac<2>[/math] . Данная функция монотонно убывающая, поэтому она достигает максимума в точке [math]x=1[/math] . Таким образом, [math]M_2=2[/math] .

2. Найдем шаг интегрирования, который обеспечивает заданную точность, из условия [math]\frac <12>(b-a)h^2 \leqslant \varepsilon[/math] . Таким образом, [math]h \leqslant \sqrt<\frac<12\cdot \varepsilon><(b-a)M_2>>= \sqrt<\frac<12\cdot 0,\!1><1\cdot 2>>= 0,\!245[/math] . Примем [math]h=0,\!2[/math] .

3. По значению [math]h=0,\!2[/math] вычислим [math]n=\frac= \frac<1><0,\!2>=5[/math] и сформируем сеточное представление функции, т.е.

4. Вычислим интеграл по формуле (5.49) при [math]n=5\colon[/math]

[math]\begin\widehat_<1,\text>^<\,2>&= \frac <2>\bigl[f_0+ 2(f_1 +f_2+ f_3+ f_4)+ f_5\bigr]=\\ &=0,\!1\cdot \bigl[1+ 2\cdot \bigl(0,\!83(3)+ 0,\!7142+ 0,\!625+ 0,\!5(5)\bigr)+ 0,\!5\bigr]= 0,\!695. \end[/math]

Точное значение интеграла [math]I_1^2= \ln2= 0,\!693147[/math] . Оценим фактическую абсолютную погрешность полученного (приближенного) значения:

Из оценки следует, что фактическая погрешность не превышает 0,01, то есть [math]\delta_<\text> .

Пример 5.9. Вычислить интеграл [math]\textstyle^ <2>\dfrac>[/math] на основе квадратурной формулы парабол с точностью [math]\varepsilon=0,\!01[/math] . Результат сопоставить с полученным в примере 5.8.

1,2. Аналогично примеру 5.8 найдем шаг интегрирования из неравенства

3. По значению [math]h=0,\!5[/math] вычислим [math]n= \frac= \frac<2-1><0,\!5>= 2=2k[/math] . Следовательно, [math]k=1[/math] . Сформируем сеточное представление функции:

Вычислим интеграл по формуле (5.50): [math]\widehat_<1,\text>^<\,2>= \frac<3>(f_0+ 4f_1+ f_2)= 0,\!6944[/math] . Фактическая абсолютная погрешность [math]\delta_<\text>= |0,\!69444-0,\!693147|= 0,\!001293[/math] . Таким образом, фактическая погрешность не превышает [math]\varepsilon=0,\!01[/math] .

Шаг интегрирования по формуле парабол получился более чем в 2 раза крупнее шага [math]h[/math] , соответствующего квадратурной формуле трапеций. Это является следствием того, что (5.50) имеет четвертый порядок аппроксимации, а (5.49) — второй. В связи с повышенной точностью квадратурной формулы парабол она нашла широкое применение в вычислительной практике.

1. Рассмотренный способ вычисления интегралов, когда с использованием оценок и точности [math]\varepsilon[/math] предварительно вычисляется шаг интегрирования [math]h[/math] , является способом с априорным определением шага [math]h[/math] . В вычислительной практике большое распространение получил и другой апостериорный способ вычисления [math]h[/math] , основанный на многократном дроблении отрезка [math][a,b][/math] . В этом случае интеграл вычисляется с помощью итерационного алгоритма методом Рунге.

2. Существуют и другие широко используемые квадратурные формулы.

Одноинтервальная формула прямоугольников (немодифицированная) на двухточечном шаблоне [math]H_<2,i>= (x_,x_)[/math] , где [math]r=0,

s=1[/math] (рис. 5.3,а): [math]\widehat_^<\,i+1>= h\cdot f_

Четырехинтервальная формула Боде на пятиточечном шаблоне [math]H_<5,i>= (x_, x_, x_,x_, x_)[/math] , где [math]r=2,

Шестиинтервальная формула Уэддля на семиточечном шаблоне [math]H_<7,i>= (x_, x_, x_, x_,x_, x_,x_)[/math] , где [math]r=3,

Формулы Ньютона-Котеса (приведем два частных случая):

– трехинтервальная формула на четырехточечном шаблоне [math]H_<4,i>= (x_, x_, x_, x_)[/math] , где [math]r=2,

s=1[/math] (формула «трех восьмых», рис. 5.3,б):

– шестиинтервальная формула на семиточечном шаблоне [math]H_<7,i>= (x_, x_, x_, x_,x_, x_,x_)[/math] , где [math]r=3,

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

Пример 5.10. Вычислить интегралы по формулам прямоугольников, трапеций, парабол с шагом [math]h=1[/math] . Найти оценки погрешностей.

Точные значения интегралов:

Для формул прямоугольников и трапеций порядок аппроксимации [math]p=2[/math] , а для формулы парабол [math]p=4[/math] . В поставленной задаче [math]a=0,

b=2[/math] . Сначала получим оценки погрешностей априорным способом.

f(x)=x^3[/math] ;
[math]M_4=24[/math] для функции [math]f(x)=x^4[/math] .

Согласно (5.51)–(5.53), справедливы оценки:

[math]|\varepsilon_<\text>| \leqslant \frac <24>(b-a)h^2;\qquad |\varepsilon_<\text>| \leqslant \frac <12>(b-a)h^2;\qquad |\varepsilon_<\text>| \leqslant \frac <180>(b-a)h^4.[/math] Оценки погрешностей формулы трапеций:
[math]|\varepsilon_<\text>|=0[/math] для [math]f(x)=x[/math] ; [math]|\varepsilon_<\text>| \leqslant \frac<2><12>\cdot2\cdot1^2= 0,\!3(3)[/math] для [math]f(x)=x^2[/math] ;
[math]|\varepsilon_<\text>| \leqslant \frac<12><12>\cdot2\cdot1^2=2[/math] для [math]f(x)=x^3[/math] ; [math]|\varepsilon_<\text>| \leqslant \frac<48><12>\cdot2\cdot1^2=8[/math] для [math]f(x)=x^4[/math] .

f(x)=x^3[/math] ;
[math]|\varepsilon_<\text>| \leqslant \frac<24><180>\cdot2\cdot1^2= 0,\!26(6)[/math] для [math]f(x)=x^4[/math] .

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

Теперь рассчитаем значения интегралов по соответствующим квадратурным формулам. При [math]h=1[/math] сеточное представление функций имеет вид

(0); &\quad & \widehat_2= 1\cdot\! \left(\frac<1><4>+ \frac<9><4>\right)=2,\!5

(0,\!5); &\quad & \widehat_4= 1\cdot\! \left(\frac<1><16>+ \frac<81><16>\right)=5,\!125

Здесь в скобках указана величина фактической ошибки.

По формуле трапеций (5.49) находим [math]\widehat_<\text>= \frac <2>[f_0+ 2f_1+ f_2][/math] , в частности:

(0); &\quad & \widehat_2= \frac<1><2>\cdot\! \left(0+2+4\right)=3

(0,\!3(3));\\ & \widehat_3= \frac<1><2>\cdot\! \left(0+2+8\right)=5

(1); &\quad & \widehat_4= \frac<1><2>\cdot\! \left(0+2+16\right)=9

По формуле парабол (5.50), учитывая, что [math]n=2k=2[/math] и, следовательно, [math]k=1[/math] , получаем [math]\widehat_<\text>= \frac <3>[f_0+ 4f_1+ f_2][/math] , в частности:

(0); &\quad & \widehat_2= \frac<1><3>\cdot\! \left(0+4+4\right)=\frac <8>

(0);\\ & \widehat_3= \frac<1><3>\cdot\! \left(0+4+8\right)=5

(0); &\quad & \widehat_4= \frac<1><3>\cdot\! \left(0+4+16\right)=\frac <20>

Очевидно, полученные фактические погрешности соответствуют вычисленным ранее оценкам.

Пример 5.11. Вычислить интеграл [math]\textstyle^<\pi\!\not<\phantom<|>>\,\,2> \sin\,dx>[/math] методами прямоугольников, трапеций, парабол, Боде, Уэддля, Ньютона-Котеса на семиточечном шаблоне, применяя разбиение на 1,2,4,8 частичных отрезка.

Точное значение интеграла: [math]I= \int\limits_<0>^<\pi\!\not<\phantom<|>>\,\,2> \sin\,dx=1[/math] . Разобьем отрезок [math]\left[0;\frac<\pi><2>\right][/math] на [math]l=1,2,4,8[/math] частичных отрезка. Применим методы прямоугольников, трапеций, парабол, используя составные формулы (5.48)–(5.50). Остальные значения интеграла подсчитаем по формуле (5.44) при [math]l=1,2,4,8[/math] , используя на каждом частичном отрезке соответствующую квадратурную формулу Боде, Уэддля, Ньютона-Котеса на семиточечном шаблоне. Результаты расчетов приведем в табл. 5.4.

5.4>>\\\hline \text & l=1 & l=2 & l=4 & l=8 \\\hline \text& 1,\!110720338230& 1,\!026171820190& 1,\!006454224265& 1,\!001607874019 \\\hline \text& 0,\!785398006439& 0,\!948059172335& 0,\!987115496263& 0,\!996784860265 \\\hline \text& 1,\!002279560960& 1,\!000134270907& 1,\!000007981598& 1,\!000000202767 \\\hline \text& 0,\!999991251569& 0,\!999999562310& 0,\!999999684178& 0,\!999999686054 \\\hline \text& 0,\!999999293425& 0,\!999999680058& 0,\!999999685980& 0,\!999999686082 \\\hline \begin&\text\\[-4pt] &\text\end& 0,\!999999711921& 0,\!999999686177& 0,\!999999686084& 0,\!999999686083 \\\hline \end[/math]

Из данного примера видно, что для повышения точности возможны два альтернативных подхода:

а) увеличивать количество частичных отрезков, и на каждом из них использовать более простую квадратурную формулу;

б) применять более точные формулы на пяти и семиточечном шаблонах при малом количестве частичных отрезков.

Формулы, полученные на основе сплайнов

Явные формулы. Пусть в общем случае на неравномерной сетке [math]\Omega_n[/math] с некоторой точностью задана сеточная функция [math]y_= f(x_),

i=\overline<0,n>[/math] . Предположим, что точность этого задания в каждом из узлов не ниже [math]O(h_^3)[/math] . В частном случае [math]f(x_)[/math] является сеточным представлением формульной функции, и тогда сетка может быть выбрана исходя из характера [math]f(x)[/math] . Например, в зонах больших изменений [math]f(x)[/math] или ее производных, которые могут быть разрывными, шаг сетки желательно уменьшать. Закон мельчения сетки можно выбирать, например, по формулам арифметической или геометрической прогрессий. В последнем случае в рассмотрение вводится параметр нерегулярности сетки [math]\delta_= \frac>>[/math] , совпадающий со знаменателем прогрессии. Поэтому данный параметр будем включать здесь в правые части некоторых квадратурных формул, и с его использованием проводить оценку погрешности аппроксимации.

Из параболических интегрально-дифференциальных сплайнов путем анализа их параметрических соотношений на трехточечном шаблоне [math]H_<3,i>= (x_,x_,x_)[/math] получаются следующие лево- и правосторонние безусловные аппроксимации интегралов [math]I_^,\, I_^[/math] четвертого порядка (одноинтервальные функциональные аппроксимации):

Предполагая, что [math]f(x)\in C_3[a,b][/math] , и используя методику оценки погрешности, для данных аппроксимаций нетрудно получить оценки, где [math]M_<3,i>= \max_<[x_,x_]> \bigl|f»'(x)\bigr|\colon[/math]

При [math]h=\text[/math] квадратурные формулы (5.54), (5.55) упрощаются (условная аппроксимация):

Их оценки при [math]\delta_=1[/math] становятся одинаковыми (с мажорантой [math]\frac<24>M_<3,i>[/math] ).

Таким образом, из оценок для [math]\widehat_^<\,i+1>,\, \widehat_^<\,i>[/math] следует, что при безусловной аппроксимации обе квадратурные формулы имеют четвертый порядок по [math]h[/math] . Однако для квадратурной формулы (5.55) этот порядок можно повысить, если применить условную аппроксимацию, в которой положить [math]\delta_ \leqslant h_[/math] . В этом случае из формулы [math]\delta_= \frac>> получаем неравенство [math]h_ устанавливающее очень малый шаг [math]h_[/math] , такой, что [math]x_[/math] отстоит от [math]x_[/math] на величину второго порядка по [math]h_[/math] . Этот тип шаблона можно назвать «сильно сжатым», так как он как бы устремляет трехточечный шаблон к двухточечному, и в связи с этим сама величина интеграла [math]\widehat_^<\,i+1>[/math] является малой.

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

Предполагая, что разбиение отрезка [math][a,b][/math] четное, т.е. [math]n=2k[/math] где [math]k[/math] — количество пар разбиения, просуммируем (5.56) по [math]k[/math] . Получаем обобщенную составную квадратурную формулу парабол

Известно, что для квадратурных формул важно, чтобы их коэффициенты были положительны. Это свойство, называемое свойством положительности коэффициентов, обеспечивает минимум погрешности вычисления интеграла по данной квадратурной формуле, а также сходимость вычислительного процесса на последовательно сгущающихся сетках. В квадратурной формуле (5.56) свойство положительности будет выполнено, если 0,»>[math]2-\delta_>0,[/math] 0″>[math]2-\frac<1><\delta_>>0[/math] . Отсюда вытекают пределы изменения параметра нерегулярности [math]\delta_\colon\, 0,\!5 .

Предполагая, что [math]f(x)\in C_3[a,b][/math] , для (5.56) получаем остаточное слагаемое и оценку погрешности:

Переходя в (5.58) к мажоранте для оценки интеграла на всем отрезке [math][a,b][/math] , имеем

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

В общем случае квадратурная формула (5.56) имеет четвертый порядок аппроксимации (безусловная аппроксимация). Возможны также условные аппроксимации двух типов:

а) при [math]\delta_=1[/math] реализуется регулярный шаблон и тогда порядок аппроксимации повышается на единицу, а (5.56) переходит в классическую формулу парабол (Симпсона);

б) при [math]|\delta_-1|\leqslant h_[/math] , реализуется нерегулярный шаблон, и тогда порядок аппроксимации может быть повышен на единицу за счет небольшого отклонения [math]\delta_[/math] от единицы: [math]h_-h_^2 \leqslant h_\leqslant h_+ h_^2[/math] . В этом случае шаблон является квазирегулярным. Такие разбиения с локальными сгущениями узлов сетки могут формироваться при численном интегрировании с помощью квадратурной формулы (5.56) сильно изменяющихся функций или функций, имеющих в некоторых окрестностях разрывы производных.

Неявные методы. Из анализа параметрических соотношений (4.101) и (4.105) двух типов параболических интегрально-дифференциальных сплайнов, изложенных в разд. 4.5.5, для внутренних узлов получаются системы линейных алгебраических уравнений относительно интегралов (при [math]i=\overline<1,n-2>[/math] ):

которые при [math]h=\text[/math] преобразуются к виду

Для вычисления интегралов на всех отрезках эти системы необходимо дополнить двумя граничными условиями. Это можно сделать двумя способами:

1. Системы дополняют значениями интегралов [math]I_0^1[/math] и [math]I_^[/math] на конечных отрезках, которые вычисляются по односторонним формулам (5.54),(5.55) соответственно.

2. Каждая из систем может быть дополнена двумя связями при [math]i=1[/math] и [math]i=n-1[/math] , которые выражаются из соотношения

1. Первая из приведенных здесь систем связывает искомые интегралы со значениями интегрируемой функции, а вторая — со значениями производных.

2. Описанные методы применимы для сеточных функций с [math]n\geqslant3[/math] .

3. Полученные системы линейных алгебраических уравнений являются трехдиагональными и решаются методом прогонки.

Применение кубических интегрально-дифференциальных сплайнов

Явные формулы. Если кубические многочлены, построенные на [math][x_,x_][/math] с помощью функциональных, дифференциальных и интегральных параметров [math]f_,\, f’_,\, I_^[/math] , использовать для вывода квадратурных формул, то получаются различные типы одноинтервальных функционально-дифференциальных квадратурных формул Эйлера пятого порядка.

Квадратурная формула Эйлера-Маклорена имеет вид

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

Подчеркнем, что производные [math]f’_[/math] для подстановки в (5.60) предварительно должны быть вычислены по вышеприведенным аппроксимационным формулам с порядком не ниже третьего. Сокращать производные в последнем слагаемом составной квадратурной формулы, как это обычно делается в учебниках по численным методам, не следует, поскольку при этом в сущности предполагается, что все производные во внутренних узлах известны точно. Это условие для сеточных функций, от которых рассчитываются интегралы, не должно использоваться, так как производные этих функций вычисляются с некоторой погрешностью, порядок которой при суммировании по всем частичным отрезкам понижается на единицу. Справедливость этого утверждения подтверждается следующим примером. После суммирования и сокращения в последней сумме квадратурной формулы (5.60) производных во всех внутренних узлах при [math]h=\text[/math] формула (5.60) приобретает вид

Если квадратурную формулу (5.61) использовать для вычисления интеграла от функции с нулевыми производными на концах или с их одинаковыми значениями [math]f’_= f’_<0>[/math] , то квадратурная формула четвертого порядка переходит в квадратурную формулу трапеций, имеющую второй порядок. Таким образом, реализуется потеря двух порядков, что свидетельствует о необоснованности проведенного сокращения.

Квадратурную формулу (5.59) целесообразно применить для получения функциональных квадратурных формул путем подстановки в них конкретных аппроксимаций [math]\widehat\,’_,\, \widehat\,’_[/math] . Подставив, например, оператор [math]\widehat\,’_= \frac-f_><2h>[/math] , а также аппроксимации производных в левой и правой крайних точках четырехточечного шаблона, получим одноинтервалъные функциональные квадратурные формулы — центральную и лево- и правостороннюю:

Последние три квадратурные формулы (5.62)–(5.64) могут быть объединены в единую составную квадратурную формулу:

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

Неявный метод. Неявный метод численного интегрирования получается путем анализа параметрических соотношений кубических интегрально-дифференциальных сплайнов. После преобразования соотношения (5.59) путем подстановки в него производных на нерегулярной сетке для внутренних узлов получается система линейных алгебраических уравнений относительно интегралов (при [math]i=1,2,\ldots,n-2[/math] )

На регулярной сетке эта система преобразуется к более простому виду:

Для вычисления интегралов на всех отрезках эти системы необходимо дополнить значениями [math]I_0^1[/math] и [math]I_^[/math] , которые могут быть вычислены по квадратурной формуле (5.59) при [math]h_i=\text[/math] или по формулам (5.63), (5.64) при [math]h=\text[/math] . Полученная в результате система линейных алгебраических уравнений решается методом прогонки.

Система (5.65) удовлетворяет условию преобладания диагональных элементов в следующем диапазоне изменения параметра нерегулярности: [math]\Delta_<\ast>^ <-1>\leqslant \delta_ \leqslant \delta_<\ast>[/math] , где [math]\delta_<\ast>\approx 1,\!722[/math] .

Применение интегрально-дифференциальных сплайнов четвертой степени

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

Одноинтервальные явные квадратурные формулы, представляющие функционально-дифференциальные аппроксимации на трехточечном шаблоне [math]H_<3,i>= (x_,x_,x_)[/math] , имеют вид

k=i,i+1[/math] , соответственно для (5.67),(5.68). На регулярном шаблоне эти формулы упрощаются:

Суммирование формул (5.69), (5.70) приводит к двухинтервальной квадратурной формуле (модифицированной формуле Эйлера):

Порядок аппроксимации в этой формуле повышается на единицу из-за симметричного вида формулы (5.71).

Подставляя в функционально-дифференциальные квадратурные формулы аппроксимации производных на пятиточечном шаблоне с порядком [math]O(h^4)[/math] , получаем функциональные квадратурные формулы шестого порядка:

Подстановка в формулы (5.69), (5.70) лево- и правосторонних аппроксимаций производных в крайних точках шаблона приводит к лево- и правосторонним квадратурным формулам шестого порядка:

Для расчета интегралов [math]\widehat_^<\,i>,\, \widehat_^<\,i+1>[/math] на нерегулярной сетке по формулам (5.67), (5.68) предварительно необходимо рассчитать производные в соответствующих точках с порядком не ниже четвертого. С этой целью необходимо воспользоваться интерполяционным многочленом четвертой степени.

Неявный метод. В заключение данного раздела приведем неявный метод вычислений интегралов на всех частичных отрезках путем решения системы линейных алгебраических уравнений трехдиагонального вида [math](h=\text)\colon[/math]

Данная система должна быть дополнена значениями интегралов [math]I_<0>^<1>,\,I_^[/math] , которые рассчитываются по вышеприведенным квадратурным формулам.

Можно показать, что неявный алгоритм вычисления интеграла на отрезке [math][a,b][/math] имеет оценку: [math]\bigl|I_^-\widehat_^<\,b>\bigr| \leqslant \frac<180>M_<5>[/math] .

Формулы на основе разложения первообразных по формуле Тейлора

Пусть [math]f(x)\in C_[a,b][/math] — формульная функция. Тогда из разложений первообразной [math]F(x)[/math] по формуле Тейлора относительно точек [math]x_[/math] и [math]x_[/math] и выражений для [math]F_= F(x_)[/math] и [math]F_= F(x_)[/math] после преобразований получаются одноинтервальные обобщенные квадратурные формулы Эйлера, которые имеют функционально-дифференциальный тип:

В скобках справа от квадратурных формул (5.75), (5.76) указаны мажоранты, характеризующие оценки их погрешностей. На рис. 5.4 показан характер аппроксимации интеграла [math]\widehat_^<\,i+1>[/math] на отрезке [math][x_,x_][/math] первыми (функциональными) слагаемыми комбинаций (5.75), (5.76). Эти аппроксимации представляют собой площади прямоугольников [math]abcd[/math] (рис. 5.4,а) и [math]ABCD[/math] (рис. 5.4,б), построенных по значениям функций [math]f(x)[/math] в левой точке [math]x_[/math] и в правой точке [math]x_[/math] соответственно.

Последующие слагаемые учитывают более высокие по порядку [math]h_[/math] величины, и они вычисляются по значению [math]h_[/math] и производным от первого до m-го порядка в точках [math]x_[/math] и [math]x_[/math] . В соответствии с этим (5.75), (5.76) можно назвать обобщенными квадратурными формулами прямоугольников.

Составляя линейную комбинацию (5.75), (5.76) (среднеарифметическое [math]\widehat_^<\,i+1>[/math] и [math]\widehat_^<\,i+1>[/math] ), получаем обобщенную квадратурную формулу трапеций:

где [math]t=m+3[/math] при четном [math]m[/math] и [math]m=0;

t=m+2[/math] при нечетном [math]m[/math] . Первое слагаемое данной квадратурной формулы соответствует площади трапеции, основания которой равны [math]f_[/math] и [math]f_[/math] , а высотой является [math]h_[/math] (см. рис. 5.2,б). Из выражения для [math]\widehat_^<\,i+1>[/math] видно, что порядок аппроксимации интеграла этим слагаемым квадратурной формулы равен трем, в то время как для (5.75), (5.76) этот порядок равен двум. Это обусловлено знаками производных в комбинациях слагаемых квадратурной формулы (5.77) (при [math]m=1[/math] во втором слагаемом производные [math]f’_[/math] и [math]f’_[/math] вычитаются и. если это слагаемое остаточное, то оно будет равно нулю, так как в этом случае обе производные записываются в одной и той же точке [math]\xi\in (x_,x_)[/math] ).

Аналогичным путем на трехточечном шаблоне [math]H_<3,i>= (x_,x_, x_)[/math] можно получить двухинтервальную квадратурную формулу Эйлера, где [math]H= \max(h_,h_)\colon[/math]

Замечание. На основе рассмотренных здесь квадратурных формул дифференциального типа можно строить вложенные алгоритмы численного интегрирования с заданной точностью, состоящие в том, что на фиксированном шаге очередное слагаемое суммы вычисляется и добавляется к ранее вычисленному результату до тех пор, пока абсолютная величина этого слагаемого не станет меньше заданной величины 0″>[math]\varepsilon>0[/math] .

Пример 5.12. Применяя вложенный алгоритм, вычислить интеграл [math]\textstyle^<\pi\!\not<\phantom<|>>\,\,2> \cos x\,dx>[/math] о с точностью [math]\varepsilon=0,\!01[/math] на основе обобщенной квадратурной формулы трапеций (5.77). Результат сравнить с его точным значением.

Для функции [math]f(x)[/math] с ограниченными производными квадратурная формула (5.77) содержит слагаемые, которые при [math]h_ уменьшаются по абсолютной величине при переходе от текущего слагаемого к последующему. Поэтому достигнутая точность может контролироваться в процессе вычислений по величине очередного слагаемого (по модулю). Как только это слагаемое по модулю становится меньше заданной точности [math]\varepsilon[/math] , то вычисления можно закончить. Проведем вычисления слагаемых: 1,2,3-го и т.д. для заданной в примере функции [math]f(x)=\cos[/math] . При этом принимается [math]x_=0;

x_= \pi\!\!\not<\phantom<|>>\,2[/math] . В результате имеем

Следующее седьмое слагаемое будет по модулю меньше заданной точности [math]\varepsilon= 0,\!01[/math] , поэтому его и последующие слагаемые можно не учитывать. Сумма полученных слагаемых с учетом их знаков и знаков линейной комбинации (5.77) дает результат:

Отличие полученного значения от точного не превышает 0,3%. Отметим, что для данной функции в качестве отрезка [math][x_,x_][/math] выбран целиком весь отрезок [math][a,b]= \left[0; \frac<\pi><2>\right][/math] . Для достижения заданной точности в примере 5.12 потребовалось вычислить шесть слагаемых. Если отрезок [math][a,b][/math] разбить на несколько частичных отрезков, например на два или три, то на каждом шаге потребовалось бы меньшее количество слагаемых и в результате общее их количество могло бы уменьшиться. В связи с этим в процессе реализации вычислительного алгоритма возникает задача выбора оптимальных шагов, которые для достижения той же точности обеспечивают минимальное общее количество слагаемых и в силу этого минимальное количество операций.