close

Вход

Забыли?

вход по аккаунту

?

Скачать оригинальный документ PDF (673.4 КБ)

код для вставкиСкачать
Ангарская государственная техническая академия
УДК 519.6(075)
Вычислительная математика. Часть первая: Учебное пособие для студентов
дневного и заочного обучения технических и химико-технологических специальностей./ В.С.Асламова, А.Г.Колмогоров, Н.Н.Ступакова. Ангарская государственная техническая академия. – Ангарск: АГТА, 2003г. - 82 с.
ISBN 5-89864-030-4
Вычислительная математика
Часть первая
В пособии рассматриваются основные положения численных методов, относящиеся к решению нелинейных уравнений, систем нелинейных и линейных алгебраических уравнений, интегрированию функций. Значительное внимание уделяется вопросам алгоритмизации методов.
Пособие может быть использовано при выполнении лабораторных, курсовых и дипломных исследовательских работ, так как содержит подробные
блок-схемы алгоритмов численных методов, с составлением которых у студентов чаще всего связаны основные трудности.
Рецензенты:
доктор физ.-мат. наук, старший научный
сотрудник Института динамики систем
и теории управления СО РАН
М.В. Булатов
кандидат физ.-мат. наук, доцент
кафедры математики АГТА
С.А.Чихачев
© Ангарская государственная техническая академия, 2003
© Кафедра автоматизации технологических процессов и производств
© В.С.Асламова, А.Г.Колмогоров, Н.Н.Сумарокова
Ангарск 2003 г.
Содержание
Введение........................................................................................................... 4
1. Элементы общей теории приближенных методов............................... 5
1.1. Источники и виды погрешности...................................................................... 5
Абсолютная и относительная погрешности ...................................................... 6
Погрешность функции......................................................................................... 7
Устойчивость, корректность, сходимость ......................................................... 7
1.2. Контрольные вопросы .................................................................................... 11
Введение
Данное учебное пособие содержит основы численных методов решения для
нелинейных уравнений,
систем нелинейных
и линейных уравнений,
дифференциальных уравнений, методы аппроксимации функций, обращения
матриц, вычисления определенных интегралов. Так же в него включена глава, в
которой рассматриваются элементы общей теории приближенных методов.
2. Методы решения нелинейных уравнений и систем.......................... 12
Основное предназначение пособия – облегчить работу преподавателя и
2.1. Графический метод отделения корней.......................................................... 12
2.2. Метод половинного деления (метод бисекции, метод дихотомии)............ 14
2.3. Метод простой итерации ................................................................................ 16
2.3.1. Решение нелинейных уравнений ............................................................ 16
2.3.2. Решение систем нелинейных уравнений методом итераций ............... 17
2.4. Метод Ньютона (метод касательных) ........................................................... 21
2.4.1. Решение нелинейных уравнений ............................................................ 21
2.5. Метод хорд ...................................................................................................... 26
2.6 Комбинированный метод секущих и хорд..................................................... 29
2.7. Контрольные вопросы .................................................................................... 32
повысить эффективность учебного процесса. Оно позволяет сформировать у
3. Методы вычисления определенных интегралов............................... 33
3.1. Методы прямоугольников............................................................................. 35
3.2. Метод трапеций. Вычисление значения интеграла с заданной точностью 38
3.3. Метод Симпсона ............................................................................................ 41
3.4. Метод Монте-Карло........................................................................................ 44
3.5. Использование сплайнов для численного интегрирования......................... 49
3.6. Погрешность численного интегрирования ................................................... 51
3.7. Контрольные вопросы .................................................................................... 55
студентов основные сведения о численных методах, необходимых для
первоначального ознакомления с предметом, привить навыки алгоритмизации
численных методов.
Пособие может быть использовано при выполнении лабораторных, курсовых
и дипломных исследовательских работ, так как содержит подробные блок-схемы,
с составлением которых у студентов чаще всего связаны основные трудности.
Каждая
глава
содержит
теоретическое
и
блок-схемы
рассматриваемых методов и завершается контрольными вопросами по данной
теме.
4. Методы решения систем линейных алгебраических уравнений .. 56
4.1. Метод Гаусса .................................................................................................. 57
4.2. Метод Гаусса с выбором главного элемента ................................................ 59
4.3. Метод Гаусса-Зейделя ................................................................................... 66
4.4. Метод прогонки.............................................................................................. 69
4.5. Метод Гаусса-Жордана................................................................................... 71
4.6. Вычисление определителя по методу Гаусса ............................................... 74
4.7. Метод итераций.............................................................................................. 77
4.9. Контрольные вопросы .................................................................................... 81
Литература .................................................................................................... 82
3
обоснование
4
1. Элементы общей теории приближенных методов
1.1. Источники и виды погрешности
На некоторых этапах решения задачи на ЭВМ могут возникать погрешности,
искажающие результаты вычислений. Оценка степени достоверности
получаемых результатов является важнейшим вопросом при организации
вычислительных работ. Это особенно важно при отсутствии опытных или других
данных для проверки адекватности модели.
Погрешность решения задачи обуславливается следующими причинами:
1. Математическое описание задачи является неточным, в частности,
неточно заданы исходные данные.
2. Применяемый для решения метод часто не является точным: получение
точного решения возникающей математической задачи требует
неограниченного или неприемлемо большого числа арифметических
операций, и, поэтому вместо получения точного решения задачи
приходится прибегать к приближенному.
3. При вводе данных в машину, при выполнении арифметических
операций и при выводе данных производятся округления. Погрешности,
соответствующие этим причинам, называют:
1) неустранимой погрешностью,
2) погрешностью метода,
3) вычислительной погрешностью.
Часто неустранимую погрешность разделяют на две части:
a) неустранимой погрешностью называют лишь погрешность,
являющуюся следствием неточности задания числовых данных, входящих в
математическое описание задачи;
в) погрешность, являющуюся следствием несоответствия
математического описания задачи реальности, называют погрешностью
математической модели.
Математическая модель, принятая для описания данного процесса или
явления, может внести существенные погрешности, если в ней не учтены какиелибо существенные черты
рассматриваемой
задачи. В частности,
математическая модель может прекрасно работать в одних условиях и быть
совершенно неприемлемой в других; поэтому важно правильно учитывать
область её применимости.
Численный метод также является источником погрешностей. Это связано,
например, с заменой интеграла суммой, усечением рядов при вычислениях
значений функций, интерполированием табличных данных и т. п. Как правило,
погрешность численного метода регулируема, т. е. она может быть уменьшена до
любого разумного значения путем изменения некоторого параметра (например,
шага интегрирования, числа членов усеченного ряда и т.п.). Погрешность метода
обычно стараются довести до величины, в несколько раз меньшей погрешности
исходных данных. Дальнейшее снижение погрешности не приведет к
5
повышению точности результатов, а лишь увеличит стоимость расчетов из-за
необоснованного увеличения объема вычислении.
При вычислениях с помощью ЭВМ неизбежны погрешности округлений,
связанные с ограниченностью разрядной сетки машины. Обычно после
выполнения операции производится не округление результата, а простое
отбрасывание лишних разрядов с целью экономии машинного времени. Правда,
в современных машинах предусмотрена свобода выбора программистом способа
округления; соответствующими средствами располагают и некоторые
алгоритмические языки.
Максимальная относительная погрешность при округлении есть δmax=0.5α1-k,
где α - основание системы счисления, k—количество разрядов мантиссы числа.
При простом отбрасывании лишних разрядов эта погрешность увеличивается
вдвое.
В современных машинах с памятью, измеряемой в байтах, принята
шестнадцатеричная система счисления, и любое число с плавающей точкой
содержит шесть значащих цифр. Следовательно, α =16, k = 6, максимальная
погрешность округления αmax= 0.5*16-5 ≈ 0.5*10-8.
Несмотря на то что при решении больших задач выполняются миллиарды
операций, это вовсе не означает механического умножения погрешности при
одном округлении на число операций, так как при отдельных действиях
погрешности могут компенсировать друг друга (например, при сложении чисел
разных знаков). Вместе с тем иногда погрешности округлений в сочетании с
плохо организованным алгоритмом могут сильно исказить результаты.
Перевод чисел из одной системы счисления в другую также может быть
источником погрешности из-за того, что основание одной системы счисления не
является степенью основания другой (например, 10 и 2). Это может привести к
тому, что в новой системе счисления число становится иррациональным.
Например, число 0.1 при переводе в двоичную систему счисления примет вид
0.1=0.00011001100… Может оказаться, что с шагом 0.1 нужно при вычислениях
пройти отрезок [0,1] от х=1 до х= 0; десять шагов не дадут точного значения х=0.
Абсолютная и относительная погрешности
Если А – точное значение некоторой величины, а а – известное приближение
к нему, то абсолютной погрешностью приближения а числа А называют
некоторую величину Δ(а) удовлетворяющую условию:
⏐А - а⏐≤ Δ(А).
Относительной погрешностью называют некоторую величину δ(А), для
которой выполняется условие:
а−А
≤ δ (а ) .
а
Относительную погрешность часто выражают в процентах.
6
Информацию о том, что а является приближенным значением числа А с
абсолютной погрешностью Δ(а), принято записывать в виде:
А = а ± Δ(а).
Числа а и Δ(а) записываются с одинаковым количеством знаков после
запятой.
Информацию о том, что а является приближенным значением числа А с
относительной погрешностью δ(а), записывают в виде:
А = а(1± δ(а)).
Погрешность функции
Пусть искомая величина Y является функцией параметров а1, а2,…, аn; т.е. Y
= Y(a), и известна область G в пространстве переменных а1, а2,…, аn, которой
принадлежат параметры. Необходимо получить приближение y к Y и оценить его
погрешность.
Если y – приближённое значение величины Y, то предельной абсолютной
погрешностью А(у) называют наилучшую при имеющейся информации оценку
погрешности величины у:
(1.1)
A ( y ) = sup Y ( a1 , a 2 ,..., a n ) − y .
a∈G
Предельной относительной погрешностью δ(у) называют величину
A( y ) .
δ ( y) =
y
Устойчивость, корректность, сходимость
Устойчивость. Рассмотрим погрешности исходных данных. Поскольку это
так называемые неустранимые погрешности и вычислитель но может с ними
бороться, то нужно хотя бы иметь представление об их влиянии па точность
окончательных результатов. Конечно, мы вправе надеяться на то, что
погрешность результатов имеет порядок погрешности исходных данных. Всегда
ли это так? К сожалению, нет. Некоторые задачи весьма чувствительны к
неточностям в исходных данных. Эта чувствительность характеризуется так
называемой устойчивостью.
Пусть в результате решения задачи по исходному значению величины х
находится значение искомой величины у. Если исходная величина имеет
абсолютную погрешность Δх, то решение имеет погрешность Δу. Если решение
непрерывно зависит от входных данных, т.е. всегда Δy → 0 при Δx → 0 , то
uxx+uyy = 0, u(x, 0) = 0, uy(x, 0) = ϕ(x).
Входными данными является ϕ(x) Если ϕ ( x) = 0 , то задача имеет только
тривиальное решение u ( x, y ) = 0 . Если же ϕ n ( x) =
1
cos nx , то решением
n
будет
1
cos nx ⋅ shny .
n2
Очевидно, ϕ n (x) равномерно сходятся к ϕ (x) при n → ∞ ; но при этом, если
y ≠ 0 , то un (x,y) неограниченно и никак не может сходиться к u ( x, y ) . Этот
u n ( x, y ) =
пример связан с физической задачей о тяжелой жидкости, налитой поверх
легкой; при этом действительно возникает релей-тейлоровская неустойчивость.
Отсутствие устойчивости обычно означает, что даже незначительные
погрешности в исходных данных приводят к большим погрешностям в решении
или вовсе к неверному результату. О подобных устойчивых задачах также
говорят, что они чувствительны к погрешностям исходных данных.
Примером такой задачи является отыскание действительных корней
уравнения вида:
tg(x) – x = 1.
Изменение аргумента х на малую величину Δх может привести к
существенному изменению функции.
Корректность. Задача называется поставленной корректно, если для любых
значений исходных данных из некоторого класса ее решение существует,
единственно и устойчиво по исходным данным.
Рассмотренные выше примеры неустойчивой задачи являются некорректно
поставленными. Применять для решения таких задач численные методы, как
правило, нецелесообразно, поскольку возникающие в расчетах погрешности
округлений будут сильно возрастать в ходе вычислений, что приведет к
значительному искажению результатов.
Вместе с тем отметим, что в настоящее время развиты методы решения
некоторых некорректных задач. Это в основном так называемые методы
регуляризации. Они основываются на замене исходной задачи корректно
поставленной задачей. Последняя содержит некоторый параметр, при
стремлении которого к нулю решение этой задачи переходит в решение
исходной задачи.
задача называется устойчивой по входным данным.
Рассмотрим классический пример неустойчивости – задачу Коши для
эллиптического уравнения в полуплоскости y≥0:
Неустойчивость методов. Иногда, даже если задача устойчива, то
численный алгоритм может быть неустойчивым. Например, если производные
заменяются разностями, то приходится вычитать близкие числа и сильно
теряется точность. Эти неточные промежуточные результаты используются в
дальнейших вычислениях, и ошибки могут сильно нарастать.
7
8
Рассмотрим еще один пример неустойчивого
численный метод вычисления интеграла
алгоритма.
1
∫
I n = x n e x −1dx ,
n=1,2,…
0
Интегрируем по частям, для этого используем следующую формулу:
b
I = ∫ uvdu
a
b
= uv a
b
− ∫ vdu
a
Получаем
1
1
∫
1
∫
I1 = xe x −1dx = xe x −1 − e x −1dx =
0
0
1
∫
0
1
,
e
1
1
∫
I 2 = x 2 e x −1dx = x 2 e x −1 − 2 xe x −1dx = 1 − 2 I1 ,
0
0
0
LLLLLLLLLLLLLLLLLLL
1
∫
1
1
∫
I n = x n e x −1dx = x n e x −1 − n x n −1e x −1dx = 1 − nI n −1.
0
0
0
Пользуясь полученным рекуррентным соотношением, вычисляем:
I1=0.367879,
I2=0.263242,
I3=0.207274,
I4=0.170904,
I5=0.145480,
I6=0.127120,
I7=0.110160,
I8=0.118720,
I9=-0.0684800
9
Построим
Значение интеграла I9 не может быть отрицательным, поскольку
подынтегральная функция x9ex-1 на всем отрезке интегрирования [0,1]
неотрицательна. Исследуем источник погрешности. Видим, что округление в I1
дает погрешность, равную примерно лишь 4.4*10-7. Однако на каждом этапе эта
погрешность умножается на число, модуль которого больше единицы (-2, -3,...,9), что в итоге дает 9!. Это и приводит к результату, не имеющему смысла. Здесь
снова причиной накопления погрешностей является алгоритм решения задачи,
который оказался неустойчивым.
Численный алгоритм (метод) называется корректным в случае
существования и единственности численного решения при любых значениях
исходных данных, а также в случае устойчивости этого решения относительно
погрешностей исходных данных.
Понятие сходимости. При анализе точности вычислительного процесса
одним из важнейших критериев является сходимость численного метода. Она
означает близость получаемого численного решения задачи к истинному
решению. Строгие определения разных оценок близости могут быть даны лишь с
привлечением аппарата функционального анализа. Здесь мы ограничимся
некоторыми понятиями сходимости, необходимыми для понимания
последующего материала.
Рассмотрим понятие сходимости итерационного процесса. Этот процесс
состоит в том, что для решения некоторой задачи и нахождения искомого
значения определяемого параметра (например, корня нелинейного уравнения)
строится метод последовательных приближений. В результате многократного
повторения этого процесса (или итераций) получаем последовательность
значении х1, х2, …, хn,… эта последовательность сходится к точному решению
х=a, если при неограниченном возрастании числа итераций предел этой
последовательности существует и равен a: lim xn = a . B этом случае имеем
n→∞
сходящийся численный метод.
Другой подход к понятию сходимости используется в методах дискретизации.
Эти методы заключаются в замене задачи с непрерывными параметрами на
задачу, в которой значения функций вычисляются в фиксированных точках. Это
относится, в частности, к численному интегрированию, решению
дифференциальных уравнений и т. п. Здесь под сходимостью метода
понимается стремление значений решения дискретной модели задачи к
соответствующим значениям решения исходной задачи при стремлении к нулю
параметра дискретизации (например, шага интегрирования).
Таким образом, для получения решения задачи с необходимой точностью, ее
постановка должна быть корректной, а используемый численный метод должен
обладать устойчивостью и сходимостью.
10
1.2. Контрольные вопросы
1. Перечислите виды погрешностей.
2. Чему равна предельная относительная погрешность произведения или
частного?
3. Назовите требования к оценкам точности алгоритма.
4. Понятие сходимости, устойчивости и корректности приближённого
метода.
5. Назовите причины возникновения погрешностей.
6. Назовите единицы измерения абсолютной и относительной погрешности.
7. Может ли погрешность быть отрицательным числом?
8. Какая погрешность позволяет судить о качестве произведенных
измерений?
2. Методы решения нелинейных уравнений и систем
Задача нахождения корней нелинейных уравнений вида f(x)=0 встречается в
различных областях научных исследований (здесь f(x) - некоторая непрерывная
функция). Нелинейные уравнения можно разделить на два класса алгебраические и трансцендентные. Алгебраическими уравнениями называются
уравнения, содержащие только алгебраические функции (целые, рациональные,
иррациональные). В частности, многочлен является целой алгебраической
функцией. Уравнения, содержащие другие функции (тригонометрические,
показательные, логарифмические и др.), называются трансцендентными.
Методы решения нелинейных уравнений делятся на прямые и итерационные.
Прямые методы позволяют записать корни в виде некоторого конечного
соотношения (формулы).
Однако встречающиеся на практике уравнения не удается решить такими
простыми методами. Для их решения используются итерационные методы, т.е.
методы последовательных приближений. Алгоритм нахождения корня уравнения
с помощью итерационного метода состоит из двух этапов:
а) отыскания приближенного значения корня или содержащего его отрезка;
б) уточнения приближенного значения любым итерационным методом до
некоторой заданной степени точности.
Приближённое значение корня (начальное приближение) может быть найдено
различными способами: из физических соображений, из решения аналогичной
задачи при других исходных данных, с помощью графических или
аналитических методов.
Если такие априорные оценки исходного приближения к корню провести не
удаётся, то находят две близко расположенные точки a и b, в которых
непрерывная функция f(x) принимает значения разных знаков, то есть f(a)*f(b)<0.
В этом случае между точками a и b есть, по крайней мере, одна точка, в которой
f(x)=0. В качестве начального приближения x0 можно принять середину отрезка
[a, b], то есть x0=(a+b)/2.
Итерационный процесс состоит в последовательном уточнении начального
приближения x0. Каждый такой шаг называется итерацией. В результате
итераций находится последовательность приближённых значений корня x1, x2,...,
xn. Если эти значения с ростом n приближаются к истинному значению корня, то
говорят, что итерационный процесс сходится.
2.1. Графический метод отделения корней
Отделение корней – отыскание начального приближения корней, сводится к
отысканию достаточно малых областей, в которых находится только один
корень.
Для отделения корней полезна известная теорема из математического
анализа.
11
12
Теорема 1. Если непрерывная функция принимает значения разных знаков на
концах отрезка [a, b], т.е. f(a)*f(b) < 0, то внутри этого отрезка содержится
по меньшей мере один корень уравнения f(x)=0.
Корень заведомо будет единственным, если производная f '(x) существует и
сохраняет постоянный знак внутри интервала (а, b).
Для отделения корней можно воспользоваться графическим путем. Для этого
представим уравнение f(x)= 0 в виде x = ϑ(x), и строим графики функций у1 = x
и y2 = ϑ(x). Тогда абсциссы точек пересечения этих кривых и будут
действительными корнями уравнения f(x) = 0.
Пример. Отделить корни уравнения
x3 - 4x + 2 = 0.
Перепишем уравнение в виде
x(x2 - 4) + 2 = 0; x = -2 / (x2 - 4).
Построим графики функций
y1 = x;
y2 = -2 / (x2 - 4).
Абсциссы точек пересечения кривых y1, y2 дают приближенные значения
корней x2 ≈ 0,5; x2 ≈1,6; x3 ≈ -2,2 (рис.1).
y1
2.2. Метод половинного деления
(метод бисекции, метод дихотомии)
Этот метод требует непрерывности функции f(x), кроме того, на концах
отрезка [a, b] функция должна иметь разные знаки.
a+b
В качестве начального приближения принимаем середину отрезка c =
.
2
Для определения отрезка, содержащего корень, проверяем условие:
f(a)⋅f(с)<0
(2.1)
Можно сравнивать значения функции на концах отрезка [с, b].
Если условие (2.1) выполняется, то корень лежит в отрезке [a, c], поэтому
нужно передвинуть правую границу b в точку c , т.е. b=c.
Если условие (2.1) не выполняется, то корень лежит в отрезке [с, b] и следует
поменять левую границу a=c.
Проверяем условие выхода из цикла: «длина отрезка меньше или равна
заданной точности?»:
|b-a|≤ε
(2.2)
f(x)
f(x)
f(b)
f(c1 )
2
0
1
2
1
1
1
2
2
a
х
c
c c
b
f(c0 )
f(a)
Рис.2 Геометрическая интерпретация метода
половинного деления
y2
Рис.1.
Таким образом, функция имеет три корня, каждый из которых нужно
уточнить выбранным численным методом. Для уточнения корня выбираем
поочерёдно три интервала x1∈[-2.5; -2.05], x2∈[0; 0.8], x3∈[1.2; 1.8].
13
Если условие (2.2) не выполняется, то снова делим отрезок пополам,
определяем отрезок, содержащий корень и т.д. После каждой итерации отрезок,
содержащий корень уменьшается вдвое, т.е. после n итераций он сократится в 2n
раз. На рис.2 приведена графическая интерпретация метода половинного
деления, а на рис.3 – блок-схема метода половинного деления.
Метод половинного деления всегда сходится, то есть при его использовании
решение получается всегда, причём с заданной точностью. Требуемое обычно
большее число итераций по сравнению с некоторыми другими методами не
является препятствием к применению этого метода, если каждое вычисление
значения функции f(x) несложно.
14
2.3. Метод простой итерации
Начало
На концах отрезка
функция имеет разные
знаки?
2.3.1. Решение нелинейных уравнений
Для использования этого метода исходное нелинейное уравнение f(x)=0
преобразуется к виду:
x=ϕ(x).
(2.3)
Пусть известно начальное приближение корня x0. Подставляя это значение в
правую часть уравнения (2.3) получаем новое приближение: х1=ϕ(х0). Далее,
подставляя каждый раз новое значение корня, получаем последовательность хn+1
= ϕ (хn), n = 1, 2, ....
Итерационный
процесс
прекращается,
если
результаты
двух
последовательных итераций близки:
⏐хn+1−хn⏐ ≤ ε.
(2.4)
Достаточным условием сходимости метода простой итерации является
выполнение условия
dϕ ( x)
(2.5)
<1
ввод
a, b, ε
f(a)*f(b)<0
нет
да
n: = 0
c : = (a + b)/2
n: = n + 1
выбор отрезка,
содержащего корень
нет
f(a)*f(c)<0
a: = c
да
b: = c
нет
проверка условия
выхода из цикла
⎟b – a⎟ ≤ ε
да
вывод n,
c, f(c)
dx
для любого х, принадлежащего отрезку [a, b], который содержит корень.
На рис.6 представлена геометрическая интерпретация метода итераций. При
начальном приближении х0 последовательность приближений хn сходится к
корню х*. В случае выполнения условия (2.5) в качестве начального
приближения можно выбрать любую точку интервала [a, b].
Если мы попытаемся уточнить корень х* и выберем в качестве начального
приближения х0, то, как видно из рис.6, последовательность приближений хn (n =
0,1,…) будет расходящейся; расстояние между двумя последовательными
приближениями будет возрастать с увеличением числа итераций n.
Если же процесс расходится, т.е. не выполняется условие (2.5), то уравнение
f(x)=0 умножаем на константу l и к левой и правой части уравнения прибавляем
x, получаем x = x + lf(x). Новая ϕ(x) = x + lf(x). Выберем константу l таким
образом, чтобы обеспечить выполнение условия (2.5). Положим ϕ '(x)=0.5.
ϕ '(x) = 1 + lf'(x) = 0.5,
отсюда l = − 0 . 5 ,
u
Конец
Рис.3. Блок-схема метода половинного деления
15
где u соответствует max {⏐f '(a)⏐,⏐f '(b)⏐}, т.е. если ⏐f '(a)⏐> ⏐f '(b)⏐ то u = f'(a),
иначе u=f '(b) (блок-схема на рис. 4).
Стоит отметить, что приведенный способ преобразования уравнения
справедлив, в основном, для монотонных функций.
В блок-схемах: с - начальное приближение корня, а в дальнейшем - результат
предыдущей итерации, х - значение корня после каждой итерации.
16
2.3.2. Решение систем нелинейных уравнений методом итераций
Начало
Пусть требуется найти корни системы нелинейных уравнений вида:
ввод
a, b, ε
x1 = f1 (x1, x2,…, xn),
x 2= f2 (x1, x2,…, xn),
..………………….
xi = fi (x1, x2,…, xn)
..………………….
xn = fn (x1, x2,…, xn).
(2.6)
Алгоритм решения этой системы методом простой итерации напоминает
метод Гаусса-Зейделя, используемый для решения систем линейных уравнений.
Пусть в результате предыдущей итерации получены значения неизвестных
x1 = a1, x2 = a2, …, xn = an.
⎟f1⎟>⎟f2⎟
нет
l: = -0,5 / f2
да
l: = -0,5/f1
n: = 0;
c: = a;
x: = b
Начало
ввод
a, b, ε
проверка сходимости
процесса итерации
⏐ϕ’(a)⎟< 1
and
⏐ϕ’(b)⎟< 1
да
нет
⏐x-c⏐> ε
⎟ x-c⎟>ε
нет
да
нет
c: = a
x: = b
c: = x
x: = c + l*f (c)
печать ‘процесс итерации
расходится, выберите
блок-схему №2 (рис.5)’
n: = n + 1
exit
печать
‘корень =’,c,
f(c), n
да
n: = n +1
f1: = f ’(a)
f2: = f ’(b)
c: = x
x: = ϕ(c)
Конец
Рис.5. Блок-схема решения нелинейного уравнения методом
простой итерации №2.
печать c, f(c),
n.
Конец
Рис.4. Блок-схема решения нелинейных уравнений методом
простой итерации №1.
17
18
Начало
Ввод
k ,m ,n ,ε
к – счетчик итераций
m – максимальное число итераций
n – число уравнений в системе
ε - точность
Ввод
массива а
y
k:= 0
ϕ(x1
1
ϕ(x2
a x*
x2
x1
ϕ(x1
x0
Предыдущее приближение
корней запоминаем в
массиве b
d:= 0
ϕ(x0)
b:=a
ϕ(x0
x[1]:=f1(a)
a[1]:=x[1]
x[2]:=f2(a)
a[2]:=x[2]
………….
x
b
x**x0 x1
Расчет следующего
приближения по формуле (2.7)
х2
Рис.6. Геометрическая интерпретация метода
итераций
Восстанавливаем массив
предыдущих приближений
a:=b
определяем
максимальное
приращение
Тогда выражения для неизвестных на следующей итерации имеют вид:
x1 = f1(a1, a2,…, an),
x 2 = f2(x1, a2,…, an),
.….......................
xi = fi(x1,…, xi-1,ai,…,an)
.….......................
xn = fn(x1,…, xn-1, an).
i := 1, n
d1 := ⎜x[i]-а[i]⎜
(2.7)
Итерационный процесс продолжается до тех пор, пока изменения всех
неизвестных в двух последовательных итерациях не станут малыми, то есть
абсолютные величины их разностей не станут меньшими заданного малого
числа.
При использовании метода простой итерации успех во многом определяется
удачным выбором начальных приближений неизвестных: они должны быть
достаточно близкими (несколько долей единицы) к истинному решению. Блоксхема метода простой итерации для системы нелинейных уравнений
представлена на рис. 7.
нет
d < d1
да
d := d1
k := k+1
Переопределяем
предыдущее
приближение
1
нет
d≤ε
да
нет
а := x
k>m
да
Печать
корней х
Печать
‘процесс
расходится’
Конец
exit
Рис.7 Блок-схема метода простой итерации для системы нелинейных
уравнений
19
20
2.4. Метод Ньютона (метод касательных)
2.4.1. Решение нелинейных уравнений
Пусть известно, что корень уравнения f(x)=0 принадлежит отрезку [a, b].
Выбирается начальное приближение с0 к корню и через точку М0 с
координатами (с0, f(c0)) проводится касательная к кривой y=f(x). Точка
пересечения с1 касательной с осью х принимается за следующее приближение к
корню. Проводя касательную через точку М1 найдём следующее приближение с2
и т.д.
Начальное приближение с0 выбирается из условия, чтобы знак функции
совпадал со знаком кривизны, определяемой второй производной, т.е.
f(c0)⋅f ′′(c0)>0
(2.8)
Условие (2.8) вытекает из теоремы 2.
Теорема 2. Если f(a)*f(b) < 0, причем f '(x) и f ''(x) отличны от нуля и
сохраняют определенный знак при а≤ х≤ b, то исходя из начального
приближения с0∈ [a, b] удовлетворяющего неравенству (2.8), можно вычислить
(n=1,2…)
методом Ньютона с = с − f (с n −1 )
n
n −1
f ′(с n −1 )
единственный корень уравнения f(x)=0 с любой степенью точности ε.
Если условие (2.8) выполняется для точки а, то с0=а; если условие (2.8)
выполняется для точки b, то с0=b.
Уравнение касательной, проведённой к кривой y=f(x) через точку М0 имеет
следующий вид:
y-f(c0)=f′(c0)(x-c0).
В точке х=с1 y(c1)=0, тогда формула для определения следующего
приближения к корню с1 через предыдущее с0 примет следующий вид:
f (с0 )
(2.9)
с =с −
1
0
f ′(с0 )
Для окончания итерационного процесса может быть использовано условие:
⎟f(c1)⎜<ε
(2.10)
или условие близости двух соседних приближений:
⎟c1 - с0⎜<ε
(2.11)
y
M0
Начало
ввод
k, m, n
i : = 1, n
ввод
x[i]
расчет первых приближений
по формулам (2.2)
x[i] = ϕi(x01,…, x0n)
да
⏐ x0i- x1i⏐≤ ε
печать
корней
нет
k: = k + 1
k≤m
нет
‘процесс
расходится’
да
i: = 1, n
x0i: = x1i
y=f(x)
0
c1
c2
c0
x
M2
Конец
Рис. 9. Блок-схема метода простой итерации для
системы нелинейных уравнений.
M1
Рис. 8. Геометрическая интерпретация метода
21
22
На каждой итерации объём вычислений в методе Ньютона больший, чем в
других методах, поскольку приходится находить не только значение функции
f(x), но и значение её производной. Однако скорость сходимости метода Ньютона
значительно выше, чем в других методах.
Трудность в применении метода Ньютона состоит в выборе начального
приближения с0. Поэтому, иногда целесообразно использовать смешанный
алгоритм. Он состоит в том, что сначала применяется всегда сходящийся метод
(например, метод половинного деления), а после некоторого числа итераций −
быстро сходящийся метод Ньютона.
Блок-схема метода Ньютона для нелинейных уравнений представлена на
рис. 10.
Задание. Модифицируйте блок-схему на рис.10 так, чтобы было
использовано условие 2.11.
Метод обладает гораздо более быстрой сходимостью, чем метод простой
итерации. В основе метода Ньютона для системы уравнений (2.6) лежит
использование разложения функций fi (x1, x2,…, xn) в ряд Тейлора, причём члены,
содержащие вторые (и более высоких порядков) производные, отбрасываются.
Пусть приближённые значения корней системы уравнений (2.6) (например,
полученные на предыдущей итерации) равны соответственно а1, а2, …, аn. Задача
состоит в нахождении приращений (поправок) к этим значениям Δx1, Δx2, …,
Δxn, благодаря которым решение системы (2.7) запишется в виде:
x1 = a1 +Δ x1,
x2 = a2 + Δx2,
(2.12)
.....................
xn = an + Δxn.
Проведем разложение левых частей уравнений первой системы в ряд
Тейлора, ограничиваясь лишь линейными членами относительно приращений:
∂f
∂f1
Δx1 + ... + 1 Δx n
∂x n
∂x1
f 2 ( xi ,..., x n ) = f 2 (ai ,..., a n ) +
∂f 2
∂f
Δx1 + ... + 2 Δx n
∂x1
∂x n
.................................................................................
∂f
∂f
f n ( xi ,..., x n ) = f n (ai ,..., a n ) + n Δx1 + ... + n Δx n
∂x1
∂x n
23
n: = 1
выбор начального
приближения
ввод
a, b, ε
да
f(a)*f’’(a) > 0
с0: = a;
n: = 0
нет
нет
2.4.2. Решение систем нелинейных уравнений методом Ньютона
f1 ( xi ,..., x n ) = f1 (ai ,..., a n ) +
Начало
f(b)*f’’(b) > 0
да
c0: = b;
n: = 0
нет
n: = 0
да
с0: = c0 – (F(c0)/F′(c0))
n: = n + 1
нет
⏐f(c0)⏐< ε
да
вывод n,
c0, F(c0)
(2.13)
Конец
Рис. 10. Блок-схема метода Ньютона для нелинейных уравнений.
24
Поскольку в соответствии с системой нелинейных уравнений левые части
этих выражений должны обращаться в нуль, то приравняем нулю и правые
части. Получим следующую систему линейных алгебраических уравнений
относительно приращений:
y
1
∂f
∂f 1
∂f
Δx1 + 1 Δx 2 + ... + 1 Δx n = − f1
∂x1
∂x 2
∂x n
∂f 2
∂f
∂f
Δx1 + 2 Δx 2 + ... + 2 Δx n = − f 2
∂x1
∂x 2
∂x n
(2.14)
..........................................................
∂f
∂f
∂f n
Δx1 + n Δx 2 + ... + n Δx n = − f n
∂x n
∂x 2
∂x1
Значения f1, f2,…, fn и их производные вычисляются при x1=a1, x2=a2,…, xn=an.
Определителем системы является якобиан
∂ f1
⎡ ∂ f1
⎢ ∂ x ... ∂ x
n
⎢ 1
⎢ ∂ f 2 ... ∂ f 2
J = ⎢ ∂ x1 ∂ x n
⎢
....
⎢ ∂f
∂f
n
... n
⎢
⎢⎣ ∂ x1 ∂ x n
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
Для существования единственного решения системы он должен быть
отличным от нуля на каждой итерации.
Таким образом, итерационный процесс решения системы уравнений методом
Ньютона состоит в определении приращений Δx1, Δx2, …,Δxn к значениям
неизвестных на каждой итерации. Счет прекращается, если все приращения
становятся малыми по абсолютной величине: max ⏐Δxi⏐< ε.
1≤ i ≤ n
1
-1
x
-1
Рис.11. Отделение корней графическим
методом.
Из графика видно, что система имеет два решения. Уточним одно из них,
принадлежащее области D: 0,3<x<0,5; -0,8<y<-0,5. За начальное приближение
принимаем х0=0,3; у0=-0,8. Каждую пару корней следует уточнять отдельно,
используя какой-либо численный метод.
На рис. 13 представлена блок-схема метода Ньютона для системы из двух
нелинейных уравнений:
f1(x, y) = 0
f2(x, y) = 0
Значение приращений Δх, Δy определялось методом Крамера. В программе
следует записать формулы для расчета производных и функций в виде
подпрограммы FUNCTION.
2.5. Метод хорд
В методе Ньютона также важен удачный выбор начального приближения для
обеспечения хорошей сходимости. Метод имеет квадратичную сходимость, т.е.
погрешность на следующей итерации пропорциональна квадрату погрешности
на предыдущей итерации [8].
Очень часто не удаётся найти решения системы нелинейных уравнений даже
в случае двух переменных. Успешный поиск зависит, прежде всего, от выбора
начальных приближений к корням. Их следует выбирать вблизи корней. Пусть
дана следующая система уравнений:
⎧sin( 2 x − y ) − 1,2 x = 0,4
⎨
2
2
⎩0,8 x + 1,5 y = 1
Найдём корни графическим методом (рис.11).
Пусть мы нашли отрезок [a, b], на концах которого функция f(x) имеет разные
знаки. Для определенности примем f(a) > 0, f(b) < 0. Для данной функции
f ′′(х)<0 для любого х (рис.12). В методе хорд одна из точек фиксируется, а
положение другой постоянно переносится в точку ci пересечения хорды с осью х.
Фиксируется та точка, для которой знак функции совпадает со знаком второй
производной. Если выполняется условие
f(a)⋅f′′(a)>0 ,
(2.15)
то фиксируется точка а , если выполняется условие
f(b)⋅f′′( b)>0 ,
(2.16)
то фиксируется точка b.
В данном методе процесс итераций состоит в том, что в качестве
приближений к корню уравнения принимаются значения с0, с1, ... точек
25
26
пересечения хорды с осью абсцисс. Последовательность с0, с1, ... будет сходиться
к корню, если первая производная f′(х) и вторая производная f′′(х) сохраняют
постоянный знак на отрезке [a,b] [7].
f(x)
f(a)
А
Начало
ввод максимального
числа итераций m
и точности ε
f ′′<0
ввод начального
числа приближений
к корням х0, у0
n: = 1
0
a
c0
det :=
c1 b
B
Рис.12. Геометрическая интерпретация метода хорд
Уравнение хорды АВ:
y − f (b )
x−b .
=
f (b ) − f ( a ) b − a
Для точки пересечения хорды с осью абсцисс (x = c0, y = 0) получим
уравнение:
(b − a) ⋅ f (b) .
(2.17)
c0 = b −
f (b) − f (a)
Далее проверяется условие выхода из итерационного цикла (разность между
двумя соседними приближениями должна быть не больше заданной точности ε):
(2.18)
c0 − b ≤ ε .
Если условие (2.18) не выполняется, то точку b переносим в точку с0 (b=c0),
снова проводим хорду и по формуле (2.17) определяем следующее приближение
к корню.
Если функция имеет вид (рис.12), то для такой функции f''(x)<0 для любого x
из отрезка [a, b]. В этом случае фиксируется точка b и формула для расчёта точек
ci (i=0,1,...) примет вид:
(b − a) ⋅ f (a) .
(2.19)
c0 = a −
f (b) − f (a)
Условие выхода из цикла:
(2.20)
c0 − a ≤ ε .
∂f 1
∂f
∂f
∂f
( x , y ) * 2 ( x , y ) − 2 ( x 0 , y 0 ) * 1 ( x0 , y 0 )
∂x 0 0 ∂y 0 0
∂x
∂y
∂f 1
∂f
( x , y ) − f 1 ( x 0 , y 0 ) * 2 ( x 0 , y 0 )) / det
∂y 0 0
∂y
∂f 2
∂f 1
Δy := ( f 1 ( x 0 , y 0 ) *
( x , y ) − f 2 ( x0 , y 0 ) *
( x , y )) / det
∂x 0 0
∂x 0 0
n: = n + 1
Δx := ( f 2 ( x 0 , y 0 ) *
уточнение
корней
x0: = x0 + Δx
y0: = y0 + Δy
(⏐Δx⏐≤ ε)
and
(⏐Δy⏐≤ ε)
да
печать корней
x0, y0, n,
f1(x0, y0), f2(x0, y0)
нет
да
n≤m
нет
печать ‘метод не
сошелся за m
итераций’
exit
Конец
Рис.13. Блок-схема метода Ньютона системы из двух уравнений
Если условие (2.20) не выполняется, то переопределяем точку а=с0, затем
ищем новую точку c0 по формуле (2.19) и т.д.
27
расчет приращений
к корням
28
Блок-схема метода хорд представлена на рис.16. Алгоритмы метода
половинного деления и метода хорд похожи, однако второй из них в ряде
случаев дает более быструю сходимость итерационного процесса. При этом
успех его применения, как и в методе половинного деления, гарантирован.
Задание. Разработайте алгоритм, в котором бы проверялось условие выхода
из цикла (2.20).
2.6 Комбинированный метод секущих и хорд
Этот метод обладает гарантированной сходимостью. Вводятся два
приближения х0 и х1. Если х0 и х1 лежат по разные стороны от корня, то
проводится хорда, в противном случае проводится секущая.
Случай, когда проводится хорда, рассмотрен на рис.14.
f(x)
D
f(x )
1
0
f ( x1 )x0 − f (x0 )x1 − f (x0 )x0 + f (x0 ) x0 x0 ( f ( x1 ) − f (x0 )) − f ( x0 )(x1 − x0 ) ,
=
f (x1 ) − f ( x0 )
f (x1 ) − f ( x0 )
тогда точка пересечения хорды с осью ОХ находится по формуле:
f ( x 0 ) ⋅ ( x1 − x 0 ) .
(2.21)
x = x0 −
f ( x1 ) − f ( x 0 )
Случай, когда проводится секущая, изображен на рис.15.
x=
y
B
f(x0)
x0
x
A
C
Рис.15
E
и
АDЕ,
вытекает пропорциональность
BC AC
соответственных сторон
, значит
=
DE AE
C
Рис.14.
Так как ΔАBC подобен ΔBDЕ, то длины соответственных сторон
пропорциональны:
AC AB ,
=
ED BE
то есть
− f ( x0 ) x − x0 .
=
f ( x1 )
x1 − x
Отсюда получаем
-f(x0)x1 + f(x0)x = f(x1)x + f(x1)x0 .
f ( x0 ) x0 − x .
=
f ( x1 ) x1 − x
Нетрудно вывести, что конечная формула для определения точки х и в этом
случае совпадает с формулой (2.21).
После определения х проверяется условие выхода из цикла:
x − x1 ≤ ε .
(2.22)
Если оно не выполняется, то в качестве точки х0 выбирается точка х1 (х0 := х1),
а в качестве точки х1 берём точку х (х1 := х). Снова проводим секущую или хорду,
рассчитываем следующее приближение х по формуле (2.21) и т.д.
Блок-схема комбинированного метода секущих и хорд представлена на
рис.17.
Выделим из этого уравнения х:
x=
x1
E
x1
B
x
Из подобия треугольников АBC
f(x0)
D
f(x1)
0
A
x0
y = f(x)
f ( x1 ) x0 − f ( x0 ) x1 .
f ( x1 ) − f ( x0 )
Чтобы упростить выражение, добавим к числителю f(x0)x0 и его же отнимем:
29
30
Начало
Начало
ввод
a, b, ε
Печать
‘ввод
ε, х1, х’
ввод
n: = 1
да
f(a)*f ’’(a)>0
нет
нет
f(b)*f (b)>0
да
d: = a;
a: = b;
b: = d
n: = 0
печать
' точка b
фиксируется’
нет
n: = 0
печать
' точка а
фиксируется’
’’
ε, х1, х
выбор
фиксированной точки
n: = 0
Переобозначение точек.
Для разработки
универсальной программы,
пригодной для любой
функции в случае
фиксации точки
а, достаточно
переобозначить точки.
x0: = x1
x1: = x
x: = x1-(f(x1)( x1- x0))/(f(x1)-f(x0));
n: = n + 1
нет
⎮ x- x1⎮≤ ε
да
печать
x, f(x), n
n: = 0
Конец
да
c: = a;
a: = b
⎟c - a⎜> ε
Рис. 17. Блок-схема комбинированного метода секущих и хорд.
2.7. Контрольные вопросы
нет
да
a: = c
c: = a –(f(a)*(b – a))/(f( b)-f( a))
n: = n + 1
печать
‘ корень
x=’,c,
‘n =’,n
‘f(x) =’,f(c)
Конец
Рис. 16. Блок-схема метода хорд.
31
1. В чём состоит отличие алгебраического уравнения от трансцендентного?
2. Сущность и физический смысл процедуры отделения корней.
3. Обладает ли метод половинного деления гарантированной сходимостью?
4. Как выбирается начальное приближение в методе Ньютона?
5. Для каких функций не рекомендуется применять метод Ньютона?
6. Модификация метода Ньютона. Его особенности и случаи применения.
7. Может ли в методе хорд интервал находиться с одной стороны от корня?
8. Назовите условие выбора интервала в методе секущих.
9. Назовите достоинства комбинированного метода секущих и хорд.
10. К какому виду нужно преобразовать уравнение для метода итераций?
11. Можно ли воспользоваться методом итераций при невыполнении условия
сходимости?
32
3. Методы вычисления определенных интегралов
Требуется вычислить определенный интеграл:
Вычисление суммы (3.2) дает простейший пример численного
интегрирования. А верхняя S2 и нижняя S1 суммы Дарбу определяют величину
погрешности S, а именно:
I − S ≤ S 2 − S1 , где
b
I =∫
f ( x)dx .
(3.1)
n −1
S1 = ∑ mi ( xi +1 − xi ); mi = min f ( x ),
a
Геометрический смысл простейшего определенного интеграла от
неотрицательной функции f(x)≥0 состоит в том, что значение I – это площадь,
ограниченная кривой y=f(x), осью абсцисс и прямыми х=а, х=b (рис. 18).
Напомним определение интеграла Римана от f(x) , записываемого в виде (3.1).
Пусть вещественная функция f(x) определена и ограничена на отрезке [a, b].
Разобьем [a, b] на n частичных отрезков [xi, xi+1], 0 ≤ i ≤ n-1, xn=b, x0=a.
На каждом из этих отрезков выберем произвольную точку ξi ( x i ≤ ξ i ≤ x i +1 )
и составим интегральную сумму:
n
S = s1 + s2 +…+ sn =
∑
f(ξi)⋅(xi+1 – xi)
(3.2)
.
i =1
Если существует предел S при стремлении длины наибольшего частичного
отрезка к нулю и произвольных ξi , то этот предел называется интегралом
Римана от f(x)
(3.3)
I =
lim
S .
max x i + 1 − x i → 0
Теорема 3. О существовании определенного интеграла. Если функция f(x)
непрерывна на [a, b], то предел интегральной суммы существует и не зависит
ни от способа разбиения отрезка [a, b] на элементарные отрезки, ни от выбора
точек ξi.
y
M1
M2
Mi
M n-1
n −1
S 2 = ∑ M i ( xi +1 − xi ); M i = max f ( x).
n −1
S = ∑ q i f (ξ i ) .
i =0
Запишем интеграл (3.1) в виде
n−1
I = ∑ qi f (ξ i ) + R.
s2
si
sn-1
Формула (3.4) называется квадратурной, R – погрешность квадратурной формулы.
В случаях, когда подынтегральная функция f(x) не допускает
непосредственного интегрирования, т. е. первообразную нельзя выразить в
элементарных функциях или при табличном способе задания функции
используют методы численного интегрирования. Сущность большинства этих
методов состоит в замене подынтегральной функции f(x) аппроксимирующей
функцией ϕ(x), для которой можно легко записать первообразную в
элементарных функциях, т.е.
∫
y = f(x)
sn
Δxi
ξ1
0
а=x0
ξ2
x1
ξn
ξi
x2
x i-1
xi
xn-2
x n-1
b=xn
Рис. 18. Геометрическая интерпретация определенного интеграла.
33
x
(3.4)
i =0
a
s1
xi ≤ x ≤ xi +1
i =0
Разнообразные формулы численного интегрирования отличаются от (3.2)
только явным указанием способа:
1) выбора xi, ξi;
2) ускорения сходимости в (3.3);
3) оценки погрешности, использующей дополнительную информацию о
поведении f(x).
Обобщим понятие интегральной суммы (3.2). Точки ξi назовем узлами, а
коэффициенты (xi+1 – xi) заменим некоторыми числами qi, не зависящими от f(x),
называемыми весами. Тогда формула (3.2) заменится следующей:
b
Mn
xi ≤ x ≤ xi +1
i =0
b
f ( x)dx = ∫ ϕ ( x)dx + R = S + R ,
a
где S – приближенное значение интеграла, R – погрешность вычисления
интеграла.
Используемые на практике методы численного интегрирования можно
сгруппировать в зависимости от способа аппроксимации подынтегральной
функции.
Методы Ньютона-Котеса основаны на полиномиальной аппроксимации
подынтегральной функции. Методы этого класса отличаются друг от друга
степенью используемого полинома, от которой зависит количество узлов, где
необходимо вычислить функцию f(x).
34
Другими методами, которые могут быть использованы для вычисления
интегралов, является представление подынтегральной функции в виде
степенного ряда (ряда Тейлора). Это позволяет свести вычисление интеграла от
сложной функции к интегрированию многочлена, представляющего первые
несколько членов ряда.
Сплайновые методы базируются на аппроксимации подынтегральной
функции сплайнами, представляющими собой кусочный полином. Методы
различаются по типу выбранных сплайнов.
В методах Монте-Карло узлы выбираются с помощью датчика случайных
чисел. Полученное значение интеграла носит вероятностный характер.
Начало
ввод
n, a, b
h: = (b-a)/n
x: = a
S: = 0
3.1. Методы прямоугольников
Простейшими методами численного интегрирования являются методы
прямоугольников. В этих методах подынтегральная функция на интервале
интегрирования заменяется полиномом нулевой степени или константой.
Подобная замена является неоднозначной, т.к. константу можно выбрать равной
значению подынтегральной функции в любой точке в интервале интегрирования.
В зависимости от этого методы прямоугольников делятся на методы левых,
правых и средних прямоугольников.
Обозначим yi=f(xi), hi=Δxi. Площадь элементарного прямоугольника
рассчитывается по формуле: Si = yi ⋅ hi . Точность метода прямоугольников
пропорциональна h (3.26).
Для метода левых прямоугольников (блок-схема на рис. 19) в качестве точек ξi
выбираются левые границы элементарных отрезков: ξi=xi (i=0,1,…n-1). Для
этого случая получаем следующую формулу:
b
1 0
a
y=S*h
печать y
Конец
При hi = h = const формула примет вид:
∫
S: = S + f(x)
x: = x + h
+ h2 y1 + ... + hn yn −1 .
∫ f ( x)dx ≈ h y
b
i: = 0, n-1
n−1
f ( x)dx ≈ h∑ f ( xi ) .
(3.5а)
Рис.19. Блок-схема метода левых прямоугольников.
i =0
a
Для метода правых прямоугольников (блок-схема на рис. 20) в качестве точек
ξi выбираются правые границы элементарных отрезков (ξi=xi+1), и формула
выглядит следующим образом:
b
∫ f ( x)dx ≈ h y
1 1
+ h2 y 2 + ... + hn y n .
a
При hi=h=const:
b
∫
a
n
f ( x)dx ≈ h∑ f ( xi ) .
(3.5б)
i =1
35
36
Начало
Начало
ввод
n, a, b
ввод
n, a, b
h: = (b-a)/n
x: = a
S: = 0
h: = (b-a)/n
x: = a-h/2
S: = 0
i: = 1, n
i: = 1, n
x: = x + h
S: = S + f(x)
x: = x + h
S: = S + f(x)
y=S*h
y=S*h
печать y
печать y
Конец
Конец
Рис.21. Блок-схема метода средних прямоугольников.
Рис.20. Блок-схема метода правых прямоугольников.
Широко распространенным и более точным (точность метода
пропорциональна h2, см. 3.28) является метод прямоугольников, использующий
значения функции в средних точках элементарных отрезков (в полуцелых узлах)
– метод средних прямоугольников (блок-схема на рис. 21):
b
n −1
a
i =0
∫ f ( x)dx ≈ ∑ h f ( x
i
i
+ hi / 2) .
(3.6)
3.2. Метод трапеций. Вычисление значения интеграла с
заданной точностью
Метод трапеций использует линейную интерполяцию, т.е. график функции
y=f(x) представляется в виде ломаной, соединяющей точки (xi, yi). В этом случае
площадь всей фигуры (криволинейной трапеции) складывается из площадей
элементарных прямолинейных трапеций (рис. 22).
При численном интегрировании с постоянным шагом hi=h=const (i=0,1,…, n-1)
формула принимает соответственно вид:
b
n −1
a
i =0
∫ f ( x)dx ≈ h∑ f ( x
37
i
+ h 2) .
(3.7)
38
f(x)
f(x)
( xi , yi )
Начало
(xi-1 , yi-1 )
N: = TRUNC ((b-a) / SQRT( ε ))
B: = FALSE
yi
yi-1
Процедура IN ( N, I 2 )
I1 : = I 2
N: =2*N
0
xi-1
hi
xi
x
Процедура IN ( N, I 2 )
Рис. 22
Площадь каждой такой трапеции равна произведению полусуммы оснований
на высоту:
y + y i +1
si = i
hi при i=0,1,…,n-1.
2
Складывая все эти равенства, получаем формулу трапеций для численного
интегрирования:
b
1 n −1
(3.8)
hi ( yi + yi +1 ) .
∫a f ( x)dx ≈ 2 ∑
i =0
При hi = h = const формула примет вид:
b
n −1
h
(3.9)
yi ) .
∫a f ( x)dx ≈ 2 ( y 0 + y n + 2∑
i =1
Точность вычисления интеграла по методу трапеций пропорциональна h2 (3.29).
Чтобы найти значение приближенной величины, например, значение
интеграла с заданной точностью, нужно вычислить интеграл I1 с шагом h1, затем
вычисляется интеграл I2 с шагом h1 = h1/2. Если ⏐I2 ⏐≥ 1, то близость
приближённых значений двух величин устанавливается по критерию:
⏐( I2-I1 )/ I2⏐≤ ε .
(3.10),
При⏐I2 ⏐< 1 близость двух приближенных величин устанавливается по критерию
⏐I2-I1⏐≤ ε.
(3.11)
Если критерий не выполняется, то вычисляем интеграл I3 с шагом h3 = h2/2 и
снова сравниваем по критериям близости значения I2 и I3.
Деление шага пополам эквивалентно двукратному увеличению числа
разбиения n отрезка интегрирования. Блок-схема программы для вычисления
интеграла методом трапеций с заданной точностью представлена на рис.23-24.
39
нет
нет
⎜I2⎜≥ 1
да
⎜I2-I 1⎜≤ ε
да
⎜(I2 -I1 )/I 2 ⎜≤ ε
нет
да
B: = TRUE
B: = TRUE
B = TRUE
нет
да
печать
I2 , N
Конец
Рис.23. Блок-схема программы для вычисления интеграла методом трапеций
с заданной точностью
40
y
Начало
Mi+1
y =ϕi(x)
h: = (b-a)/n
S: = 0
x:=a
Mi-1
Mi
i: = 1, n-1
si
x:=x+h
S: = S + f(x)
2h
0
S: = (2⋅S+f(a)+f(b))⋅0.5⋅h
Конец
Рис.24.Блок-схема процедуры IN(N:integer; var S: real);
для вычисления интеграла методом трапеций.
3.3. Метод Симпсона
Разобьем отрезок интегрирования [a, b] на четное число n равных частей с
шагом h. На каждом отрезке [x0,x2], [x2,x4],…, [xi-1,xi+1],…, [xn-2,xn]
подынтегральную функцию f(x) заменим интерполяционным многочленом
второй степени:
f ( x ) ≈ ϕ i ( x ) = a i x 2 + b i x + c i при x i −1 ≤ x ≤ x i +1 .
Коэффициенты этих квадратных трехчленов могут быть найдены из условий
равенства многочлена в точках xi соответствующим табличным данным yi=f(xi).
В качестве ϕi(x) можно принять интерполяционный многочлен Лагранжа второй
степени, проходящей через точки Mi-1(xi-1,yi-1), Mi(xi,yi), Mi+1(xi+1,yi+1) (рис. 24):
(x − xi )(x − xi+1 )
(x − xi−1 )(x − xi+1 )
(x − xi−1 )(x − xi )
yi−1 +
yi +
yi+1 .
(xi−1 − xi )(xi−1 − xi+1 )
(xi − xi−1 )(xi − xi+1 )
(xi+1 − xi−1 )(xi+1 − xi )
Элементарная площадь si может быть вычислена с помощью определенного
интеграла. Учитывая равенства xi+1-xi = xi-xi-1 = h, получаем:
ϕi (x) =
xi +1
si = ∫ ϕ i ( x)dx ≈
xi −1
1
2h 2
xi +1
∫ [( x − x )( x − x
i
i +1
) yi −1 − 2( x − xi −1 )( x − xi +1 ) yi + ( x − xi −1 )( x − xi ) yi +1 ] dx =
xi −1
h
= ( yi −1 + 4 yi + yi +1 ).
3
41
xi-1
xi
xi+1
x
Рис. 25
Проведя такие вычисления для каждого элементарного отрезка [xi-1,xi+1],
просуммируем полученные выражения:
S =
h
( y 0 + 4 y1 + 2 y 2 + 4 y 3 + 2 y 4 + ... + 2 y n − 2 + 4 y n −1 + y n ).
3
Данное выражение для S принимается в качестве значения определенного
интеграла:
b
h
∫ f ( x)dx ≈ 3 [y
0
+ 4( y1 + y3 + ... + y n−1 ) + 2( y 2 + y 4 + ... + y n− 2 ) + y n ] =
a
=
n
h ⎧
⎫
⋅ ⎨ y 0 + y n + ∑ (c + 3) ⋅ y i ⎬ ,
3 ⎩
i =1
⎭
(3.10)
-1 при i нечетном
1 при i четном
Полученное соотношение называется формулой Симпсона.
Эту формулу можно получить и другими способами, например
комбинированием формул прямоугольников и трапеций или двукратным
применением метода трапеций при разбиении отрезка [a, b] на части с шагами h
и 2h. При этом важно добиться, чтобы главные части погрешностей этих методов
при сложении уничтожались.
Число узлов в методе Симпсона обязательно должно быть четным, так как
каждый третий узел используется два раза – он служит окончанием одной
параболы и началом следующей. Точность метода Симпсона пропорциональна h4
(3.30).
Блок-схема программы для вычисления интеграла методом Cимпсона
представлена на рис.26.
где
с=
42
3.4. Метод Монте-Карло
Начало
ввод
a, b, ε
n : = TRUNC ((b-a)/SQRT(SQRT( ε)))
n должно быть
четным
Этот метод, называемый еще методом статистических испытаний, появился в
начале 60-х годов и получил свое название от города Монте-Карло, на весь мир
знаменитого своими игорными домами. Рулетка, как датчик случайных чисел,
явилась
прототипом для
создания
генератора случайных чисел,
необходимых для работы по этому методу.
Идея принципа достаточно проста. Допустим, нужно определить площадь
внутри заданного сложного замкнутого контура (рис.27).
нет
n mod 2 < > 0
да
n:=n+1
h : = (b-a)/n
x:=a
Рис. 27
c:=1
S:=0
с = 1 при i– нечётном
-1 при i – чётном
i : = 1, n -1
S : = h * (S + f(a) + f(b)) / 3
x:=x+h
S: = S + (c + 3) * f(x)
вывод
S
c : = -c
Конец
Рис.26. Блок-схема метода Cимпсона.
43
Известно много способов определения площади, есть различного типа
приборы (планиметры, интеграторы) и т.д., если известно функциональное
уравнение фигуры, можно его проинтегрировать и т. д. Но также можно
воспользоваться методом Монте-Карло.
Рисунок представляется мишенью. Проводится серия N независимых
испытаний – выстрелов, причем плотность вероятности выстрела в
прямоугольник ABCD имеет равномерное распределение.
Тогда при большом количестве выстрелов N в ABCD отношение числа
выстрелов, попавших в измеряемый контур, к общему числу выстрелов
будет равно доле площади измеряемой фигуры от ABCD.
S
N попаданий
≈ изм . фигуры .
N выстрелов
S ABCD
Чем больше N (N≥11), тем точнее производится расчет площади измеряемой
фигуры.
Следующим логическим шагом в создании работоспособного алгоритма
метода Монте-Карло является разбиение прямоугольника ABCD сеткой на
квадраты. Если мы пронумеруем все деления, то каждая клетка будет иметь
определенный адрес (номер строки, номер столбца), то есть фактически мы
получаем матрицу (рис.28).
44
B 1 2 3 …….
Вычисление определенных интегралов по методу Монте-Карло
NС
Рассмотрим интеграл
1
1
I =
f ( x ) dx .
(3.17)
0
2
1
∫
Функция f(x) должна быть такой, что
3
:
:
f
2
( x ) dx тоже существует.
0
Пусть η - равномерно распределенная на отрезке [0, 1] случайная величина,
то есть ее плотность распределения задается соотношением:
M
А
D
Рис. 28
P (η ) =
Следовательно, задача статистических испытаний сводится к задаче
получения набора случайных чисел, распределенных равномерно на наборе
[1…M, 1…N]. Задается случайная величина ξ, математическое ожидание которой
равно искомой площади фигуры z, т. е.
z=М(ξ).
(3.11)
Осуществляется серия n независимых испытаний, в результате которых
получается (генерируется) последовательность n случайных чисел ξ1, ξ2,…, ξn,
по совокупности этих значений приближенно определяется искомая величина
ξ + ξ 2 + ... + ξ n
(3.12)
ξ = 1
≈ z .
1 n
1 n
n∗z
M [ξ n ] = M [ ∑ ξ i ] = ∑ M [ξ i ] =
=z.
n i =1
n i =1
n
то есть неравенство z − ξ < 3 σ
(3.15)
0,997.
45
.
0
0
(3.18)
n
таким образом (рис.29), I ≈ ξ n = 1 ∑ f (ηi ) = f ( x ) ,
n i =1
где ηi - независимые реализации случайной величины η.
(3.19)
f(x)
0
1
x
Рис. 29.
Поскольку наряду с (3.17) по предположению существует
1
∫f
2
( x)dx , то
0
σ 2 = D [ f (η )] = M [ f 2 (η )] − ( M [ f (η )]) 2 .
На
n
σ =
1
f(x)
практике, если σ неизвестна, но во всяком случае, конечна, то ее оценивают
при n ≥11 по формуле:
1 n
∑ (ξ i − ξ ) 2
n − 1 i =1
1
(3.14)
) ≈ 0 . 997 ,
n
выполняется с вероятностью
η ∉ [0,1]
y
Причем в силу центральной предельной теоремы распределение случайной
величины ξ асимптотически нормально. Поэтому при достаточно большом n
(практически при n ≥11) согласно (3.13), (3.14) и известному правилу "трех
сигм" имеем:
р( z − ξ < 3
0,
(3.13)
n ∗σ 2 σ 2
1 n
.
D
ξ
[
]
=
=
∑ i
n 2 i =1
n2
n
σ
0 ≤η ≤1
M [ξ ] = ∫ f ( x) Pη ( x)dx = ∫ f ( x)dx = I ,
Если дисперсия конечна D[ξ] = σ 2, то
D [ξ n ] =
1,
ξ=f(η) также будет случайной
Тогда кусочно-непрерывная функция
величиной, и ее математическое ожидание равно
n
Следовательно,
∫
Так как
D x = M [ X − M [ x ]] 2 = M [ x 2 − 2 x ⋅ M [ x ] + ( M [ x ]) 2 ] =
= M ( x 2 ) − 2 ⋅ ( M [ x ]) 2 + ( M [ x ]) 2 = M ( x 2 ) − ( M [ x ]) 2
(3.16)
46
следовательно
σ
1
2
=
∫
0
1
f 2 ( x ) dx − ( ∫ f ( x ) dx ) 2 < ∞ .
(3.20)
0
Величину σ оценивают либо по формуле (3.16), либо исходя из формулы (3.20),
1
оценивая
∫
f 2 ( x ) dx
сверху, а
0
1
( ∫ f ( x ) dx ) 2
снизу. Погрешность
I −ξn
0
оценивается по формуле (3.15) с вероятностью 0.997.
Случайные величины, равномерно распределенные на отрезке [0, 1] в
персональных компьютерах задаются с помощью физических датчиков или
программ. При применении эти числа называются псевдослучайными, так как
практически они обладают статическими характеристиками, но, строго говоря,
не являются случайными.
Если нужно вычислить интеграл с точностью ε, то есть, чтобы выполнялось
условие
b
∫
f ( x)dx − ξ ≤ ε , то из формулы (3.15) следует, что
3
n
][
Рассмотрим
I=
∫ f ( x ) dx = f ⋅ ( d − a ) ,
то есть площадь криволинейной
a
трапеции, ограниченной графиком функции, отрезком оси х, прямыми х = a и
x = d, равна площади прямоугольника abcd, ордината которого равна среднему
значению функции f (x) на отрезке [а, d] , а другая – длине отрезка (рис.30).
Функция f ( x ) = 1
n
n
∑
i =1
f (η i ) , где ηi - значения случайной величины η.
y
f(x)
∫∫∫
f ( x , y , z ) dxdydz = ( x b − x a )( y b − y a )( z b − z a )
xa y a z a
1
n
n
∑
i =1
f ( xi , y i , z i ),
то выбирают точки со случайными координатами внутри прямоугольного
параллелепипеда.
Преимущества применения метода Монте-Карло
=ε .
2
Таким образом, число независимых испытаний равно n = ⎤ (9σ ) ⎡ , где за
⎥ 2 ⎢
⎦ ε ⎣
обозначена целая часть числа.
При вычислении кратных интегралов метод Монте-Карло часто дает
лучшие результаты, чем другие численные методы.
d
xb y b z b
1.
a
σ
Если имеется функция двух независимых переменных и необходимо
вычислить двойной интеграл по Y области в плоскости ху, то выбирают точки со
случайными координатами внутри области интегрирования, вычисляются
значения функции в этих точках, определяют среднее значение, которое
умножают на площадь области интегрирования.
Если следует вычислить
2.
3.
В многомерном случае при применении метода Монте-Карло число
вычисляемых значений подынтегральной функции растет значительно
медленнее относительно 1/ ε (ε - точность вычисления интеграла), чем
в квадратичных формулах.
Точность метода Монте-Карло не зависит от гладкости
подынтегральной функции.
Простая приспособляемость к форме области интегрирования:
∫∫
Ω
f ( P ) dP =
∫∫
R2
n
f ∗ ( P ) dP =
R
n
∑
i =1
f (ξ i , x i ) ,
где R- сторона квадрата, содержащего заданную область Ω (рис.31), а
f*(P) =
f ( P ),
0,
P∈Ω
P∉Ω
.
R
Ω
c
b
f(x)
Рис. 31.
а
Рис. 30
47
d
x
Блок-схема программы для вычисления интеграла методом Монте-Карло
представлена на рис.32-33.
Оценка погрешности приближенного значения интеграла имеет вид [1]:
48
S n ( f ) − I ( f ) ≤ const ⋅ D( f ) / n ,
Н ачало
где I(f) - точное значение интеграла, Sn(f) – приближенное значение интеграла,
D(f) – дисперсия подынтегральной функции. Из формулы видно, что точность
оценки интеграла обратно пропорциональна квадратному корню из числа
испытаний n. Это обстоятельство обуславливает сравнительно медленную
сходимость метода Монте-Карло: например, чтобы уменьшить погрешность
результата в 10 раз, число испытаний необходимо увеличить в 100 раз.
I: = (n 0 -1)*I
j: = n 0 , n
x: = random (b-a)+a
Начало
R [j]: = f(x)
I: = I + R [j]
Ввод ε
n: = 11
I: = 0
n 0: = 1
I: = I/n
D : =0
расчет
дисперсии
Процедура
INTEG
j: = 1, n
да
(3*sig)/√n ≤ ε
К онец
нет
n 0: = n
n: =TRUNC((9*sig 2 )/ε 2 )
Рис.33. Блок-схема процед уры IN T E G (n 0 , n: integer; var sig, I: real);
печать
n, I*(b-a)
Процедура
INTEG
метода М онте-К арло.
Используя выражение (3.21), в результате вычислений находим
n
1
1
1
3
4
I ≈
(a h + b h 2 + c h + d h ) .
Конец
∑
i =1
Рис.32. Блок-схема основной программы метода Монте-Карло.
3.5. Использование сплайнов для численного интегрирования
Одним из методов численного интегрирования, особенно эффективным при
строго ограниченном числе узлов, является метод сплайнов, использующий
интерполяцию сплайнами (часть II, глава 1.4).
Разобьем отрезок интегрирования [a, b] на n частей точками xi. Пусть xi xi-1 = hi (i=1,2,…, n). На каждом элементарном отрезке интерполируем
подынтегральную функцию f(x) с помощью кубического сплайна:
ϕ i ( x ) = a i + bi ( x − x i −1 ) + c i ( x − x i −1 ) 2 + d i ( x − x i −1 ) 3 ,
i=1,2,…,n.
Выражение для интеграла представим в виде:
b
I =
∫
a
f ( x ) dx =
n
xi
∑ ∫
i =1 xi −1
(3.21)
n
f ( x ) dx ≈ ∑
49
sig: = SQ R T (D /(n-1))
D : = (R [j]-I) 2 + D
i
2
i
i
i
( x ) dx .
3
i
i
4
i
i
(3.22)
Коэффициенты ai, bi, ci, di находятся по формулам (1.27)-(1.29) (часть II,
глава 1), учитывая, что ai = yi-1. Для практических расчетов формула (3.22)
может быть представлена в виде:
1 n
1 n 3
(3.23)
I ≈ ∑ h i ( y i −1 + y i ) −
∑ h i ( c i −1 + c i ) .
2 i =1
12 i =1
Анализ этой формулы показывает, что первый член в правой части совпадает
с правой частью формулы (3.8) для метода трапеций. Следовательно, второй
член характеризует поправку к методу трапеций, которую дает метод сплайнов.
Как следует из формулы (3.21) коэффициенты ci выражаются через вторые
производные ϕ i ″ ( x ) :
1
1 ″
″
(ϕ i ( x i − 1 ) ≈ y i .
2
2
Это позволяет оценить второй член правой формулы (3.23):
ci =
xi
∫ϕ
i
i =1 xi −1
50
hi3
h3
( ci −1 + c i ) ≈ i y ∗″ ,
12
12
Начало
где y ∗′′ -вторая производная в некоторой внутренней точке. Полученная оценка
показывает, что добавка к формуле трапеций, которую дает использование
сплайнов, компенсирует погрешность самой формулы трапеций.
Во всех предыдущих методах формулы численного интегрирования можно
условно записать в виде линейной комбинации табличных значений функции:
b
∫
n
ввод
a, b, n
h: = (b-a)/n
x: = a - h
i : = 0, n
f ( x ) dx = ∑ α i y i .
вычисление
подынтегральной
функции в узлах
i =0
a
При использовании сплайнов такое представление невозможно, поскольку
сами коэффициенты αi зависят от всех значений yi.
Блок-схема для интегрирования методом сплайнов представлена на рис. 34. В
программе следует записать формулы для расчета подынтегральной функции в
узлах интерполяции в виде подпрограммы FUNCTION.
При вычислении интеграла для накопления суммы можно было использовать
формулу (3.23):
x: = x + h
y[i]: = f(x)
i : = 1, n
a[i]: = y[i-1]
3
In : = In + 0,5*h*( y[i-1] + y[i] – h (c[i-1] + c0))/12 .
В этом случае следует опустить в программе расчет коэффициентов сплайна
bi, di.
3.6. Погрешность численного интегрирования
В общем случае погрешность Rn численного значения Sn равна:
b
R n = ∫ f ( x )dx − S n .
In: = 0
a
Она зависит от шага разбиения, и её можно представить в виде Rn = O(hk). В
случае переменного шага можно принять h = max hi .
1≤ i ≤ n
Из этого представления погрешности численного интегрирования следует,
что при h → 0 (n → ∞) значение интеграла, получаемое путем численного
интегрирования, сходится к его точному значению. Заметим, что это имеет
место, если подынтегральная функция на конечном отрезке [a, b] достаточно
гладкая.
Оценим погрешность интегрирования для метода левых прямоугольников.
Для этого разложим подынтегральную функцию f(x) в ряд Тейлора около
«левой» точки хk-1 :
f ( x) = f ( xk −1 ) + ( x − xk −1 ) f ′( xk −1 ) +
51
( x − xk −1 )
f ′′( xk −1 ) + ...
2!
Расчет коэф-тов
трехдиагональной матрицы и
коэф-тов сплайна
c[i], d[i], b[i]
i: = 1,.., n (см. часть II рис. 7)
вычисление
интеграла
i : = 1, n
In : = In + а[i-1] *h + 0,5*b[i] *h 2 + 1/3*c[i] *h 3 + 1/4*d[i]*h 4
печать In
Конец
Рис.34. Блок-схема для интегрирования методом сплайнов.
2
(3.24)
52
Для метода Симпсона:
Интегрируя данное разложение почленно на интервале [xk-1; xk-1+h] получим:
I точн = f ( xk −1 )h +
= f ( xk −1 )h +
( x − xk −1 ) 2
2
xk −1 + h
xk −1
⋅ f ′( xk −1 ) + ... =
h2
f ′( xk −1 ) + ... = S k + Rk ,
2
(3.25)
где Sk – интегральная сумма, Rk – погрешность вычисления интеграла на
интервале [xk-1; xk-1+h].
При малой величине шага h основной вклад в погрешность будет вносить
2
слагаемое h f ′( xk −1 ) , которое называется главным членом погрешности Rок
2
вычисления интеграла на интервале [xk-1, xk-1+h]:
Главный член полной погрешности для интеграла на всем интервале [x0, xn]
определится путем суммирования погрешностей на каждом частичном интервале
[xk-1, xk-1+h]:
h n
hb
Ro = ∑ Rok = ∑ h ⋅ f ′( xk −1 ) = ∫ f ′( x)dx ,
2 k =1
2a
k =1
b−a
Ro =
h ⋅ f ′(ξ ) , где a ≤ ξ ≤ b .
2
n
(3.26)
Данная формула представляет собой теоретическую оценку погрешности
вычисления интеграла методом левых прямоугольников. Эта оценка является
априорной, т.к. не требует знания значения вычисляемого интеграла. Эта оценка
не удобна для практического вычисления погрешности, но полезна для
установления структуры главного члена погрешности. Степень при шаге h в
формуле главного члена погрешности, называется порядком точности метода
интегрирования. Чем больше степень при h, тем точнее метод. Таким образом,
метод левых прямоугольников имеет первый порядок точности.
Далее без вывода приводятся формулы главных членов погрешности для
методов, рассматриваемых в данной главе.
Для метода правых прямоугольников:
Ro = −
b−a
h ⋅ f ′(ξ ) , где a ≤ ξ ≤ b .
2
(3.27)
b−a 2
h ⋅ f ′′(ξ ) , где a ≤ ξ ≤ b .
24
Для метода трапеций:
Ro = −
b−a 2
h ⋅ f ′′(ξ ) , где a ≤ ξ ≤ b .
12
53
b − a 4 IV
h ⋅ f (ξ ) , где a ≤ ξ ≤ b .
180
(3.28)
(3.29)
(3.30)
Анализируя данные формулы, можно сделать вывод, что самым точным из
этих методов является метод Симпсона, а точность метода трапеций в два раза
меньше точности метода средних прямоугольников.
Апостериорные оценки погрешностей по методу Рунге и
процедуре Эйткена
Априорные оценки главного члена погрешности методов интегрирования
можно записать в виде:
(3.31)
Ro = A ⋅ h p ,
где А – коэффициент, зависящий от метода интегрирования и вида
подынтегральной функции; h – шаг интегрирования; р - порядок точности
метода.
Пусть вычисляется значение некоторой переменной I с шагом h, тогда:
(3.32)
I = I h + A ⋅ h p + O(h p +1 ) ,
где Ih – приближенное значение I; A⋅hp – главный член погрешности; O(h p +1 ) пренебрежимо малая величина.
Вычислим ту же самую величину I с шагом kh:
(3.33)
I = I kh + A ⋅ (kh) p + O ((kh) p +1 ) ,
где коэффициент пропорциональности k может быть как больше, так и меньше 1.
Коэффициент А в выражениях (3.32) и (3.33) будет одинаковым, т.к. вычисляется
одна и та же величина одним и тем же методом, а от величины шага h значение А
не зависит.
Пренебрегая бесконечно малыми величинами, приравняем правые части
соотношений (3.32) и (3.33) и с учетом формулы (3.31) получим:
I h + Ro = I kh + k p Ro ,
откуда найдем Ro:
Ro =
Для метода средних прямоугольников:
Ro =
Ro = −
I h − I kh
.
k p −1
(3.34)
Формула (3.34) называется первой формулой Рунге и позволяет путем
двойного пересчета величины I с шагами h и kh оценить погрешность. Т.к.
оценка осуществляется после вычисления I, то она является апостериорной.
После определения Ro можно вычислить уточненное значение искомой
величины:
54
I уточн = I h + Ro .
Последнее соотношение называется второй формулой Рунге. К сожалению,
погрешность уточненного значения остается неопределенной, хотя, как правило,
она меньше значения Ro.
Процедура Эйткена дает возможность оценить погрешность метода O(hp) и
указывает алгоритм уточнения результатов. Расчет проводится последовательно
три раза при различных шагах разбиения h1, h2, h3, причем их отношения
постоянны: h2/h1 = h3 / h2 = q (например, при делении шага пополам q = 0.5).
Пусть в результате численного интегрирования получены значения интеграла I1,
I2, I3. Тогда уточненное значение интеграла вычисляется по формуле:
(I1 − I 2 ) 2 ,
I=I −
1
а
порядок
I1 − 2I 2 + I 3
точности
используемого метода численного интегрирования
определяется соотношением p = 1 ln I 3 − I 2 .
ln q I 2 − I 1
Задание. Разработайте алгоритмы вычисления интеграла по методам
прямоугольников, трапеций, Симпсона, в котором приближенное значение
интеграла уточняется с помощью метода Рунге и процедуры Эйткена.
3.7. Контрольные вопросы
1. Геометрический смысл определенного интеграла.
2. Какой зависимостью связан шаг интегрирования с количеством
интервалов?
3. Какой из рассматриваемых методов является самым точным, и как это
определяется?
4. От чего зависит точность получаемого результата интегрирования?
5. Возможно ли получение точного значения результата методом трапеций
для линейной подынтегральной функции?
6. Основной член погрешности методов интегрирования.
7. Почему для метода Симпсона число интервалов должно быть четным?
8. Что такое апостериорная оценка погрешности результата?
9. Может ли значение интеграла получиться отрицательным числом?
10. Чему равен шаг при вычислении интеграла с заданной точностью?
11. Что дает процедура Эйткена?
55
4. Методы решения систем линейных алгебраических
уравнений
Методы решения систем линейных уравнений делятся на две группы –
прямые и итерационные. Прямые методы используют конечные соотношения
(формулы) для вычисления неизвестных. Они дают решение после выполнения
заранее известного числа операций. Эти методы сравнительно просты и наиболее
универсальны, т.е. пригодны для решения широкого класса систем линейных
уравнений.
Вместе с тем прямые методы имеют и ряд недостатков. Как правило они
требуют хранения в оперативной памяти ЭВМ сразу всей матрицы, и при
больших значениях n расходуется много места в памяти. Далее, прямые методы
обычно не учитывают структуру матрицы – при большом числе нулевых
элементов в разряжённых матрицах (например, клеточных или ленточных) эти
элементы занимают место в памяти машины, и над ними проводятся
арифметические действия. Существенным недостатком прямых методов является
также накапливание погрешностей в процессе решения, поскольку вычисления
на любом этапе используют результаты предыдущих операций. Это особенно
опасно для больших систем, когда резко возрастает общее число операций, а
также для плохо обусловленных систем, весьма чувствительных к погрешностям.
В связи с этим прямые методы используются обычно для сравнения небольших
(n < 200) систем с плотно заполненной матрицей и не близким к нулю
определителем.
Отметим ещё, что прямые методы решения линейных систем иногда
называют точными, поскольку решение выражается в виде точных формул через
коэффициенты системы. Однако точное решение может быть получено лишь при
выполнении вычислений с бесконечным числом разрядов (разумеется, при
точных значениях коэффициентов системы). На практике при использовании
ЭВМ вычисления проводятся с ограниченным числом знаков, определяемым
разрядностью машины. Поэтому неизбежны погрешности в окончательных
результатах за счёт округления.
Итерационные методы – это методы последовательных приближений. В них
необходимо задать некоторое приближённое решение – начальное приближение.
После этого с помощью некоторого алгоритма проводится один цикл
вычислений, называемый итерацией. В результате итерации находят новое
приближение. Итерации проводятся до получения решения с требуемой
точностью. Алгоритмы решения линейных систем с использованием
итерационных методов обычно более сложные по сравнению с прямыми
методами. Объём вычислений заранее определить трудно.
Тем не менее итерационные методы в ряде случаев предпочтительнее.
Погрешности окончательных результатов не накапливаются, поскольку точность
вычислений в каждой итерации определяется лишь результатами предыдущей
56
итерации и практически не зависит от ранее выполненных вычислений. Таким
образом, итерационные методы являются особенно полезными в случае
большого числа уравнений, а также плохо обусловленных систем. Следует
отметить, что при этом сходимость итераций может быть очень медленной;
поэтому ищутся эффективные пути её ускорения.
Итерационные методы могут использоваться для уточнения решений,
полученных с помощью прямых методов. Такие смешанные алгоритмы обычно
довольно эффективны, особенно для плохо обусловленных систем.
Наиболее распространёнными среди прямых методов являются метод
исключения Гаусса и его модификации.
4.1. Метод Гаусса
Пусть линейное алгебраическое уравнение имеет вид:
А⋅x=b,
(4.1)
где А – вещественная квадратная матрица порядка n, b – заданный вектор, х –
искомый вектор. Если матрица А неособенная, то есть определитель матрицы
отличен от нуля, тогда для каждого вектора b система (4.1) имеет единственное
решение.
Запишем систему (4.1) в развёрнутом виде:
a11x1+a12x2+…+a1nxn=b1
a21x1+a22x2+…+a1nxn=b2
……………………………
(4.2)
an1x1+an2x2+…+annxn=bn
Метод Гаусса основан на приведении матрицы к треугольному виду, что
достигается последовательным исключением неизвестных х1, х2, …, хn из данной
системы. Метод состоит из двух этапов – прямого и обратного хода.
Прямой ход. Предположим, что а11 не равно нулю. Поделив первое уравнение
на а11 получим:
x1 + c12x2 +…+ c1nxn = y1 ,
(4.3)
где
bi(1) = b j − y1 ai1 , i,j=2,3,…,n.
Матрица системы (4.5) имеет вид:
⎡ 1 c12 K c1n
⎢ 0 c (1) K c (1)
22
2n
⎢
⎢ KKKKKK
⎢
(1)
(1)
⎣ 0 cn 2 K cnn
a1 j
ai1x1 + ai2x 2+…+ ainxn = bi,
i=2,3,…,n.
(4.4)
Умножим (4.3) на аi1 и вычтем полученное уравнение из i-го уравнения
системы (4.4). В результате получим следующую систему уравнений:
x1 + c12 x2 + K + c1 j x j + c1n xn = y1 ,
(1)
c22
x2 + K + c2(1j) x j + K + c2(1n) xn = y2(1) ,
……………………………………….
(1)
nj
(1)
nn n
57
(4.7)
Исключая таким же образом неизвестные х2, х3,…,хn, придём окончательно к
системе уравнений вида:
x1+c12x2+…+c1nxn = y1,
x2+…+c2nxn = y2,
………………………
(4.8)
xn-1+…+cn-1xn = yn-1,
xn = yn,
эквивалентной исходной системе (4.2).
Матрица этой системы имеет треугольный вид. На этом заканчивается
прямой ход метода Гаусса.
Заметим, что в процессе исключения неизвестных приходится выполнять
операции деления на коэффициенты а11, а22 и т. д. Поэтому они должны быть
отличными от нуля; в противном случае необходимо соответственным образом
переставить уравнения системы. Перестановка уравнений должна быть
предусмотрена в вычислительном алгоритме при его реализации на ЭВМ.
Обратный ход. Он заключается в нахождении неизвестных х1, х2, …, хn из
системы (4.8). Поскольку матрица системы имеет треугольный вид, можно
последовательно, начиная с хn, найти все неизвестные. Общие формулы
обратного хода имеют вид:
n
xi = yi − ∑cij x j ,
i = n-1,…,1,
xn = yn .
(4.9)
На рис. 35 приведена блок-схема решения методом Гаусса системы n
линейных уравнений вида (4.2). При вводе коэффициентов системы (4.2) должна
быть сформирована расширенная матрица коэффициентов А размерностью
n+1×n, в n+1 столбец которой записываются правые части коэффициентов bi
(i=1,…,n). Левая часть блок-схемы соответствует прямому ходу. Поясним смысл
используемых переменных: ε - малая положительная величина (10-5÷10-7),
0, если все коэффициенты столбца i матрицы А равны нулю
Р=
номеру строки, в которой стоит первый, отличный от нуля
коэффициент в i-ом столбце;
(4.5)
c x +K+ c xj +K+ c x = y .
(1)
n2 2
⎤
⎥
⎥.
⎥
⎥
⎦
(4.6)
j =i+1
b
, j=2,…, n,
y1 = 1 .
a 11
a 11
Рассмотрим теперь оставшиеся уравнения системы (4.2):
c1 j =
Здесь обозначено:
aij(1) = aij − c1 j a i1 ,
(1)
n
58
i – номер исключаемой переменной,
k – номер уравнения, из которого исключается переменная xi,
c – вспомогательная переменная,
s – переменная для накопления суммы из уравнения (4.9).
При прямом ходе матрица А преобразуется в треугольную матрицу, с
единицами на главной диагонали.
Операция перестановки уравнений (т. е. перестановки соответствующих
коэффициентов) служит для предотвращения деления на нулевой элемент.
Правая часть блок-схемы описывает процесс обратного хода.
Метод Гаусса целесообразно использовать для решения систем с плотно
заполненной матрицей. Все элементы матрицы и правые части системы
уравнений находятся в оперативной памяти машины. Объём вычислений
определяется порядком системы n: число арифметических операций примерно
равно (2/3)n3.
Видоизмените алгоритм, чтобы в треугольной матрице на главной диагонали
стояли преобразованные коэффициенты Аii.
4.2. Метод Гаусса с выбором главного элемента
Этот метод позволяет избежать указанных выше трудностей. Основная идея
метода состоит в том, чтобы на очередном шаге исключать не следующее по
номеру неизвестное, а то неизвестное, коэффициент при котором является
наибольшим по модулю. Таким образом, в качестве ведущего элемента здесь
выбирается главный, значение которого наибольшее по модулю. Тем самым,
если определитель матрицы А не равен нулю, то в процессе вычислений не будет
происходить деление на нуль.
Различные варианты метода Гаусса с выбором главного элемента
проиллюстрируем на примере системы из двух уравнений:
⎧a11 x1 + a12 x 2 = b1 .
⎨
⎩ a 21 x1 + a 22 x 2 = b2
Применяется также метод Гаусса с выбором главного элемента по столбцу.
Предположим, что a 21 > a11 . Перепишем систему (4.9) в виде:
⎧a 21 x1 + a 22 x 2 = b2
,
⎨
⎩ a11 x1 + a12 x 2 = b1
и к новой системе применим на первом шаге обычный метод Гаусса. Таким
образом, метод Гаусса с выбором главного элемента по столбцу эквивалентен
применению обычного метода Гаусса к системе, в которой на каждом шаге
исключения поводится соответствующая перенумерация уравнений.
Иногда применяется и метод Гаусса с выбором главного элемента по всей
матрице, когда в качестве ведущего выбирается максимальный по модулю
элемент среди всех элементов матрицы системы.
Блок-схема метода Гаусса с выбором главного элемента в столбце приведена
на рисунке 36. Здесь введены следующие индексы: max – максимальный по
модулю элемент, стоящий в столбце i; k – номер строки, в которой стоит max.
Заметим, что диагональные элементы матрицы называются ведущими
элементами; ведущий элемент aii – это коэффициент при i-ом неизвестном в i-ом
уравнении на i-ом шаге исключения.
Благодаря выбору наибольшего по модулю ведущего элемента уменьшаются
множители, используемые для преобразования уравнений, что способствует
снижению погрешностей вычислений. Поэтому метод Гаусса с выбором
главного элемента обеспечивает приемлемую точность решения для
сравнительно небольшого числа (n ≤ 100) уравнений. И только для плохо
обусловленных систем решения, полученные по этому методу, ненадёжны.
Задание. Модифицируйте алгоритм так, чтобы главные элементы выбирались
по всей матрице.
Предположим, что a12 > a11 . Тогда на первом шаге будем исключать
переменное х2. Такой
переписывается в виде:
приём
эквивалентен
тому,
что
система
(4.2)
⎧a12 x 2 + a11 x1 = b1 ,
(4.10)
⎨
a
x
a
x
b
+
=
21 1
2
⎩ 22 2
и к (4.10) применяется первый шаг обычного метода Гаусса.
Указанный способ исключения называется методом Гаусса с выбором
главного элемента по строке. Он эквивалентен применению обычного метода
Гаусса к системе, в которой на каждом шаге исключения проводится
соответствующая перенумерация переменных.
59
(4.11)
60
Начало
3
ввод N, ε,
расширенной
матрицы А
2
Обратный ход
Прямой ход
X[N] := A[N, N+1]
3
i := 1, N
р:=0
j := i
i := N-1
Поиск ненулевого элемента
в i –ом столбце
нет
Модуль значения элемента
отличен от нуля?
да
да
р := j; j := N
k := i+1, N
S := S+A[i, k]*X[k]
j := j+1
нет
Меняем местами
строки i и р
i<>р
Печать ‘все
нули в
столбце’,i,’нет
единственного
решения’
да
нет
exit
да
j := i, N+1
C := A[р, j]
A[р, j] := A[i, j]
A[i, j] := C
да
S := 0
нет
abs(A[j,i])> ε
j<>р
нет
i≥1
j≤N
Деление строки
i на A[i,i ]
X[i] := A[i, N+1]-S
i := i -1
C := A[i, i]
j := i, N+1
Печать
корней
A[i, j] := A[i, j]/C
k := i+1, N
2
Получение нулей под
главной диагональю
в i-ом столбце
Проверка
правильности
получения корней
C:= -A[k, i]
j := i, N+1
Конец
A[k, j]:=A[k, j]+C*A[i, j]
Рис. 35. Блок схема метода Гаусса
61
62
Начало
ввод N, ε, расширенная матрица А
Прямой ход
2
3
i := 1, N
1
j := i
max := 0
Деление строки
i на A[i,i ]
C := A[i, i]
Поиск максимального по
модулю элемента в i –ом столбце
j≤N
нет
max < ε
да
нет
нет
max < abs(A[j,i])
да
k := j; max:=abs(A[j,i])
да
печать
“Бесчисленное
множество
решений”
Exit
j := j+1
Меняем местами
строки i и k
Получение единицы
на главной
A[i, j] := A[i, j]/C
k := i+1, N
2
C := -A[k, i]
Ведущий элемент не на
главной диагонали?
j <> k
j := i, N+1
j := i, N+1
нет
A[k, j] := A[k, j]+C*A[i, j]
да
1
j := i, N+1
C := A[k, j]
A[k, j] := A[i, j]
A[i, j] := C
63
64
4.3. Метод Гаусса-Зейделя
3
Обратный ход
X[N] := A[N, N+1]
i := N-1
i ≥1
нет
Одним из самых распространённых итерационных методов, отличающийся
простотой и лёгкостью программирования, является метод Гаусса-Зейделя.
Проиллюстрируем сначала этот метод на примере решения системы:
⎧ a11x1 + a12x2 + a13x3 = b1
⎪
.
(4.12)
⎨a21x1 + a22x2 + a23x3 = b2
⎪a x + a x + a x = b
⎩ 31 1 32 2 33 3 3
Предположим, что диагональные элементы а11, а22, а33 отличны от нуля (в
противном случае можно переставить уравнения). Выразим неизвестные х1, х2, х3
соответственно из первого, второго и третьего уравнений системы (4.12):
1
(4.13)
x =
⋅ (b − a x − a x ) ,
1
да
S := S+A[i, k]*X[k]
X[i] := A[i, N+1]-S
1
12
2
13
3
1
⋅ (b2 − a 21 x 2 − a 23 x 3 ) ,
a 22
1
x3 =
⋅ (b3 − a 31 x 2 − a 33 x 3 ) .
a 33
некоторые начальные приближения
(4.14)
x2 =
S := 0
k := i+1, N
a 11
Зададим
(4.15)
значений
неизвестных
x 10 , x 20 , x 30 .
Подставляя эти значения в правую часть уравнения (4.13), найдём новое
(первое) приближение для x1:
1
x11 =
⋅ ( b1 − a 12 x 20 − a 13 x 30 ) .
a 11
Используя это значение для x1 и приближение x 30 для х3, находим из
i := i-1
уравнения (4.14) первое приближение для x2:
1
x 12 =
⋅ ( b 2 − a 21 x 11 − a 23 x 30 ) .
a 22
Печать
корней
И наконец, используя вычисленные значения x1 , x 2 , находим с помощью
выражения (4.15) первое приближение для x3:
1
x 31 =
⋅ (b3 − a 31 x11 − a 32 x 12 ) .
a 33
Первая итерация решения системы (4.12) на этом заканчивается.
Используя найденные значения x 11 , x 12 , x 31 можно таким же способом
Проверка
правильности
получения корней
Конец
Рис. 36. Блок-схема метода Гаусса с выбором
главного элемента в столбце.
65
1
1
провести вторую итерацию, в результате которой будут найдены вторые
приближения к решению x 12 , x 22 , x 32 . Приближение с номером k можно
представить в виде:
66
1
⋅ ( b1 − a12 x 2k −1 − a13 x 3k −1 ) ,
a 11
1
x 2k =
⋅ ( b 2 − a 21 x1k − a 23 x 3k −1 ) ,
a 22
x1k =
x 31 =
Начало
1
⋅ ( b3 − a 31 x1k − a 32 x 2k ) .
a 33
Ввод ε, M, N,
матрицы А,
вектора В
Итерационный процесс продолжается до тех пор, пока значения x 1k , x 2k , x 3k
k −1
k −1
i := 1, n
k −1
не станут близкими с заданной погрешностью ε к значениям x 1 , x 2 , x 3 .
Рассмотрим теперь систему, состоящую из n линейных уравнений с n
неизвестными. Представим её в следующем виде:
ai ,1 x1 + ai , 2 x2 + K + ai ,i xi + ai ,i +1 xi +1 + K + ai ,n xn = bi ,
i=1,2,…,n.
Предположим также, что все диагональные элементы отличны от нуля. Тогда
в соответствии с методом Гаусса-Зейделя k-е приближение к решению можно
представить в виде:
X[i] = 0
k := 0
4
δ := 0
1
x =
⋅ (bi − ai ,1 x1k − K − ai ,i −1 xik−1 − ai ,i +1 xik+−11 − K − ai ,n xnk −1 ) ,
a i ,i
k
i
k := k+1
i=1,2,…,n.
Причём при i=1 отсутствует часть α = a i ,1 x1k − K − a i ,i −1 x ik−1 .
S := 0
Итерационный процесс продолжается до тех пор пока все значения x
станут близкими к x
k −1
i
k
i
не
с заданной точностью. Близость этих значений можно
характеризовать максимальной величиной их разности δ. Тогда при заданной
допустимой погрешности ε > 0 критерий окончания итерационного процесса
можно представить в виде:
δ = max x ik − x ik −1 < ε
1≤ i ≤ n
.
(4.16)
Это критерий по абсолютным отклонениям при x ik ≤ 1 . Можно заменить его
2
i :=1, n
3
j :=1, n
a[i , i] > ε
S := S-a[ i, j]·x[j]
да
1≤ i≤ n
x ik − x ik − 1
<ε .
x ki
xx:= (b[ i]+S)/a[i,i]
(4.17)
d := abs(xx-x[i])
1
67
Печать
‘a’,i, i,’=’,
нулю
Exit
критерием по относительным разностям, т.е. условие окончания итерационного
процесса можно записать в виде (при x ik >> 1 ):
δ = max
нет
68
3
1
δ<ε
да
нет
4
нет
Вывод
корней
x[i];
i=1,n
δ<d
да
δ=d
нет
x[i] := xx
K > M-1
да
2
Нет сходимости
за М итераций
Конец
Рис. 37. Блок-схема метода Гаусса-Зейделя.
Для сходимости итерационного метода достаточно, чтобы модули
диагональных коэффициентов для каждого уравнения системы были не меньше
сумм модулей всех остальных коэффициентов:
,
i=1,2,…,n.
(4.18)
a
≥
a
ii
∑
i≠ j
ij
При этом хотя бы для одного уравнения неравенство (4.18) должно
выполняться строго. Эти условия являются достаточным условием для
сходимости метода, но не являются необходимыми, т.е. для некоторых систем
итерации сходятся и при нарушении этих условий.
В качестве исходных данных вводятся коэффициенты и правые части
уравнений системы, погрешность ε , допустимое число итераций М, а также
начальные приближения переменных хi (i=1,2,…,n). Отметим, что начальные
приближения можно не вводить в ЭВМ, а полагать их равным некоторым
значениям (например, нулю).
Блок-схема данного метода представлена на рис. 37.
4.4. Метод прогонки
Метод прогонки является модификацией метода Гаусса для частного случая
разреженных систем – системы уравнений с трёхдиагональной матрицей. Такие
системы получаются при моделировании некоторых инженерных задач, а также
при численном решении краевых задач для дифференциальных уравнений.
Запишем систему уравнений в виде:
69
= d1
⎧b1x1 + c1x2
⎪
= d2
a2 x1 + b2 x2 + c2 x3
⎪
= d3
a3 x2 + b3 x3 + c3 x4
(4.19)
⎪⎪
⎨
⎪ KKKKKKKKKKKKKKKKKKKKKKK
⎪
an−1xn−2 + bn−1xn−1 + cn−1xn = dn−1
⎪
⎪⎩
an xn−1 + bn xn = dn
На главной диагонали матрицы этой системы стоят элементы b1, b2,…,bn, над
ней – элементы c1, c2,…,cn-1, под ней – элементы a1, a2,…,an. При этом обычно все
коэффициенты bi не равны нулю.
Метод прогонки состоит из двух этапов – прямой прогонки (аналога прямого
хода метода Гаусса) и обратной прогонки (аналога обратного хода метода
Гаусса). Прямая прогонка состоит в том, что каждое неизвестное xi выражается
через xi+1 с помощью прогоночных коэффициентов AAi, BBi:
xi =AAi xi+1+BBi,
i=1,2,…, n-1.
(4.20)
Из первого уравнения системы (4.19) найдём
c
d
x1 = − 1 ⋅ x 2 + 1 .
b1
b1
С другой стороны, по формуле (4.20) x1=AA1x2+BB1. Приравнивая
коэффициенты в обоих выражениях для x1, получаем
c
d
(4.21)
AA 1 = − 1 ,
BB 1 = 1 .
b1
b1
Из второго уравнения системы (4.19) выразим x2 через x3, заменяя x1 по
формуле (4.20):
a2 ⋅ ( AA1 x2 + BB1 ) + b2 x2 + c2 x3 = d 2 .
Отсюда найдём
или
− c 2 x 3 + d 2 − a 2 BB 1
,
a 2 AA 1 + b 2
x2 = AA2 x3 + BB2 ,
d 2 − a2b1
c
, e 2 = a 2 AA1 + b2 .
AA2 = − 2 , BB2 =
e2
e2
x2 =
Аналогично можно вычислить прогоночные коэффициенты для любого
номера i :
AA i = −
ci ,
ei
BBi =
ei = a i AA i −1 + bi ,
70
di − ai bi−1
,
ei
i=2,3,…, n-1.
(4.22)
Обратная прогонка состоит в последовательном вычислении неизвестных xi.
Сначала нужно найти xn. Для этого воспользуемся выражением (4.20) при i=n-1 и
последним уравнением системы (4.19). Запишем их:
xn−1 = AAn−1 xn + BBn−1 ,
Начало
ввод n
a n x n −1 + bn x n = d n .
Отсюда
i := 1, n
xn =
d n − an BBn −1
.
an AAn −1 + bn
ввод массивов
A[i], B[i],
C[i], D[i]
(4.23)
Далее, используя формулы (4.20) и выражения для прогоночных
коэффициентов (4.21), (4.22), последовательно вычисляем все неизвестные
xi для i=n-1, n-2,…,1.
При анализе алгоритма метода прогонки надо учитывать возможность
деления на нуль в формулах (4.22). Можно показать, что при выполнении
условия преобладания диагональных элементов, т.е. если
bi ≥ a i + ci ,
Расчет прогоночных
коэффициентов
AA[1] := -C[1]/B[1]
BB[1] :- D[1]/B[1]
i := 2, n-1
причём хотя бы для одного значения i имеет место строгое неравенство, деления
на нуль не возникает, и система (4.19) имеет единственное решение.
Приведённое условие преобладания диагональных элементов обеспечивает
также устойчивость метода прогонки относительно погрешностей округлений.
Последнее обстоятельство позволяет использовать метод прогонки для решения
больших систем уравнений. Заметим, что данное условие устойчивости прогонки
является достаточным, но не необходимым. В ряде случаев для хорошо
обусловленных систем вида (4.19) метод прогонки оказывается устойчивым даже
при нарушении условия преобладания диагональных элементов. Блок-схема
метода прогонки приведена на рис. 38.
E := A[i]*AA[i-1]+B[i]
AA[i] := -C[i]/E
BB[i] := (D[i]- A[i]*BB[i-1])/E
4.5. Метод Гаусса-Жордана
X[i]: =AA[i]*X[i+1]+BB[i]
Этот метод приводит системы линейных уравнений (4.2) к диагональному
виду.
Для лучшего понимания метода решим сначала систему линейных уравнений
вручную. Пусть даны три уравнения с тремя неизвестными:
4x + 3 y + z = 13
2x − y − z = −3
7 x + y − 3z = 0
71
Расчет
корней
X[n] :=
D[n] - A[n] ⋅ BB[n - 1]
B[n] + A[n] ⋅ AA[n - 1]
i := n-1,1,-1
вывод
корней
Конец
Рис. 38. Блок-схема метода прогонки.
72
Если написать эту систему уравнений, отбросив символы «х», «у», «z» и «=»,
то получится так называемая расширенная матрица системы:
⎛ 4 + 3 + 1 13 ⎞
⎜
⎟
⎜ 2 − 1 − 1 − 3⎟
⎜7 +1 − 3 0 ⎟
⎝
⎠
Решение системы линейных уравнений не изменится, если умножить правую
и левую части каждого уравнения на любое число кроме нуля.
Чтобы первый коэффициент первого уравнения стал равен единице, разделим
первое уравнение на 4:
⎛ 1 + 3 / 4 + 1/ 4 13/ 4⎞
⎜
⎟
−1
−3 ⎟
⎜ 2 −1
⎜ 7 +1
−3
0 ⎟⎠
⎝
[2] − 2 ⋅ [1]
[3] − 7 ⋅ [1]
Поскольку уравнения системы можно умножать на любые числа, а также
складывать и вычитать, коэффициенты, стоящие в первом столбце, кроме
первого можно обратить в нуль. Для этого умножим строку [1] на 2 и
произведение вычтем из второй строки. Первую строку умножим на 7 и
произведение вычтем из третьей строки:
⎛ 1 + 3 / 4 + 1/ 4 13/ 4 ⎞
⎜
⎟
⎜ 0 − 5 / 2 − 3 / 2 − 19 / 2 ⎟ /(−5 / 2)
⎜ 0 − 17 / 4 − 19 / 4 − 91/ 4 ⎟
⎝
⎠
Теперь надо обратить в единицу диагональный элемент второй строки. Для
этого разделим второе уравнение на этот диагональный элемент (-5/2) и получим:
⎛ 1 + 3 / 4 − 1/ 4 13 / 4 ⎞ [1] − 3 / 4 ⋅ [2]
⎟
⎜
1
3/ 5
19 / 5 ⎟
⎜0
⎜ 0 − 17 / 4 − 19 / 4 − 91/ 4 ⎟ [3] + 17 / 4 ⋅ [2]
⎠
⎝
⎛1 0 −1/ 5 2/ 5 ⎞ [1] +1/ 5⋅[3]
⎜
⎟
⎜0 1 3/ 5 19/ 5⎟ [2] − 3/ 5⋅[3]
⎜0 0 1
3 ⎟⎠
⎝
Прибавляя умноженную на соответствующие коэффициенты третью строку к
первой и второй, можно исключить недиагональные элементы третьего столбца,
что соответствует исключению переменной строке х3 из 1 и 2 уравнения:
⎛1 0 0 1⎞
⎜
⎟
⎜0 1 0 2⎟
⎜0 0 1 3⎟
⎝
⎠
Если теперь дописать «забытые» символы, то получится следующая система
уравнений:
1x + 0 y + 0 z = 1
0 x − 1y − 0 z = 2
0 x + 0 y − 1z = 3
Из этой системы уравнений непосредственно получается искомое решение:
х = 1;
у = 2;
z = 3.
Преимущество метода состоит в том, что исключается обратный ход в методе
Гаусса. Его недостатком является увеличение объёма вычислений.
Можно показать, что наибольшая точность достигается тогда, когда ведущий
элемент имеет наибольшее значение. Поэтому строку с нулевым или малым по
абсолютной величине ведущим элементом надо заменить на ту из стоящих под
ней строк, в которой в том же столбце стоит элемент, имеющий наибольшее по
модулю значение.
Блок-схема метода Гаусса-Жордана приведена на рис. 39.
Задание. Преобразуйте алгоритм в метод Гаусса-Жордана с выбором главного
элемента в ведущем столбце.
Теперь попытаемся исключить x2 из первого и третьего уравнения. Для этого
вычтем из первого уравнения умноженное на 3/4 второе уравнение, а третьму
уравнению прибавим второе, умноженное на 17/4:
2/5 ⎞
⎛ 1 0 − 1/ 5
⎜
⎟
19 / 5 ⎟
⎜0 1 3/ 5
⎜ 0 0 − 11/ 5 − 33 / 5⎟ /(−11/ 5)
⎝
⎠
Для того чтобы диагональный элемент третьей строки стал равен единице,
разделим третье уравнение на диагональный элемент третьей строки:
Непосредственное нахождение определителя требует большого объёма
вычислений. Вместе с тем легко вычисляется определитель треугольной
матрицы: он равен произведению её диагональных коэффициентов.
Пусть требуется вычислить определитель треугольной матрицы А1
размерностью n x n. Для приведения матрицы к треугольному виду можно
воспользоваться методом Гаусса или методом Гаусса с выбором главного
элемента.
73
74
4.6. Вычисление определителя по методу Гаусса
Начало
1
ввод N, ε
C := 1/A[i, i]
i := 1, N
j := i, N+1
j := 1, N
A[i, j] := A[i, j]*C
Деление i-ой строки на
диагональный элемент
ввод
коэффициентов
A[i, j]
Получение нулей в
верхнем столбце
j := 1, N
Печать ‘ввод
правой части
уравнения’, i
да
C := -A[j, i]
k := i, N+1
3
i := 1, N
Поиск не нулевого
элемента в i –ом столбце
j := 1, N
да
abs(A[i, j])≤ε
“нет
единственного
решения”
нет
i=j
Меняем местами
строки i и k
i=j
нет
ввод
A[i, N+1]
2
2
да
A[j, k]:=A[j, k]+C*A[i, k]
3
Проверка правильности
полученных корней
i := 1, N
нет
Exit
k := i, N+1
X[i]:=A[i, N+1]
B := A[i, k];
A[i, k] := A[j, k];
A[j, k] := B
Печать
X’(‘,i ,’)= ‘, X[i]
Конец
Рис. 39. Блок-схема метода Гаусса-Жордана.
1
75
76
Знак определителя меняется на противоположный при перестановке его
столбцов или строк. Следовательно, значение определителя после приведения
матрицы А к треугольному виду вычисляется по формуле:
n
det A = u ⋅ ∏ a kk ,
(4.26)
a11ξ1( 0 ) + a12ξ 2( 0 ) + K + a1nξ n( 0 ) = r1( 0) ,
a21ξ1( 0 ) + a22ξ 2( 0 ) + K + a2 nξ n( 0 ) = r2( 0 ) ,
KKKKKKKKKKKKKKK
an1ξ1( 0 ) + an 2ξ 2( 0 ) + K + annξ n( 0 ) = rn( 0 ) .
k =1
где
u=
Решая эту систему, находим значения погрешностей
+1, если число перестановок четное;
-1, если число перестановок нечетное.
аkk – диагональные элементы преобразованной матрицы.
Этот метод позволяет вычислять определители 100-го и большего порядков.
Блок-схема данного метода приведена на рис. 40.
(4.30)
ξ i( 0) ,
которые
используем в качестве поправок к решению. Следующие приближения
неизвестных имеют вид:
x 1(1 ) = x 1( 0 ) + ξ 1( 0 ) , x2(1) = x2( 0 ) + ξ 2( 0 ) , …, x n(1 ) = x n( 0 ) + ξ n( 0 ) .
Таким же способом можно найти новые поправки к решению ξ i(1) и
следующие приближения переменных xi( 2 ) = xi(1) + ξ i(1) и т.д., пока все очередные
4.7. Метод итераций
значения погрешностей (поправок)
Этот метод позволяет уточнить решение, полученное с помощью прямого
метода.
Найдём решение системы линейных уравнений:
a11x1+a12x2+…+a1nxn=b1
a21x1+a22x2+…+a1nxn=b2
……………………………
(4.27)
an1x1+an2x2+…+annxn=bn
Пусть с помощью некоторого прямого метода вычислены приближённые
значения неизвестных x1( 0 ) , x 2( 0 ) , K , x n( 0 ) . Подставляя это решение в левые части
(0)
системы (4.27), получаем векторные значения bi
(0)
11 1
(0)
21 1
a x
a x
+ a12 x
+ a 22 x
(0)
2
(0)
2
+K+ a x
+K+ a x
(0)
1n n
(0)
2n n
, отличные от bi (i=1,2,…,n):
= b1( 0 ) ,
= b2( 0 ) ,
(4.28)
K K K K K K K K K K K KK K K
a n1 x1( 0 ) + a n 2 x 2( 0 ) + K + a nn x n( 0 ) = bn( 0 )
ri ( 0 ) = bi − bi( 0 ) ,
i=1,2,…,n.
(4.29)
Вычитая каждое уравнение системы (4.28) из соответствующего уравнения
системы (4.27), с учётом обозначений (4.29) получаем:
77
не станут достаточно малыми.
Заметим, что для нахождения очередного приближения, т.е. на каждой
итерации, решаются системы уравнений вида (4.30) с одной и той же матрицей,
являющейся матрицей исходной системы (4.27) при разных правых частях. Это
позволяет строить экономичные алгоритмы. Например, при использовании
метода Гаусса сокращается объём вычислений на этапе прямого хода.
Решение системы уравнений с помощью рассмотренного метода, а также при
использовании других итерационных методов сводится к следующему. Вводятся
исходные данные, например коэффициенты уравнений и допустимое значение
погрешности. Необходимо также задать начальные приближения значений
неизвестных. Они либо вводятся в ЭВМ, либо вычисляются каким-либо
способом (в частности, путём решения системы уравнений с помощью прямого
метода). Затем организуется циклический вычислительный процесс, условием
окончания которого является
max ξ i ≤ ε ,
(4.31)
j ≤ i ≤1
Введём обозначения:
ξ i( 0 ) - погрешности значений неизвестных,
ri( 0) - невязки, т.е.
ξ i( 0 ) = x i − x i( 0 ) ,
ξi
где ε- заданная малая положительная величина (точность).
Для предотвращения непроизводительных затрат машинного времени в
алгоритм вводят счётчик числа итераций и при достижении им некоторого
заданного значения счет прекращают.
Блок-схема метода итераций приведена на рис. 41, фигурирующая в
алгоритме переменная р имеет следующий смысл, если р=0, то производится
приведение матрицы А к треугольному виду и определяются коэффициенты
преобразования правых частей системы уравнений.
78
Начало
1
Получение нулей ниже
главной диагонали
ввод n,
матрицы А
u := 0
2
2
k := i+1, n
Приведение матрицы
к треугольному виду
C := -A[k, i]/A[i, i]
3
i: = 1, n
j := i, n
j: = i; k:=0
A[k, j] := A[k, j]+C*A[i, j]
нет
j≤ n
3
да
A[i , j ] > 1E − 6
да
нет
det := 1
j := j+1
k: = j
j: = n+1
i := 1, n
k=0
да
det := det*A[i, i]
нет
нет
k< >i
Перестановка
строк i и k
да
Печать ‘в
столбце’,i,’все
элементы=0’
u mod 2<>0
Exit
det := -det
j: = i, n
C := A[k, j]
A[i, j] := A[k, j]
A[k, j] := C
нет
да
Печать
‘определитель=’, det
конец
u: = u+1
Рис. 41. Блок-схема вычисления определителя матрицы по методу
Гаусса с выбором главного элемента в столбце.
1
79
80
Литература
4.9. Контрольные вопросы
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
Чем отличаются прямые методы от итерационных?
К какому виду приводится матрица коэффициента в прямом ходе метода
Гаусса?
В каком случае нельзя применить метод Гаусса?
В каком порядке определяются неизвестные в обратном ходе метода Гаусса?
Какой элемент является главным в столбце матрицы?
В чём состоит преимущество метода Гаусса с выбором главного элемента в
столбце?
К какому виду приводится матрица в методе Гаусса-Жордана?
Нужен ли обратный ход в методе Гаусса-Жордана?
Для каких систем применителен метод прогонки?
С каким методом схож метод прогонки?
Что является определителем треугольной матрицы?
Что нужно предусмотреть при использовании метода Гаусса?
Каково условие прекращения итерации в итерационных методах?
Основные достоинства метода Гаусса-Зейделя перед методом простых
итераций?
Как проверить являются ли полученные корни истинными или ложными?
81
1.
Бахвалов Н.С., Жидков Н.П., Кобельков Г.Н. Численные методы –
М.:Лаборатория базовых знаний, 2001.
2.
Бахвалов Н.С. Численные методы (часть 1). – М.: Наука, 1983.
3.
Боглаев Ю.П. Вычислительная математика и программирование. - М.:
Высшая школа, 1990.
4.
Вержбицкий В.М. Численные методы. - М.: Высшая школа, 2000.
5.
Волков Е.А. – Численные методы. - М.: Наука, 1982.
6.
Воробьёва Г.Н., Данилова Д. Н. Практикум по вычислительной
математике. - М.: Высшая школа, 1990 .
7.
Демидович Б.П., Марон И.А. и др. Численные методы анализа. - М.:
Наука, 1985.
8.
Самарский А.А. Введение в численные методы. - М.: Наука, 1989.
9.
Самарский А.А., Гулин А.В. Численные методы. - М.: Наука, 1989.
10.
Турчак Л.И. Основы численных методов. - М.: Наука, 1987.
82
Документ
Категория
Информатика
Просмотров
115
Размер файла
673 Кб
Теги
1/--страниц
Пожаловаться на содержимое документа