Формула прямоугольников

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

,

и на отрезке интеграл заменяется выражением

. (7.5)

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

Рис. 7.1. Схема численного интегрирования методом прямоугольников

Воспользуемся для представления функции f(x) вблизи точки формулой Тейлора

.

Определим погрешность вычисления интеграла на отрезке :

.

Полученное выражение позволяет оценить погрешность:

(7.6)

Здесь обозначено: .

Формула (7.6) устанавливает, что при ограниченности второй производной на рассматриваемом отрезке погрешность формулы прямоугольников имеет третий порядок. Для всего отрезка интегрирования [a, b] получаем

, (7.7)

где .

Иными словами, для всего интервала [a, b] погрешность рассмотренного способа интегрирования имеет второй порядок.

Аналогичным способом можно проверить также часто применяемую на практике формулу интегрирования

,

геометрический смысл которой пояснен на рисунке 7.2.

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

,

.

xk-1 xk x

Рис. 7.2. Схема численного интегрирования методом прямоугольников

Для всего интервала [a, b] погрешность интегрирования составляет

.

На рис. 7.3 приведены графики, отражающие сходимость процесса приближенного вычисления определенного интеграла с помощью формул метода прямоугольников с центральной точкой (формула (7.5)), “левой” точкой (рис. 7.2) и “правой” точкой .

Рис. 7.3. Значения интеграла , вычисленные точно (-------) и по формулам метода прямоугольников с центральной (- o -), “левой” (- D -) и “правой” (- à -) точками на сетках W_n

Формула трапеций

Заменим функцию f(x) на отрезке линейным приближением

.

Это означает, что в разложении (7.2) удерживаются две функции .

Тогда весовые коэффициенты

Отсюда вытекает формула метода трапеций (рис. 7.4):

. (7.8)

Воспользуемся формулами Тейлора

,

.

Рис. 7.4. Схема численного интегрирования методом трапеций

Оценим погрешность представления (7.8) на отрезке :

В силу того, что

,

,

вычисляемая погрешность

.

Получим оценку погрешности:

.

Погрешность для всего отрезка интегрирования [a, b] имеет второй порядок. (7.9)

Формула Симпсона

Заменим на отрезке функцию f(x) полиномом Лагранжа. В частности, для трех точек полином второй степени имеет вид

Это означает, что в разложении (7.2) оставлены три функции:

,

,

.

Для определения коэффициентов вычислим интегралы:

.

Формула приближенного интегрирования Симпсона (рис. 7.4):

. (7.10)

f(x)

x_k-1 x_k-1/2 x_k x

Рис. 7.6. Схема численного интегрирования методом Симпсона

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

,

,

Здесь принято, что .

Подсчитаем интеграл в левой части формулы (7.10):

.

Подсчитаем выражение в правой части (7.10):

.

Определим величину погрешности формулы Симпсона:

.

Модуль погрешности на отрезке :

.

Для всего отрезка [a, b] интегрирования погрешность

(7.11)

имеет четвертый порядок.

В последних выражениях использованы обозначения

, .

Метод Эйлера

Пусть для отрезка [a, b], на котором ищется решение дифференциального уравнения (2.1), построена сеточная область с постоянным шагом h. Для построения решения уравнения (2.1) воспользуемся разложением искомой функции y(x) в ряд Тэйлора вблизи произвольной точки :

Учитывая, что согласно уравнению (2.1) , это разложение решения можно записать в виде

С помощью полученного выражения построим вычислительный процесс

. (2.6)

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

Пример 2.3. Решим методом Эйлера уравнение .

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

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

,

,

,

,...,

.

Для оценки точности получаемого численного решения производятся вычисления с различными шагами интегрирования h (табл. 2.1). Точное решение , например, для x=10 дает значение y(10)=0,453999297610^-4.

Таблица 2.1.

Результаты численного решения y_n методом Эйлера обыкновенного дифференциального уравнения с начальным условием y(0)=1.

Верные три значащие цифры получены для шага h = 0,0001, при этом требуется выполнить 100000 “шагов” по аргументу x, то есть выполнить вычислений в соответствии с формулой (2.6).

Для оценки точности решения, как и ранее, представим погрешность в виде

,

определяющем отклонение численного решения по схеме Эйлера от точного.

Соотношение (2.6) запишем в виде разностного аналога дифференциального уравнения (2.1):

. (2.7)

Подставляя в выражение (2.7) и , получаем

.

Преобразуем последнюю формулу, добавляя в правую часть и одновременно вычитая ,

,

(2.8)

Выражение в квадратных скобках представляет собой разностный аналог (2.7), в котором приближенные значения y_k искомой функции заменены на точные значения y(x_k). Это позволяет проверить, как точно разностный аналог (2.7) аппроксимирует исходное дифференциальное уравнение (2.1).

Величина

(2.9)

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

в формулу (2.9):

Учтем, что согласно исходному уравнению . Тогда

.

При условии, что значения ограничены на отрезке [a, b], погрешность аппроксимации оказывается величиной, пропорциональной первому порядку шага интегрирования.

В общем случае порядок аппроксимации равен p, если имеет место соотношение вида

.

Вновь обратимся к формуле (2.8). Второе слагаемое согласно теореме Лагранжа о среднем можно представить в виде:

.

Тогда выражение (2.8) преобразуется следующим образом:

,

.

Для того, чтобы при выполнении вычислений погрешность не возрастала, потребуем выполнения условия:

,

,

. (2.10)

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

,

.

Учитывая условие устойчивости решения, , получаем:

.

Теперь из неравенства (2.10) можно получить:

.

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

Для построения разностной схемы интегрирования воспользуемся разложением решения задачи (2.1) в ряд Тэйлора:

Заменим вторую производную в этом разложении выражением

,

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

Для исходного уравнения (2.1) построим вычислительную схему:

,

которую преобразуем к виду

Введем обозначения:

,

которые позволяют записать предыдущее выражение в форме

.

Очевидно, что все введенные коэффициенты зависят от величины Dx и могут быть определены через a, который в этом случае играет роль параметра,

.

Окончательно схема Рунге-Кутты принимает вид

. (2.11)

Та же схема в форме разностного аналога уравнения (2.1):

При a= 0 получаем как частный случай уже известную схему Эйлера

.

При a= 1 выражение (2.11) записывается в форме

.

В этом случае проведение расчетов на очередном шаге интегрирования можно рассматривать как последовательность следующих операций (рис. 2.2):

1. Вычисляется выражение

,

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

2. Для той же промежуточной точки определяется приближенное значение производной

.

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

.

Геометрические построения показывают, что получаемое в такой последовательности решение лежит “ближе” к истинному, чем вычисляемое по схеме Эйлера, то есть следует ожидать более высокой точности решения, получаемого методом Рунге-Кутты.

y

y_k+1

y_k y_k+1/2

x_k x_k+1/2 x_k+1

Рис.2.2. Схема интегрирования обыкновенного дифференциального уравнения методом Рунге-Кутты с параметром a = 1.

Теперь рассмотрим схему, получающуюся из выражения (2.11) при (геометрическая интерпретация результата приведена на рис. 2.3):

.

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

.

2. Для этой же точки вычисляется приближенное значение производной

.

3. Находится среднее значение двух производных, определенных на концах отрезка

.

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

.

Геометрически понятно, что получаемый указанным способом результат также должен быть “ближе” к истинному решению, чем получаемый по схеме Эйлера.

y

y_k+1

y_k

x_k x_k+1

Рис.2.3. Схема интегрирования обыкновенного дифференциального уравнения методом Рунге-Кутта с параметром a= 0,5

Не нашли, что искали? Воспользуйтесь поиском: