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

ТАГАНРОГСКИЙ ГОСУДАРСТВЕННЫЙ РАДИОТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

Чистяков А.Е., Савицкий О.А., Чистякова Т.А

Кафедра высшей математики

Руководство к лабораторным работам по курсу «Численные методы» Таганрог 2009

УДК

Составители: Чистяков А.Е., Савицкий О.А., Чистякова Т.А.

Руководство к лабораторным работам по курсу «Численные методы». Таганрог: ТТИ ЮФУ, 2009. 88 с.

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

Целью работы является обучение студентов работе с задачами, требующими большого объема вычислительной работы, с использованием универсальных решающих программ типа MathCad.

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

Рецензент: Н.Е. Сапунцов, канд. техн. наук, доцент кафедры высшей математики ТТИ ЮФУ

Ó Таганрогский технологический институт,

южного федерального университета 2009

Введение

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

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

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

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

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

Учебно-практическое пособие по курсу «Численные методы» содержит краткое описание наиболее широко используемых на практике методов решения систем линейных алгебраических уравнений, элементы теории интерполирования, численного интегрирования.

Численное интегрирование

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

Пользуясь малостью h, заменим интеграл (1.3) выражением , где . = x_i – 0, 5h.

Тогда получим формулу

(1.4)

которая называется формулой прямоугольников на частичном отрезке [x_i-1, x_i].

Погрешность формулы (1.4) определяется величиной

которую легко оценить с помощью формулы Тейлора. Действительно, запишем y_i в виде

(1.5)

и воспользуемся разложением

где z_i = z_i(x) Î [x_i-1, x_i]. Тогда из (1.5) получим

Обозначая , оценим y_i следующим образом:

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

(1.6)

т.е. формула имеет погрешность О(h³) при h ® 0.

Заметим, что оценка (1.6) является неулучшаемой, т.е. существует функция f(x), для которой (1.6) выполняется со знаком равенства. Действительно, для f(x) = (x – x_i-1/2)² имеем М_{2, i} = 2,

f(x_i-1/2) = 0 и

Суммируя равенства (1.4) по i от 1 до N, получим составную формулу прямоугольников ( центральных прямоугольников )

(1.7)

Погрешность этой формулы

равна сумме погрешностей по всем частичным отрезкам,

Отсюда, обозначая , получим

(1.8)

т.е. погрешность формулы прямоугольников на всем отрезке есть величина О(h²).

В этом случае говорят, что квадратурная формула имеет второй порядок точности.

Замечание. Можно также использовать формулы прямоугольников при ином расположении узлов, например, такие формулы (формулы левых и правых прямоугольников соответственно):

Однако, из-за нарушения симметрии погрешность таких формул является величиной О(h).

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

На частичном отрезке эта формула имеет вид

(1.9)

и получается путем замены подынтегральной функции f(x) интерполяционным многочленом первой степени, построенным по узлам x_i-1, x_i, т.е. функцией

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

Отсюда получим

и, следовательно,

(1.10)

Оценка (1.10) неулучшаема, так как в ней достигается равенство, например, для f(x) = (x – x_i)².

Составная формула трапеций имеет вид

(1.11)

где f_i = f(x_i), i = 0, 1, …, N, x_i = a + ih, hN = b – a.

Погрешность этой формулы оценивается следующим образом:

Таким образом, формула трапеций имеет, так же как и формула прямоугольников, второй порядок точности, y = О(h²), но ее погрешность оценивается величиной в два раза большей (см. (1.8)).

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

При аппроксимации интеграла (1.3) заменим функцию f(x) параболой, проходящей через точки (x_j, f(x_j)), j = i – 1, i – 0, 5, i, т.е. представим приближенно f(x) в виде

f(x) » L_{2, i}(x), x Î [x_i-1, x_i],

где L_{2, i}(x) – интерполяционный многочлен Лагранжа второй степени,

(1.12)

Проводя интегрирование, получим

Таким образом, приходим к приближенному равенству

(1.13)

которое называется формулой Симпсона или формулой парабол.

На всем отрезке [a, b] формула Симпсона имеет вид

Чтобы не использовать дробных индексов, можно обозначить

x_i = a + 0, 5hi, f_i = f(x_i), i = 0, 1, …, 2N, hN = b - a

и записать формулу Симпсона в виде

(1.14)

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

если, f(x) = a₀ + a₁x + a₂x²+ a₃x³. Это утверждение нетрудно проверить непосредственно, что и предоставляется сделать читателю.

Для оценки погрешности формулы Симпсона построим многочлен третьей степени Н₃(х) такой, что

H₃(x_i-1) = f(x_i-1), H₃(x_i-1/2) = f(x_i-1/2),

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

(1.15)

Представим теперь f(x) в виде

f(x) = H₃(x) + r_i(x), x Î [x_i-1, x_i],

(1.16)

где r_i(x) – погрешность интерполирования многочленом Н₃(х). Интегрируя (1.16) и учитывая (1.15) получим

(1.17)

Имеем

поэтому для погрешности y_i получаем оценку

где .

Вычисляя интеграл, приходим окончательно к оценке

(1.18)

Погрешность составной формулы Симпсона (1.14) оценивается так:

Отсюда видно, что формула Симпсона существенно точнее, чем формулы прямоугольников и трапеций. На частичном отрезке она имеет точность О(h⁵), а на всем отрезке – О(h⁴).

Пример выполнения лабораторной работы №1

Применим методы численного интегрирования для вычисления интеграла .

Задаем число разбиений

Устанавливаем пределы интегрирования

Вычисляем шаг сетки

Вводим подынтегральную функцию

Рассчитываем точное значение интеграла

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

Выводим полученное значение

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

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

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

Рассчитываем значение интеграла и погрешности методом трапеции

Рассчитываем значение интеграла и погрешности методом Симпсона

Варианты заданий к лабораторной работе №1

Примените методы численного интегрирования для вычисления следующих заданий.

1. ; 6. ;

2. ; 7. ;

3. ; 8. ;

4. ; 9. ;

5. ; 10. ;

Содержание отчета

Отчет должен содержать:

1. титульный лист

2. постановку задачи (согласно варианту)

3. краткое описание методов численного интегрирования

4. программную реализацию данных методов

5. выводы о проделанной работе.

Контрольные вопросы и задания

1. Какие методы численного интегрирования вы знаете?

2. Какой из методов численного интегрирования, в вашем случае, оказался наиболее точным и наименее точным?

3. Чему равна погрешность численного интегрирования для выше изложенных методов?

4. Запишите формулы для приближенного вычисления определенных интегралов.

5. Вычислите определенный интеграл при помощи методов численного интегрирования.

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

7. Сравните погрешность методов трапеции и центральных прямоугольников.

8. Как еще называется формула Симпсона и почему?

9. Запишите формулу для расчета погрешности.

10.* Запишите формулу Симпсона через линейную комбинацию формул трапеции и центральных прямоугольников.

Теорема Больцано–Коши.

Если непрерывная на отрезке функция на концах имеет противоположные знаки, т.е.

то на интервале она хотя бы один раз обращается в ноль.

Метод половинного деления

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

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

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

Рис.2.1. Последовательное деление отрезка пополам и приближение к корню

Поступая тем же образом и далее, получаем, что после делений длина отрезка, на котором лежит корень, сокращается в раз и становится равной (если корень не был точно определён на каком-то предыдущем этапе, то есть не совпал с при некотором ). Пусть – заданная точность, с которой требуется отыскать корень. Процесс деления отрезков следует остановить, как только станет верным неравенство . Очевидно, что если при этом положить в качестве корня

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

Метод секущих

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

и на каждой итерации нужно один раз вычислить значение функции .

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

Найдём точку пересечения этой прямой с осью из уравнения

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

Рис.2.3 Последовательные итерации метода секущих

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

Метод Ньютона

Рассмотрим эффективный метод решения нелинейных уравнений, носящий имя Ньютона. Вначале приведем некоторые наводящие рассуждения. Пусть функция y = F(x), корень которой ищется, имеет производные до 2-го порядка в окрестности корня - точки . Пусть уже найдено приближение номера n к корню (n-ая итерация) и требуется найти приближение номера n+1. По формуле Тейлора имеем

F(x_n+1) = F(x_n) + F’(x_n) × (x_n+1 – x_n) + O(x_n+1 – x_n)².

Пренебрежем остаточным членом порядка O(x_n+1 – x_n)² в правой части формулы и, будем считать, что x_n+1 » , т.е. приближение номера n+1 найдено столь точно, что F(x_n+1) » 0.

Тогда имеем приближенное равенство

0 » F(x_n) + F'(x_n) (x_n+1 – x_n).

Выражая отсюда x_n+1 при условии F'(x_n) ¹ 0, и, переходя от приближенного равенства к точному, получим

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

Предположим, что уравнение

F(x) = 0

(2.1)

имеет простой вещественный корень x_*, т.е.

F(x_*) = 0,

Будем предполагать, что F(x) дважды дифференцируема в некоторой окрестности точки x_*, т.е. для всех х принадлежащих некоторому интервалу (x_* - r₁, x_* + r₁), где r₁ > 0, причем F" (x) непрерывна на отрезке [x_* - r, x_* + r], 0 < r £ r₁.

Исследуем сходимость метода Ньютона

(2.2)

Теорема 1. Пусть x_* - простой вещественный корень уравнения (4.1) и пусть F'(x) ¹ 0 в окрестности точки. x_*

= {x: ½ x - x_*½ < r}.

Пусть, что F" (x) непрерывна на отрезке [x_*-r, x_*+r] Ì Ì , причем

(2.3)

Тогда, если и

(2.4)

то метод Ньютона (2.2) сходится, и для погрешности справедлива оценка

(2.5)

Замечания.

Метод Ньютона имеет квадратичную сходимость, т.е. он сходится быстрее метода простой итерации, который имеет линейную сходимость. Однако, метод Ньютона требует задания достаточно близкого к корню x_* начального приближения, удовлетворяющего неравенству (2.4) при соблюдении соотношений (2.3).

Пример выполнения лабораторной работы №2

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

Вводим функцию в исходном уравнении F(x)=0

Строим график данной функции

Из графика видно, что корень находится на интервале (0, 1). Вводим концы интервала.

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

Метод половинного деления.

Вводим функцию расчета нелинейных уравнений методом половинного деления

Вызываем данную функцию

Выводим найденное при помощи метода половинного деления приближенное значение корня и количество итераций

Считаем значение функции в данной точке (оно должно быть близким к нулю)

На графике функции отмечаем значение корня уравнения

Метод хорд.

Вводим функцию расчета нелинейных уравнений методом хорд

Вызываем данную функцию

Выводим найденное при помощи метода хорд приближенное значение корня и количество итераций

Считаем значение функции в данной точке

Метод секущих.

Вводим функцию расчета нелинейных уравнений методом секущих

Вызываем данную функцию

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

Считаем значение функции в данной точке

Метод Ньютона.

Считаем производную функции F(x)

Вводим функцию расчета нелинейных уравнений методом Ньютона

Вызываем данную функцию

Выводим найденное при помощи метода Ньютона приближенное значение корня и количество итераций

Считаем значение функции в данной точке

Варианты заданий к лабораторной работе №2

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

1. 6.

2. 7.

3. 8.

4. 9.

5. 10.

Содержание отчета

Отчет должен содержать:

1. титульный лист;

2. постановку задачи (согласно варианту);

3. краткое описание методов расчета нелинейных уравнений;

4. программную реализацию данных методов;

5. выводы о проделанной работе.

Контрольные вопросы и задания

1. Какие методы решения нелинейных уравнений вы знаете?

2. Какой из методов решения нелинейных уравнений, в вашем случае, оказался наиболее быстрым и медленным?

3. Дайте описание метода половинного деления.

4. Запишите расчетную формулу метода хорд.

5. Запишите расчетную формулу метода секущих.

6. Запишите расчетную формулу метода Ньютона.

7. Сформулируйте теорему Больцано–Коши.

8. Решите нелинейное уравнение.

9. Запишите формулу для расчета погрешности метода половинного деления.

10. Запишите формулу для расчета погрешности метода Ньютона.

Алгоритм LU-разложения.

Данный алгоритм можно рассматривать как конкретную форму метода Гаусса. Алгоритм LU-разложения используется не только для решения СЛАУ, но и также для обращения матрицы, т.е. вычисления матрицы, обратной данной.

Пусть и будем искать представление A в виде:

(3.10)

где L и U – соответственно нижняя и верхняя треугольные матрицы вида

Известно, что если все угловые миноры матрицы А отличны от нуля, т.е.

то разложение вида (3.10) существует и единственно. Для того чтобы получить расчётные формулы, поступим следующим образом. Обозначим a_ij произведение i-ой строки матрицы L на j-ый столбец матрицы R, причём будем считать вначале, что i< j.

Тогда .

Выразим из последней формулы u_ij.

(3.11)

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

В случае i = j имеем

Учитывая, что u_ii º 1 и, выражая из последнего соотношения l_ii, получаем:

(3.12)

Наконец, при i > j получаем

откуда, с учетом того, что u_jj º 1 приходим к формуле

(3.13)

Итак, расчетные формулы (3.11) – (3.13) получены. Для того чтобы при их применении не использовались неизвестные (не вычисленные) величины, необходимо выбрать соответствующий порядок вычисления элементов матриц L и U.

Например, можно рекомендовать порядок расчета элементов матриц L и U, схематически изображенный на рис. 3.1. На нем цифры слева для матрицы L и сверху - для матрицы U означают, что на первом шаге рассчитывается l₁₁ по формуле (3.12), затем вычисляется элемент u₁₂ по формуле (3.11).

Далее (3 шаг) определяются элементы второй строки матрицы L в порядке, указанном стрелкой: l₂₁ и l₂₂ (по формулам (3.13) и (3.12) соответственно).

На 4 шаге выполняется расчет элементов 3 столбца матрицы U в порядке, обозначенном стрелкой: u₁₃, u₂₃ (формулы (3.11)) и т.д.

Рис. 3.1.

Рассмотрим теперь применение LU-разложения для решения СЛАУ вида

Ах = b,

где

A = LU.

Введем вспомогательный вектор у,

у = U × x.

(3.14)

Тогда исходную систему можно записать так

L × Ux = b ó Ly = b.

(3.15)

В силу формул (3.14) и (3.15) решение исходной СЛАУ сводится к последовательному решению систем (3.15) и (3.14) соответственно с верхней и нижней треугольной матрицами.

Метод прогонки.

Метод прогонки представляет собой вариант метода Гаусса, примененный к специальным системам линейных алгебраических уравнений, и учитывающий ленточную структуру матрицы системы. Пусть имеем СЛАУ со специальной трехдиагональной формой матрицы

c₀y₀ – b₀y₁ = f₀, -a_iy_i-1 + c_iy_i – b_iy_i+1 = f_i, 1 £ i £ N-1; -a_Ny_N-1 + c_Ny_N = f_N,

(3.16)

или в матричной форме: AY = F, где Y = (y₀, y₁, ..., y_n)^T - вектор неизвестных; F = (f₀, f₁,..., f_n)^T - вектор правых частей; А - квадратная (N+1)´ (N+1) матрица

Системы вида (3.16) возникают при конечно-разностной аппроксимации краевых задач математической физики, описываемых обыкновенными дифференциальными уравнениями второго порядка с постоянными и переменными коэффициентами, а также уравнениями в частных производных. Ставится задача разработать экономичные методы решения задач вида (3.16), число арифметических операций для которых пропорционально числу неизвестных. Таким методом для системы (3.16) является метод прогонки. Специфика матрицы А состоит в расположении ненулевых элементов, матрица А - разреженная матрица, из (N+1)² элементов которой ненулевыми являются не более 3N+1 элементов. Это позволяет получить для решения СЛАУ простые расчетные формулы.

Будем искать решение (3.16) в виде

y_i = a_i+1y_i+1 + b_i+1, i = N-1, N-2, ..., 0

(3.17)

c неопределенными коэффициентами a_i, b_i. Выражение y_i-1 = a_iy_i + b_i подставим в (3.16)

(с_i - a_ia_i)y_i - b_iy_i+1 = f_i +a_ib_i

c учетом (3.17) имеем

[(с_i - a_ia_i)a_i+1 - b_i]y_i+1 + (c_i - a_ia_i)b_i+1 - a_ib_i = f_i.

Это равенство имеет место для любых y_i, если

(с_i - a_ia_i)a_i+1 - b_i = 0, (с_i - a_ia_i)b_i+1 - a_ib_i = f_i.

Отсюда получаем рекуррентные формулы для определения a_i+1, b_i+1

;	(3.18)
.	(3.19)

Коэффициенты a_i, b_i, i = 1, 2, ..., N называются прогоночными.

Если коэффициенты a_i и b_i известны, а также известно y_N то, двигаясь справа налево (от i+1 к i) последовательно определяем все y_i. Задача нахождения a_i, b_i по формулам (3.18), (3.19) решается слева направо (от i к i+1). Начальные значения прогоночных коэффициентов a_i, b_i можно определить следующим образом. Полагаем в формуле (3.17) i=0, имеем y₀=a₁y₁+b₁, а из первого уравнения (3.16) c₀y₀ - b₀y₁ = f₀, откуда

12 3 4 5 6