Учебное пособие: Вычислительная математика

Учебное пособие: Вычислительная математика

Содержание

Введение

Тема 1. Решение задач вычислительными методами. Основные понятия

1.1 Погрешность

1.2 Корректность

1.3 Вычислительные методы

Тема 2. Решение нелинейных уравнений

2.1 Постановка задачи

2.2 Основные этапы отыскания решения

2.3 Метод деления отрезка пополам (метод дихотомии, метод бисекции)

2.4 Метод простых итераций

2.5 Метод Ньютона (метод касательных)

2.6 Метод секущих (метод хорд)

2.7 Метод ложного положения

Тема 3. Решение систем линейных алгебраических уравнений

3.1 Постановка задачи

3.2 Метод исключения Гаусса. Схема единственного деления

3.3 Метод исключения Гаусса с выбором главного элемента по столбцу

3.4 Вычисление определителя методом исключения Гаусса

3.5 Вычисление обратной матрицы методом исключения Гаусса

3.6 Метод простой итерации Якоби

3.7 Метод Зейделя

Тема 4. Приближение функций

4.1 Постановка задачи

4.2 Приближение функции многочленами Тейлора

4.3 Интерполяция функции многочленами Лагранжа

4.4 Аппроксимация функций. Метод наименьших квадратов

Тема 5. Численное интегрирование функций одной переменной

5.1 Постановка задачи численного интегрирования

5.2 Метод средних прямоугольников

5.3 Метод трапеций

5.4 Метод Симпсона (метод парабол)

5.5 Правило Рунге практической оценки погрешности

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

6.1 Постановка задачи Коши

6.2 Метод Эйлера

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

6.4 Метод Рунге – Кутты

Контрольные задания по курсу “Вычислительные методы”

Указания к выполнению лабораторных работ

Указания к выполнению курсовых работ

Краткие сведения о математиках

Список литературы

Введение

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

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

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

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

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

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

Тема 1. Решение задач вычислительными методами.

Основные понятия

1.1 Погрешность

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

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

2. Исходные данные. Исходные данные, как правило, содержат погрешности, так как они либо неточно измерены, либо являются результатом решения некоторых вспомогательных задач. Например, масса снаряда, производительность оборудования, предполагаемая цена товара и др. Во многих физических и технических задачах погрешность измерений составляет 1 – 10%. Погрешность исходных данных так же, как и погрешность математической модели, считается неустранимой и в дальнейшем учитываться не будет.

3. Метод вычислений. Применяемые для решения задачи методы как правило являются приближенными. Например, заменяют интеграл суммой, функцию – многочленом, производную – разностью и т. д. Погрешность метода необходимо определять для конкретного метода. Обычно ее можно оценить и проконтролировать. Следует выбирать погрешность метода так, чтобы она была не более, чем на порядок меньше неустранимой погрешности. Большая погрешность снижает точность решения, а меньшая требует значительного увеличения объема вычислений.

4. Округление в вычислениях. Погрешность округления возникает из-за того, что вычисления производятся с конечным числом значащих цифр (для ЭВМ это 10 – 12 знаков). Округление производят по следующему правилу: если в старшем из отбрасываемых разрядов стоит цифра меньше пяти, то содержимое сохраняемых разрядов не изменяется; в противном случае в младший сохраняемый разряд добавляется единица с тем же знаком, что и у самого числа. При решении больших задач производятся миллиарды вычислений, но так как погрешности имеют разные знаки, то они частично взаимокомпенсируются.

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

D(а*) = | а  –  а*|

называют абсолютной погрешностью числа а*, а величину

d(а*) =

– его относительной погрешностью.

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

1.2 Корректность

Определим вначале понятие устойчивости решения.

Решение задачи y* называется устойчивым по исходным данным x*, если оно зависит от исходных данных непрерывным образом. Это означает, что малому изменению исходных данных соответствует малое изменение решения. Строго говоря, для любого e > 0 существует d = d(e) > 0 такое, что всякому исходному данному x*, удовлетворяющему условию |x  - x*| < d, соответствует приближенное решение y*, для которого  |y – y*| < e.

Говорят, что задача поставлена корректно, если выполнены следующие три условия:

1. Решение существует при любых допустимых исходных данных.

2. Это решение единственно.

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

Если хотя бы одно из этих условий не выполнено, задача называется некорректной.

Пример 1.1.

Покажем, что задача вычисления определенного интеграла  I = корректна. Пусть f*(x) – приближенно заданная функция и I* = . Очевидно, приближенное решение I* существует и единственно. Определим абсолютную погрешность f* с помощью равенства  D(f*) = |f(x) – f*(x)|. Так как

D(I) = |I – I*| = || £  (b – a)D(f*),

то для любого e > 0 неравенство D(I) < e  будет выполнено, если будет выполнено условие D(f*) < d, где d  = e /(b – a).

Таким образом, решение I* устойчиво. Все три условия корректности задачи выполнены.

Пример 1.2.

Покажем, что задача вычисления производной u(x) = f '(x) приближенно заданной функции некорректна.

Пусть f*(x) – приближенно заданная на отрезке [a, b] непрерывно дифференцируемая  функция и u*(x) = (f*(x))'. Определим абсолютные погрешности следующим образом: D(f*) = |f(x) – f*(x)|, D(u*) = |u(x) – u*(x)|.

Возьмем, например, f*(x) = f(x) + a sin(x/a2), где 0 < a < 1. Тогда, u*(x) = u(x) + a-1cos(x/a2), D(u*) = a-1, т. е. погрешность задания функции равна a, а погрешность производной равна a-1. Таким образом, сколь угодно малой погрешности задания функции f может отвечать сколь угодно большая погрешность производной f '.

1.3 Вычислительные методы

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

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

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

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

Оценки погрешности приближения, полученные до вычислений, называют априорными оценками (от лат. a'priori – "до опыта"), а  соответствующие оценки, полученные в ходе вычислений называют апостериорными оценками (от лат. a'posteriori – "после опыта").

Важной характеристикой итерационных методов является скорость сходимости метода. Говорят, что метод имеет p-ый порядок  сходимости если

|xn+1 - x*| = C|xn - x*|p,

где xn и xn+1 – последовательные приближения, полученные в ходе итерационного процесса вычислений, x* – точное решение, C – константа, не зависящая от n. Говорят, что метод сходится со скоростью геометрической прогрессии со знаменателем q < 1, если для всех n справедлива оценка:

|xn - x*| £  Cqn.

Итерационный процесс называется одношаговым, если для вычисления очередного приближения xn+1 используется только одно предыдущее приближение xn и k –шаговым, если для вычисления xn+1 используются k предыдущих приближений xn-k+1, xn-k+2, …, xn.

Тема 2. Решение нелинейных уравнений

2.1 Постановка задачи

Пусть дана некоторая функция f(x) и требуется найти все или некоторые значения x, для которых

f(x) = 0.                                                                                 (2.1)

Значение x*, при котором f(x*) = 0, называется корнем (или решением) уравнения (2.1).

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

Корень x* уравнения (2.1) называется простым, если первая производная функции f(x) в точке x* не равна нулю, т. е. f '(x*)  0. Если же f '(x*) = 0, то корень x* называется кратным корнем.

Геометрически корень уравнения (2.1) есть точка пересечения графика функции y = f(x) с осью абсцисс. На рис. 2.1 изображен график функции y = f(x), имеющей четыре корня: два простых (xи x) и два кратных (xи x).

Рис. 2.1.

Большинство методов решения уравнения (2.1) ориентировано на отыскание простых корней уравнения (2.1).

2.2 Основные этапы отыскания решения

В процессе приближенного отыскания корней уравнения (2.1) обычно выделяют два этапа: локализация (или отделение) корня и уточнение корня.

Локализация корня заключается в определении отрезка [a, b], содержащего один и только один корень. Не существует универсального алгоритма локализации корня. В некоторых случаях отрезок локализации может быть найден из физических соображений. Иногда удобно бывает локализовать корень с помощью построения графика или таблицы значений функции y = f(x). На наличие корня на отрезке [a, b] указывает различие знаков функции на концах отрезка. Основанием для этого служит следующая теорема математического анализа.

Теорема 2.1. Если функция f непрерывна на отрезке [a, b] и принимает на его концах значения разных знаков, так, что f(a)f(b) < 0, то отрезок [a, b] содержит по крайней мере один корень уравнения f(x) = 0.

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

На этапе уточнения корня вычисляют приближенное значение корня с заданной точностью e > 0. Приближенное значение корня уточняют с помощью различных итерационных методов. Суть этих методов состоит в последовательном  вычислении  значений x0, x1, …, xn, …, которые являются приближениями к корню x*.

2.3 Метод деления отрезка пополам (метод дихотомии, метод бисекции)

Метод деления отрезка пополам является самым простым и надежным способом решения нелинейного уравнения.

Пусть из предварительного анализа известно, что корень уравнения (2.1) находится на отрезке [a0, b0], т. е. x*[a0, b0], так, что f(x*) = 0.

Пусть функция f(x) непрерывна на отрезке [a0, b0] и принимает на концах отрезка значения разных знаков, т.е.

f(a0)f(b0) < 0.                                                                                    (2.2)

Разделим отрезок [a0, b0] пополам. Получим точку x0 = . Вычислим значение функции в этой точке: f(x0). Если  f(x0) = 0, то x0 – искомый корень, и задача решена. Если f(x0)0, то f(x0) – число определенного знака: f(x0) > 0, либо f(x0) < 0. Тогда либо на концах отрезка [a0, x0], либо на концах отрезка [x0, b0] значения функции f(x) имеют разные знаки. Обозначим такой отрезок  [a1, b1]. Очевидно, что x*[a1, b1], и длина отрезка [a1, b1] в два раза меньше, чем длина отрезка [a0, b0]. Поступим аналогично с отрезком [a1, b1]. В результате получим либо корень x*, либо новый отрезок [a2, b2], и т.д. (рис. 2.2).

Рис. 2.2

Середина n-го отрезка xn = . Очевидно, что длина отрезка [an, bn] будет равна , а т. к. x*[an, bn], то

| xn – x*| £ £ .                                                              (2.3)

Погрешность метода. Оценка (2.3) характеризует погрешность метода деления отрезка пополам и указывает на скорость сходимости: метод  сходится со скоростью геометрической прогрессии, знаменатель которой q = 1/2. Заметим, что оценка (2.3) является априорной.

Критерий окончания. Из соотношения (2.3) следует, что при заданной точности приближения e  вычисления заканчиваются, когда будет выполнено неравенство  bn –  an < 2e или неравенство  n > log2((b0  –  a0)/e) – 1. Таким образом, количество итераций можно определить заранее. За приближенное значение корня берется величина xn.

Пример 2.1.

Найдем приближенно x =  с точностью = 0.01. Эта задача эквивалентна решению уравнения x5 – 2 = 0, или нахождению нуля функции f(x) = x5 – 2. В качестве начального отрезка [a0, b0] возьмем отрезок [1, 2]. На концах этого отрезка функция принимает значения с разными знаками: f(1) < 0, f(2) > 0.

Найдем число n делений отрезка  [1, 2], необходимых для достижения требуемой точности. Имеем:

| xn – x*| £    =  £  10-2,

n6.

Следовательно, не позднее 6-го деления найдем  с требуемой точностью,  » 1.1484. Результаты вычислений представлены в таблице 2.1.

Таблица 2.1

n	0             1             2             3             4             5            6
an	1.0000   1.0000   1.0000    1.1250   1.1250    1.1406 1.1406
bn	2.0000  1.5000    1.2500    1.2500   1.1875   1.1875  1.1562
xn	1.5000  1.2500   1.1250    1.1875    1.1406    1.1562  1.1484
Зн f(an)	-              -             -              -              -             -             -
Зн f(bn)	+             +            +             +             +             +            +
f(xn)	5.5938  0.7585  -0.2959   0.1812   -0.0691    0.0532  -0.0078
bn – an	1.0000  0.5000   0. 2500   0.1250    0.0625    0.0312   0.0156

2.4 Метод простых итераций

Пусть уравнение (2.1) можно заменить эквивалентным ему уравнением

x = j(x).                                                                                           (2.4)

Например, уравнение   –  0.5 = 0 можно заменить эквивалентным ему уравнением x = 0.5sinx.

Выберем каким-либо образом начальное приближение x0. Вычислим значение функции j(x) при x = x0 и найдем уточненное значение x1 = j(x0). Подставим теперь x1 в уравнение (2.4) и получим новое приближение x2 = j(x1) и т. д. Продолжая этот процесс неограниченно, получим последовательность приближений к корню:

xn+1 = j(xn).                                                                                        (2.5)

Формула (2.5) является расчетной формулой метода простых итераций.

Если последовательность {xn} сходится при n®, т. е. существует

x* =  xn ,                                                                                        (2.6)

и функция j(x) непрерывна, то, переходя к пределу в (2.5) и учитывая (2.6), получим:

x* = xn  = j(x n -1) = j(xn -1) = j(x*).

Таким образом, x* = j(x*), следовательно, x* –  корень    уравнения (2.4).

Сходимость метода. Сходимость метода простых итераций устанавливает следующая теорема.

 

Теорема 2.2. Если в интервале, содержащем корень x*  уравнения (2.4), а также его последовательные приближения x0, x1, …, xn, …, вычисляемые по формуле (2.5), выполнено условие:

|j'(x)| £ q < 1,                                                                              (2.7)

то x* = xn.

т. е. итерационный процесс сходится и справедлива следующая оценка погрешности:

|xn – x*| £  qn|x0 – x*|                                                                          (2.8)

Оценка (2.8) является априорной. Она показывает, что метод простой итерации сходится со скоростью геометрической прогрессии с знаменателем q. Чем меньше q, тем выше скорость сходимости.

Как следует из теоремы 2.2, условие (2.7) является достаточным для сходимости метода простых итераций. Его выполнение гарантирует сходимость процесса (2.5), но невыполнение условия (2.7), вообще говоря, не означает, что итерационный процесс будет расходиться.

На рис. 2.3 – 2.6 показаны четыре случая взаимного расположения линий y = x и y = j(x) и соответствующие итерационные процессы.

Рис. 2.3 и 2.4 соответствуют случаю |j'(x)| < 1, и итерационный процесс сходится. При этом, если j'(x) > 0 (рис. 2.3), сходимость носит односторонний характер, а если j'(x) < 0 (рис. 2.4), сходимость носит двусторонний, колебательный характер. Рис. 2.5 и 2.6 соответствуют случаю |j'(x)|  > 1 – итерационный процесс расходится. При этом может быть односторонняя (рис. 2.5) и двусторонняя (рис 2.6) расходимость.

 

Рис. 2.3                          Рис. 2.4                          Рис. 2.5

Рис. 2.6

 

Погрешность метода. Если известна величина q в условии (2.7), то применима следующая апостериорная оценка погрешности:

|xn – x*| £  |xn – xn – 1|, n > 1.                                                      (2.9)

Критерий окончания. Из оценки (2.9) вытекает следующий критерий окончания итерационного процесса. Вычисления следует продолжать до выполнения неравенства

|xn – xn – 1| < e.

Если это условие выполнено, то можно считать, что xn является приближением к  x* с точностью e.

Если q £ 0.5, то можно пользоваться более простым критерием окончания:

|xn – xn – 1| < e.                                                                                  (2.10)

Пример 2.2.

Используем метод простой итерации для решения уравнения f(x) = sin x – x2 = 0 с точностью e = 0.001.

Преобразуем уравнение к виду (2.4):

x = , т. е. j(x)= .

Нетрудно убедиться, что корень уравнения находится на отрезке [p/6, p/3]. Например, вычислив значения f(x) на концах отрезка, получим: f(p/6)> 0, а f(p/3)< 0, т. е. функция на концах отрезка имеет разные знаки, что в соответствии с теоремой 2.1 указывает на то, что внутри отрезка есть корень. Расположение корня наглядно иллюстрирует рис.2.7.

Рис. 2.7

Подсчитаем, первую и вторую производные функции j(x):

j '(x) = , j"(x) = .

Так как j"(x) > 0 на отрезке [p/6, p/3], то производная j '(x) монотонно возрастает на этом отрезке и принимает максимальное значение на правом конце отрезка, т. е. в точке p/3. Поэтому, справедлива оценка:

|j '(x)| £ |j '(p/3)| » 0.312.

Таким образом, условие (2.7) выполнено, q < 0.5, и можно воспользоваться критерием окончания вычислений в виде (2.10). В табл. 2.2 приведены приближения, полученные по расчетной формуле (2.5). В качестве начального приближения выбрано значение x0 = 1.

Таблица 2.2

n	xn
0 1 2 3 4 5	1 0.8415 0.8861 0.8742 0.8774 0.8765

Критерий окончания выполняется при n = 5, |x5 – x4| < 0.001. Сходимость  двусторонняя, качественный характер такой сходимости представлен на рис. 2.4. Приближенное значение корня  с требуемой точностью x* 0.8765.

2.5 Метод Ньютона (метод касательных)

Метод Ньютона является наиболее эффективным методом решения нелинейных уравнений.

Пусть корень x* Î [a, b], так, что f(a)f(b) < 0. Предполагаем, что функция f(x) непрерывна на отрезке [a, b] и дважды непрерывно дифференцируема на интервале (a, b).  Положим x0 = b. Проведем касательную к графику функции y = f(x) в точке B0 = (x0, f(x0)) (рис. 2.8).

Рис. 2.8

Уравнение касательной будет иметь вид:

y – f(x0) = f '(x0)(x – x0).                                                                    (2.11)

Первое пересечение получим, взяв абсциссу точки пересечения этой касательной с осью OX, т. е. положив в (2.11) y = 0, x = x1:

x1 = x0  – .                                                                               (2.12)

Аналогично поступим с точкой B1(x1, f(x1)), затем с точкой B2(x2, f(x2)), и т. д. в результате получим последовательность приближений x1, x2, …, xn , …,причем

xn +1 = xn  – .                                                                             (2.13)

Формула (2.13) является расчетной формулой метода Ньютона.

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

j(x) = x - .                                                                                (2.14)

 

Сходимость метода. Сходимость метода Ньютона устанавливает следующая теорема.

Теорема 2.3.  Пусть x*  – простой корень уравнения  f(x) = 0, и в некоторой окрестности этого корня функция f дважды непрерывно дифференцируема. Тогда найдется такая малая s-окрестность корня x*, что при произвольном выборе начального приближения x0 из этой окрестности итерационная последовательность, определенная по формуле (2.13) не выходит за пределы этой окрестности и справедлива оценка:

|xn + 1 – x*| £ C |xn – x*|2, n 0,                                                           (2.15)

где С = s -1. Оценка (2.15) означает, что метод сходится с квадратичной скоростью.

Сходимость метода Ньютона зависит от того, насколько близко к корню выбрано начальное приближение. Неудачный выбор начального приближения может дать расходящуюся последовательность. Полезно иметь в виду следующее достаточное условие сходимости метода. Пусть [a, b] – отрезок, содержащий корень. Если в качестве начального приближения x0 выбрать тот из концов отрезка, для которого

f(x)f"(x) ³  0,                                                                                        (2.16)

то итерации (2.13) сходятся, причем монотонно. Рис. 2.8 соответствует случаю, когда в качестве начального приближения был выбран правый конец отрезка: x0 = b.

Погрешность метода. Оценка (2.15) является априорной и неудобна для практического использования. На практике удобно пользоваться следующей апостериорной оценкой погрешности:

 

|xn – x*| £  |xn – xn – 1|.                                                                           (2.17)

 

Критерий окончания. Оценка (2.17) позволяет сформулировать следующий критерий окончания итераций метода Ньютона. При заданной точности e > 0 вычисления нужно вести до тех пор, пока не будет выполнено неравенство

|xn – xn – 1| < e.                                                                                     (2.18)

Пример 2.3.

Применим метод Ньютона для вычисления . где a > 0, p – натуральное число. Вычисление  эквивалентно решению уравнения xp = a. Таким образом, нужно найти корень уравнения f(x) = 0, f(x) = xp – a, f '(x) = pxp – 1. Итерационная формула метода (2.13) примет вид:

xn +1 = xn  –  =  xn  + .                                      (2.19)

Используя формулу (2.19), найдем  с точностью e = 10-3.

xn +1 =  xn  + .

Простой корень уравнения x3 – 7 = 0 расположен на отрезке [1, 2]. Действительно, на концах отрезка [1, 2] функция f(x) = x3 – 7 принимает разные знаки, f (1) < 0,  f (2) > 0. Кроме того, при x = 2 выполнено достаточное условие сходимости (2.16): f (2)f" (2) ³ 0.

Поэтому в качестве начального приближения можно взять x0 = 2.

Результаты приведены в табл. 2.3.

Таблица 2.3

n	xn
0 1 2 3 4 5	2 0.8415 0.8861 0.8742 0.8774 0.8765

2.6 Метод секущих (метод хорд)

В этом и следующем разделе рассмотрим модификации метода Ньютона.

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

f '(xn) » ,

то вместо формулы (2.13) получим

xn +1 = xn –. .     (2.20)

Это означает, что касательные заменены секущими. Метод секущих является двухшаговым методом, для вычисления приближения xn +1  необходимо вычислить два предыдущих приближения  xn и xn – 1 , и, в частности, на первой итерации надо знать два начальных значения x0 и x1.

Формула (2.20) является расчетной формулой метода секущих. На рис. 2.9 приведена геометрическая иллюстрация метода секущих.

Рис. 2.9

Очередное приближение xn +1 получается как точка пересечения с осью OX секущей, соединяющей точки графика функции f(x)  с координатами (xn -1, f(xn - 1)) и (xn , f(xn)).

Сходимость метода. Сходимость метода секущих устанавливает следующая теорема.

 

Теорема 2.4 Пусть x*  – простой корень уравнения  f(x) = 0, и в некоторой окрестности этого корня функция f дважды непрерывно дифференцируема, причем f"(x) ¹ 0. Тогда найдется такая малая s-окрестность корня x*, что при произвольном выборе начальных приближений x0 и x1 из этой окрестности итерационная последовательность, определенная по формуле (2.20) сходится и справедлива оценка:

|xn + 1 – x*| £ C |xn – x*| p,    n ³  0,  p =  » 1.618.                       (2.21)

Сравнение оценок (2.15) и (2.21) показывает, что p < 2, и метод секущих сходится медленнее, чем метод Ньютона. Но в методе Ньютона на каждой итерации надо вычислять и функцию, и производную, а в методе секущих – только функцию. Поэтому при одинаковом объеме вычислений в методе секущих можно сделать примерно вдвое больше итераций и получить более высокую точность.

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

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

|xn – xn – 1| < e.                                                                                           (2.22)

Пример 2.4.

Применим метод секущих для вычисления положительного корня уравнения 4(1 – x2) – ex = 0 с точностью e = 10-3.

Корень этого уравнения находится на отрезке [0, 1], так как f (0) = 3 > 0, а f (1) =  –e  < 0. Подсчитаем вторую производную функции: f "(x) = –8 – ex. Условие f(x)f " (x) ³ 0 выполняется для точки b = 1. В качестве начального приближения возьмем x0 = b = 1. В качестве второго начального значения возьмем x1 = 0.5. Проведем вычисления по расчетной формуле (2.20). Результаты приведены в табл. 2.4.

Таблица 2.4

n	xn
0 1 2 3 4 5	1.0000 0.5000 0.6660 0.7093 0.7033 0.7034

2.7 Метод ложного положения

Рассмотрим еще одну модификацию метода Ньютона.

Пусть известно, что простой корень  x*  уравнения  f(x) = 0 находится на отрезке [a, b] и на одном из концов отрезка выполняется условие f(x)f"(x) ³ 0. Возьмем  эту точку в качестве начального приближения. Пусть для определенности это будет b. Положим x0 = a. Будем проводить из точки B = (b, f(b)) прямые через расположенные на графике функции точки Bn с координатами (xn, f(xn), n = 0, 1, … . Абсцисса точки пересечения такой прямой с осью OX есть очередное приближение xn+1.

Геометрическая иллюстрация метода приведена на рис. 2.10.

Рис. 2.10

Прямые на этом рисунке заменяют касательные в методе Ньютона (рис. 2.8). Эта замена основана на приближенном равенстве

f '(xn) » .                                                                        (2.23)

Заменим в расчетной формуле Ньютона (2.13) производную f '(xn) правой частью приближенного равенства (2.23). В результате получим расчетную формулу метода ложного положения:

xn +1 = xn –..                                                                   (2.24)

Метод ложного положения обладает только линейной сходимостью. Сходимость тем выше, чем меньше отрезок [a, b].

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

|xn – xn – 1| < e.                                                                                     (2.25)

Пример 2.5.

Применим метод ложного положения  для вычисления корня уравнения x3 + 2x – 11 = 0 с точностью e = 10-3.

Корень этого уравнения находится на отрезке [1, 2], так как f (1) =  –8 < 0, а f (2) = 1 > 0. Для ускорения сходимости возьмем более узкий отрезок [1.9, 2], поскольку f (1.9)  < 0, а f (2) > 0. Вторая производная функции f (x) = x3 + 2x – 11 равна 6x. Условие f(x)f"(x) ³ 0 выполняется для точки b = 2. В качестве начального приближения возьмем x0 = a = 1.9. По формуле (2.24) имеем

x1 = x0 –. = 1.9 +  » 1.9254.

Продолжая итерационный процесс, получим результаты, приведенные в табл. 2.5.

Таблица 2.5

n	xn
0 1 2 3	1.9 1.9254 1.9263 1.9263

Тема 3. Решение систем линейных алгебраических уравнений

3.1 Постановка задачи

Требуется найти решение системы линейных уравнений:

a11x1 + a12 x2 + a13x3 + … + a1nxn = b1

a21x1 + a22 x2 + a23x3 + … + a2nxn = b2

a31x1 + a32 x2 + a33x3 + … + a3nxn = b3                                              (3.1)

.

an1x1 + an2 x2 + an3x3 + … + annxn = bn

или в матричной форме:

Ax = b,                                                                                               (3.2)

где

a11 a12 a13 …                    a1n                         x1                b1

a21 a22 a23 …                    a2n                          x2                          b2

A = a31 a32  a33 …           a3n               x =x3  ,   b =b3

an1         an2         an3                         ann                         xn                           bn

По правилу Крамера система n линейных уравнений имеет единственное решение, если определитель системы отличен от нуля (det A 0) и значение каждого из неизвестных определяется следующим образом:

xj = ,   j = 1, …, n,                                                                           (3.3)

где det Aj – определитель матрицы, получаемой заменой j-го столбца матрицы A столбцом правых частей   b.

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

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

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

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

Среди прямых методов наиболее распространенным является метод исключения Гаусса и его модификации, Наиболее распространенными итерационными методами является метод простых итераций Якоби и метод Зейделя.

Эти методы будут рассмотрены в следующих разделах.

3.2 Метод исключения Гаусса. Схема единственного деления

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

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

Прямой ход состоит из n – 1 шагов. На первом шаге исключается переменная x1 из всех уравнений, кроме первого. Для этого нужно из второго, третьего, …, n-го уравнений вычесть первое, умноженное на величину

m = , i = 2, 3, …, n.                                                                         (3.4)

При этом коэффициенты при x1 обратятся в нуль во всех уравнениях, кроме первого.

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

a = aij – ma1j ,  b= bi – mb1.                                                          (3.5)

Легко убедиться, что для всех уравнений, начиная со второго, a= 0, i = 2, 3, …, n. Преобразованная система запишется в виде:

 a11x1 + a12 x2 + a13x3 + … + a1nxn = b1

ax2 + ax3 + … + axn = b

a x2 + ax3 + … + axn = b                                                      (3.6)

ax2 + ax3 + … + axn = b

Все уравнения (3.6), кроме первого, образуют систему (n – 1)-го порядка. Применяя к ней ту же процедуру, мы можем исключить из третьего, четвертого, …, n-го уравнений переменную x2. Точно так же исключаем переменную x3 из последних n – 3 уравнений.

На некотором k-ом шаге в предположении, что главный элемент k-ого шага a0, переменная xk исключается с помощью формул:

m = ,

a = a – ma ,

b= b – mb, i, j = k + 1, k + 2, …, n.                                      (3.7)

Индекс k принимает значения 1, 2, …, n – 1.

При k = n – 1 получим треугольную систему:

a11x1 + a12 x2 + a13x3 + … + a1nxn = b1

ax2 + ax3 + …+ axn = b

ax3 + …+ axn = b                                                      (3.8)

axn = b

с треугольной матрицей An.

Приведение системы (3.1) к треугольному виду (3.8) составляет прямой ход метода Гаусса.

При использовании метода Гаусса нет необходимости в предварительном обосновании существования и единственности решения (т. е. доказательства, что det A ¹ 0). Если на  k-ом шаге все элементы a (i = k, k + 1, …, n) окажутся равными нулю, то система (3.1) не имеет единственного решения.

Обратный ход состоит в вычислении переменных. Из последнего уравнения (3.8) определяем xn... Подставляя его в предпоследнее уравнение, находим xn-1, и т. д. Общие формулы имеют вид:

xn  = ,

xk  = (b- a xk+1 - a xk+2 - … - a xn),  k = n – 1, n – 2, …, 1 (3.9)

Трудоемкость метода. Для реализации метода исключения Гаусса требуется примерно 2/3n3 операций для прямого хода и n2 операций для обратного хода. Таким образом, общее количество операций составляет примерно 2/3n3 + n2.

Пример 3.1.

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

2.0x1 + 1.0x2  –  0.1x3 +  1.0x4  =  2.7

0.4x1 + 0.5x2 +  4.0x3  –  8.5x4 =  21.9

0.3x1  – 1.0x2 + 1.0x3 +  5.2x4  =  – 3.9                                            (3.10)

1.0x1 + 0.2x2 + 2.5x3  –  1.0x4  =  9.9

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

Прямой ход. 1-ый шаг. Вычислим множители:

m =  =  = 0.2; m =  =  = 0.15;  m =  =  = 0.5.

Вычитая из второго, третьего и четвертого уравнений системы (3.10) первое уравнение, умноженное соответственно на m, m, m, получим новую систему:

2.0x1 + 1.0x2  –  0.1x3 + 1.0x4  =  2.7

0.3x2 +  4.02x3  –  8.70x4  =  21.36

–1.15x2 + 1.015x3 + 5.05x4 = – 4.305                                                    (3. 11)

– 0.30x2 + 2.55x3  – 1.50x4  =  8.55

2-ой шаг. Вычислим множители:

m =  =  =  – 3.83333;     m =  =  = –1.0.

Вычитая из третьего и четвертого уравнений системы (3.11) второе уравнение, умноженное соответственно на m и m, приходим к системе:

2.0x1 + 1.0x2  –  0.1x3 + 1.0x4  =  2.7

0.3x2 +  4.02x3  –  8.70x4  =  21.36

16. 425x3  –  28.300x4  = 77.575                                                       (3.12)

6.570x3  – 10.200x4  =  29.910

3-ий шаг. Вычислим множитель:

m =  =  = 0.4.

Вычитая из четвертого уравнения системы (3.12) третье, умноженное на m, приведем систему к треугольному виду:

2.0x1 + 1.0x2  –     0.1x3 +        1.0x4  =  2.7

0.3x2 +  4.02x3  –      8.70x4  =  21.36

16. 425x3  –  28.300x4  = 77.575                                                        (3.13)

1.12x4  =  –1.12

Обратный ход. Из последнего уравнения системы (3.13) находим x4  =  1.000. Подставляя значение x4   в третье уравнение, получим x3 = 2.000. Подставляя найденные значения x4  и x3 во второе уравнение, найдем x2 = 3.000. Наконец, из первого уравнения, подставив в него найденные значения x4, x3 и x2, вычислим x1 = –1.000.

Итак система (3.10) имеет следующее решение:

x1 = 1.000, x2 = 2.000, x3 = 3.000, x4  =   – 1.000.

3.3 Метод исключения Гаусса с выбором главного элемента по столбцу

Хотя метод Гаусса является точным методом, ошибки округления могут привести к существенным погрешностям результата. Кроме того исключение по формулам (3.7) нельзя проводить, если элемент главной диагонали a равен нулю. Если элемент a мал, то велики ошибки округления при делении на этот элемент. Для уменьшения ошибок округления применяют метод исключения Гаусса с выбором главного элемента по столбцу. Прямой ход так же, как и для схемы единственного деления, состоит из n – 1 шагов. На первом шаге прежде, чем исключать переменную x1, уравнения переставляются так, чтобы в левом верхнем углу был наибольший по модулю коэффициент ai1, i =  1, 2, …, n. В дальнейшем, на k-м шаге, прежде, чем исключать переменную xk, уравнения переставляются так, чтобы в левом верхнем углу был наибольший по модулю коэффициент aik, i = k, k + 1, …, n. После этой перестановки исключение переменной xk производят, как в схеме единственного деления.

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

Пример 3.2.

Применим метод исключения Гаусса с выбором главного элемента по столбцу для решения системы уравнений (3.10) из примера 3.1. Прямой ход. 1-ый шаг. Так как коэффициент a11 = 2.0 наибольший из коэффициентов первого столбца, перестановки строк не требуется и 1-ый  шаг полностью совпадает с 1-ым шагом примера 3.1.  Из второго, третьего и четвертого уравнений исключается переменная x1 и система приводится к виду (3.11).

2-ой шаг. Наибольший по модулю коэффициент при x2 в системе (3.11) a = –1.15. Поэтому переставим уравнения следующим образом:

2.0x1 + 1.0x2  –      0.1x3 +      1.0x4  =  2.7

–1.15x2 + 1.015x3 + 5.05x4 = – 4.305                                             (3.14)

0.3x2 +   4.02x3  –   8.70x4  =  21.36

– 0.30x2 +   2.55x3  –   1.50x4  =  8.55

Вычислим множители:

m =  =  =  –0.26087     m =  =  = 0.26087.

Вычитая из третьего и четвертого уравнений системы (3.14) второе уравнение, умноженное соответственно на m и m, приходим к системе:

2.0x1 + 1.0x2  –  0.1x3 + 1.0x4  =  2.7

–1.15x2 + 1.015x3 + 5.05x4 = – 4.305                                             (3.15)

4.28478x3  –  7.38261x4  =  20.23696

2.28522x3  –  2.81739x4  =  9.67305

3-ий шаг. Вычислим множитель:

m =  =  = 0.53333.

Вычитая из четвертого уравнения системы (3.15) третье, умноженное на m, приведем систему к треугольному виду:

2.0x1 + 1.0x2  – 0.1x3 + 1.0x4  =  2.7

–1.15x2 + 1.015x3 + 5.05x4 = – 4.305                                            (3.16)

4.28478x3  –  7.38261x4  =  20.23696

1.11998x4  =  –1.11998

Обратный ход. Обратный ход полностью совпадает с обратным ходом примера 3.1. Решение системы имеет вид:

x1 = 1.000, x2 = 2.000, x3 = 3.000, x4  =   – 1.000.

3.4 Вычисление определителя методом исключения Гаусса

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

det A = (–1)s det An,

где s – число перестановок строк, (s = 0, если использовался метод Гаусса по схеме единственного деления).Таким образом,

det A = (–1)s a11 aa …a                                                           (3.17)

Итак, для вычисления определителя det A необходимо выполнить процедуру прямого хода в методе Гаусса для системы уравнений Ax = 0, затем найти произведение главных элементов, стоящих на диагонали треугольной матрицы и умножить это произведение на (–1)s, где s – число перестановок строк.

Пример 3.3.

Вычислим определитель det A =

2.0    1.0    0.1    1.0

0.4    0.5    4.0    8.5

0.3    1.0    1.0    5.2

1.0    0.2    2.5    1.0

Данный определитель совпадает с определителем системы, рассмотренной в примере 3.1. Он равен произведению диагональных элементов треугольной матрицы (3.13):

det A = 2.0 × 0.30 × 16.425 × 1.12 = 11.0376.

Если же обратиться к примеру 3.2, то, учитывая, что была одна перестановка строк, т.е. s = 1, получим:

det A = (–1) × 2.0 × (–1.15)  × 4.28478 × 1.11998 = 11.0375.

3.5 Вычисление обратной матрицы методом исключения Гаусса

Обратной матрицей к матрице A называется матрица A-1, для которой выполнено соотношение:

A A-1 = E,                                                                                          (3.18)

где E – единичная матрица:

1   0   0  …  0     

0   1   0  …  0

E =    0   0   1  …  0   .                                                                       (3.19)

0   0   0  …  1  

Квадратная матрица A называется невырожденной, если det A ¹ 0. Всякая невырожденная матрица имеет обратную матрицу.

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

Пусть  A – квадратная невырожденная матрица порядка n:

          a11   a12   a13  …  a1n     

          a21   a22   a23  …  a2n     

A =    a31   a32   a33  …  a3n    

          an1   an2   an3  …  ann    

и A-1 – ее обратная матрица:

           x11   x12   x13  …  x1n     

           x21   x22   x23  …  x2n     

A-1 =   x31   x32   x33  …  x3n    

           xn1   xn2   xn3  …  xnn    

Используя соотношения (3.18), (3. 19) и правило умножения матриц, получим систему из n2 уравнений с n2 переменными xij, i, j = 1, 2, …, n. Чтобы получить первый столбец матрицы E, нужно  почленно умножить каждую строку матрицы A на первый столбец матрицы  A-1 и приравнять полученное произведение соответствующему элементу первого столбца матрицы E. В результате получим систему уравнений:

a11x11 + a12 x21 + a13x31 + … + a1nxn1 = 1

a21x11 + a22 x21 + a23x31 + … + a2nxn1 = 0

a31x11 + a32 x21 + a33x31 + … + a3nxn1 = 0                                       (3.20)

an1x11 + an2 x21 + an3x31 + … + annxn1 = 0

Аналогично, чтобы получить второй столбец матрицы E, нужно  почленно умножить каждую строку матрицы A на второй столбец матрицы  A-1 и приравнять полученное произведение соответствующему элементу второго столбца матрицы E. В результате получим систему уравнений:

a11x12 + a12 x22 + a13x32 + … + a1nxn2 = 0

a21x12 + a22 x22 + a23x32 + … + a2nxn2 = 1

a31x12 + a32 x22 + a33x32 + … + a3nxn2 = 0                                          (3.21)

an1x12 + an2 x22 + an3x32 + … + annxn2 = 0

и т. д.

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

Пример 3.4.

Вычислим обратную матрицу A-1 для матрицы

A =    1.8   –3.8    0.7    –3.7

0.7     2.1  –2.6    –2.8

7.3     8.1    1.7    –4.9

1.9   –4.3  –4.3    –4.7

По формулам (3.7) за три шага прямого хода преобразуем матрицу A в треугольную матрицу

1.8              –3.8               0.7              –3.7

0                   3.57778    –2.87222       –1.36111

0                   0               17.73577       19.04992

0                   0                  0                   5.40155

Далее, применим процедуру обратного хода четыре раза для столбцов свободных членов, преобразованных по формулам (3.7) из столбцов единичной матрицы:

1           0             0           0

0           1             0           0

0  ,        0  ,          1  ,        0

0           0             0           1

Каждый раз будем получать столбцы матрицы A-1. Опустив промежуточные вычисления, приведем окончательный вид матрицы A-1:

–0.21121   –0.46003    0.16248     0.26956

–0.03533     0.16873    0.01573   –0.08920

  0.23030     0.04607  –0.00944   –0.19885 .

–0.29316   –0.38837    0.06128     0.18513

3.6 Метод простой итерации Якоби

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

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

Ax = b                                                                                               (3.22)

с квадратной невырожденной матрицей A привести к виду

x = Bx + c,                                                                                       (3.23)

где B – квадратная невырожденная матрица с элементами bij, i, j = 1, 2, …, n, x – вектор-столбец неизвестных xi, c – вектор-столбец с элементами ci, i = 1, 2, …, n.

Существуют различные способы приведения системы (3.22) к виду (3.23). Рассмотрим самый простой. Представим систему (3.22) в развернутом виде:

a11x1 + a12 x2 + a13x3 + … + a1nxn = b1

a21x1 + a22 x2 + a23x3 + … + a2nxn = b2

a31x1 + a32 x2 + a33x3 + … + a3nxn = b3                                            (3.24)

an1x1 + an2 x2 + an3x3 + … + annxn = bn

Из первого уравнения системы (3.24) выразим неизвестную x1:

x1 = a(b1 – a12x2  – a13x3 – … –  a1nxn),

из второго уравнения – неизвестную x2:

x2 = a(b2 – a21x1  – a23x3 – … –  a2nxn),

и т. д. В результате получим систему:

x1 =                b12 x2 + b13x3 + … + b1,n-1xn-1 + b1nxn + c1

x2 = b21x1                         + b23x3 + … + b2,n-1xn-1 + b2nxn  + c2

x3 = b31x1 +  b32 x2+ …                + b3,n-1xn-1 +  b3nxn  + c3                  (3.25)

.

xn=  bn1x1 +  bn2 x2 + bn3x3               + bn,n-1xn-1                       + cn

Матричная запись системы (3.25) имеет вид (3.23). На главной диагонали матрицы B находятся нулевые элементы, а остальные элементы вычисляются по формулам:

bij = ,  ci = , i, j = 1,2, …n, i j.                                           (3.26)

Очевидно, что  диагональные элементы матрицы A должны быть отличны от нуля.

Выберем произвольно начальное приближение Обычно в качестве первого приближения берут x= ci или x= 0. Подставим начальное приближение в правую часть (3.25). Вычисляя левые части, получим значения x, x, …, x. Продолжая этот процесс дальше, получим последовательность приближений, причем (k + 1)-е приближение строится следующим образом:

x  =                b12 x  + b13 x + … + b1,n-1 x  + b1n x + c1

x = b21 x1                         + b23 x  + … + b2,n-1 x   + b2n x   + c2

x= b31 x +  b32 x  + …           + b3,n-1 x   +  b3n x   + c3  …    (3.27)

x=  bn1x +  bn2 x   + bn3 x              + bn,n-1 x                             + c.n

Система (3.27) представляет собой расчетные формулы метода простой итерации Якоби.

Сходимость метода простой итерации. Известно следующее достаточное условие сходимости метода простой итерации Якоби.

Если элементы матрицы  A удовлетворяют условию:

|aii| > , i = 1, 2, …, n.                                                            (3.28)

то итерационная последовательность xk  сходится к точному решению x*.

Условие (3.28) называют условием преобладания диагональных элементов матрицы  A, так как оно означает, что модуль диагонального элемента i-ой строки больше суммы модулей остальных элементов этой строки, i = 1, 2, …, n.

Необходимо помнить, что условие сходимости (3.28) является лишь достаточным. Его выполнение гарантирует сходимость метода простых итераций, но его невыполнение, вообще говоря, не означает, что метод расходится.

Справедлива следующая апостериорная оценка погрешности:

max|x - x| £ max|x– x|,  i = 1, 2, …, n,                            (3.29)

где b  = max |bij| i, j  = 1, 2, …, n.

Правую часть оценки (3.29) легко вычислить после нахождения очередного приближения.

Критерий окончания. Если требуется найти решение с точностью e, то в силу (3.29) итерационный процесс следует закончить как только на (k+1)-ом шаге выполнится неравенство:

max|x– x| < e,  i = 1, 2, …, n.                                              (3.30)

Поэтому в качестве критерия окончания итерационного процесса можно использовать неравенство

max|x– x| < e1,  i = 1, 2, …, n.                                                       (3.31)

где e1 = e.

Если выполняется условие b £ , то можно пользоваться более простым критерием окончания:

max|x– x| < e,  i = 1, 2, …, n.                                                         (3.32)

В других случаях использование критерия (3.32) неправомерно и может привести  к преждевременному окончанию итерационного процесса.

Пример 3.5.

Применим метод простой итерации Якоби для решения системы уравнений

20.9x1 +   1.2 x2 +   2.1x3 +   0.9x4 = 21.70

1.2x1 + 21.2 x2 +   1.5x3 +   2.5x4 = 27.46

2.1x1 +   1.5 x2 + 19.8x3 +  1.3x4 = 28.76                                         (3.33)

0.9x1 +   2.5 x2 +   1.3x3 + 32.1x4 = 49.72

Заметим, что метод простой итерации сходится, т. к. выполняется условие преобладания диагональных элементов (3.28):

|20.9| > |1.2 + 2.1 + 0.9|,

|21.2| > |1.2| + |1.5| + |2.5|,

|19.8| > |2.1| + |1.5| + |1.3|,

|32.1| > |0.9| + |2.5| + |1.3|.

Пусть требуемая точность e = 10-3. Вычисления будем проводить с четырьмя знаками после десятичной точки.

Приведем систему к виду (3.25):

x1  =                   – 0.0574 x2  – 0.1005x3  –  0.0431x4 + 1.0383

x2  = –0.0566x1                      – 0.0708x3 –  0.1179x4 + 1.2953

x3  = –0.1061x1 – 0.0758 x2                                 –  0.0657x4 + 1.4525        (3.34)

x4  = –0.0280x1 – 0.0779 x2   – 0.0405x3                     + 1.5489

Величина b  = max |bij|, i, j  = 1, 2, 3,4 равна 0.1179, т. е. выполняется условие b £ , и можно пользоваться  критерием окончания итерационного процесса (3.32).

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

x = 1.0383, x = 1.2953, x = 1.4525, x = 1.5489.                      (3.35)

Вычисления будем вести до тех пор, пока все величины |x– x|, i = 1, 2, 3, 4, а следовательно, и max|x– x| не станут меньше e = 10-3.

Последовательно вычисляем:

при k = 1

x   =                    – 0.0574x  – 0.1005x  –  0.0431x + 1.0383 = 0.7512

x  = –0.0566x                       – 0.0708x –  0.1179x + 1.2953 = 0.9511

x  = –0.1061x – 0.0758 x                                  –  0.0657x + 1.4525  = 1.1423

x = –0.0280x – 0.0779x   – 0.0405x                        + 1.5489 = 1.3601

при k = 2

x= 0.8106,    x= 1.0118,     x= 1.2117,   x= 1.4077.

при k = 3

x= 0.7978,    x= 0.9977,     x= 1.1975,   x= 1.3983.

при k = 4

x= 0.8004,    x= 1.0005,     x= 1.2005,    x = 1.4003.

Вычисляем модули разностей значений xпри k = 3 и k = 4:

| x– x| = 0.026,  | x– x| = 0.028,  | x– x| = 0.0030,  | x– x| = 0.0020.

Так как все они больше заданной точности e = 10-3, продолжаем итерации.

При k = 5

x= 0.7999, x= 0.9999, x= 1.1999, x = 1.3999.

Вычисляем модули разностей значений xпри k = 4 и k = 5:

| x– x| = 0.0005, | x– x| = 0.0006, | x – x| = 0.0006, | x–  x| = 0.0004.

Все они меньше заданной точности e = 10-3, поэтому итерации заканчиваем. Приближенным решением  системы являются следующие значения:

x1  0.7999,   x2  0.9999,   x3  1.1999,   x4  1.3999.

Для сравнения приведем точные значения переменных:

x1 = 0.8,   x2 = 1.0,   x3 = 1.2,   x4 = 1.4.

3.7 Метод Зейделя

Модификацией метода простых итераций Якоби можно считать метод  Зейделя.

В методе Якоби на (k+1)-ой итерации значения x, i = 1, 2, …, n. вычисляются подстановкой в правую часть (3.27) вычисленных на предыдущей итерации значений x. В методе Зейделя при вычислении xиспользуются значения x, x, x, уже найденные на (k+1)-ой итерации, а не  x, x, …, x, как в методе Якоби, т.е. (k + 1)-е приближение строится следующим образом:

 x  =                   b12 x         + b13 x    + … + b1,n-1 x  + b1n x  + c1

x = b21 x                                       + b23 x    + … + b2,n-1 x  + b2n x   + c2

x= b31 x +   b32 x                        + … + b3,n-1 x  +  b3n x   + c3 (3.36)

x= bn1 x +  bn2 x x  + bn3 x x+ … + bn,n-1 x                          + c.n

Формулы (3.36) являются расчетными формулами метода Зейделя.

Введем нижнюю и верхнюю треугольные матрицы:

           0      0      0  … 0                                     0      b12     b13  …  b1n

          b21    0      0  … 0                                      0       0      b23  …  b2n

B1 =  b31   b32     0  … 0         и    B2  =   0       0          0  …   b3n     .

          bn1   bn2   bn3  …0                                      0       0       0  …    0

Матричная запись расчетных формул (3.36) имеет вид:

xk+1= B1xk+1+ B2xk+ c.                                                                         (3.37)

Так как B = B1+ B2, точное решение x* исходной системы удовлетворяет равенству:

 

x*= B1x*+ B2x*+ c.                                                                             (3.38)

Сходимость метода Зейделя. Достаточным условием сходимости метода Зейделя является выполнение неравенства:

b = max |bij|,< 1,  i, j  = 1, 2, …, n.                                                 (3.39)

Неравенство (3.39) означает, что для сходимости метода Зейделя достаточно, чтобы максимальный по модулю элемент матрицы B был меньше единицы.

Если выполнено условие (3.39), то справедлива следующая апостериорная оценка погрешности:

max|x - x| £ max|x– x|  i = 1, 2, …, n,                           (3.40)

где b – максимальный элемент матрицы B, b2 – максимальный элемент матрицы B2.

Правую часть оценки (3.40) легко вычислить после нахождения очередного приближения.

Критерий окончания. Если требуется найти решение с точностью e, то в силу (3.37) итерационный процесс следует закончить как только на (k+1)-ом шаге выполнится неравенство:

max|x– x| < e,  i = 1, 2, …, n.                                           (3.41)

Поэтому в качестве критерия окончания итерационного процесса можно использовать неравенство

max|x– x| < e1,  i = 1, 2, …, n.                                                  (3.42)

где e1 = e.

Если выполняется условие b £ , то можно пользоваться более простым критерием окончания:

max|x– x| < e,  i = 1, 2, …, n.                                                      (3.43)

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

Пример 3.6.

Применим метод  Зейделя для решения системы уравнений (3.33) из примера 3.5. Первые шаги полностью совпадают с процедурой решения по методу Якоби, а именно: система приводится к виду (3.34), затем в качестве начального приближения выбираются элементы столбца свободных членов (3.35). Проведем теперь итерации методом Зейделя.

При k = 1

x   = – 0.0574x  – 0.1005x  –  0.0431x + 1.0383 = 0.7512

При вычислении xиспользуем уже полученное значение x:

x  = –0.0566 x   – 0.0708x –  0.1179x + 1.2953 = 0.9674

При вычислении x используем уже полученные значения x и x:

x  = –0.1061 x – 0.0758 x     –  0.0657x + 1.4525  = 1.1977

При вычислении x используем уже полученные значения x, x, x:

x = –0.0280 x – 0.0779 x  – 0.0405x x  + 1.5489 = 1.4037

Аналогичным образом проведем вычисления при k = 2 и  k = 3. Получим:

при k = 2

x= 0.8019,    x= 0.9996,     x= 1.9996,   x= 1.4000.

при k = 3

x= 0.80006,    x= 1.00002,     x= 1.19999,   x= 1.40000.

Известны точные значения переменных:

x1 = 0.8,   x2 = 1.0,   x3 = 1.2,   x4 = 1.4.

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

Тема 4. Приближение функций

4.1 Постановка задачи

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

1. Функция задана таблицей в конечном множестве точек, а вычисления нужно произвести в других точках.

2. Функция задана аналитически, но  ее вычисление по формуле затруднительно.

При решении задачи поиска приближенной функции возникают следующие проблемы.

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

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

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

4.2 Приближение функции многочленами Тейлора

Пусть функция y = f(x) определена в окрестности точки a и имеет в этой окрестности n + 1 производную. Тогда в этой окрестности справедлива формула Тейлора:

f(x) = c0 + c1(x – a) + c2(x – a)2 + … + cn(x – a )n + Rn(x) = Tn(x) + Rn(x),

где

ck =

Tn(x) – многочлен Тейлора:

Tn(x)= c0 + c1(x – a) + c2(x – a)2 + … + cn(x – a )n,                                 (4.1)

Rn(x) – остаточный член формулы Тейлора. Его можно записать различными способами, например, в форме Лагранжа:

Rn(x)= , a £ x £ x.

Многочлен Тейлора (4.1) обладает свойством, что в точке x = a все его производные до порядка n включительно совпадают с соответствующими производными функции f, т. е.

T(a)=  f(k)(a),    k =  0, 1, …,  n.

В этом легко убедиться, дифференцируя Tn(x). Благодаря этому свойству многочлен Тейлора хорошо приближает функцию f в окрестности точки a. Погрешность приближения составляет

|f(x) –  Tn(x)| = |Rn(x)|,

т. е. задавая некоторую точность e > 0, можно определить окрестность точки a и значение n из условия:

|Rn(x)| =  < e.                                                                 (4.2)

Пример 4.1.

Найдем приближение функции y = sinx многочленом Тейлора в окрестности точки a = 0. Воспользуемся известным выражением для k-ой производной функции sinx:

(sinx)(k) = sin  x + k                                                                        (4.3)

Применяя последовательно формулу (4.3), получим:

f(0) = sin0 = 0;

f '(0) = cos(0) = 1;

f"(0) = –sin0 = 0;

f(2k-1)(0) =  sin (2k – 1) = (–1)k – 1 ;

f(2k)(0) = 0;

f(2k+1)(x) = (–1)kcosx.

Следовательно, многочлен Тейлора для функции y = sinx для n = 2k имеет вид:

sinx = x – + … + (–1)k – 1 + R2k(x),

R2k(x) = (–1)k  .

Зададим e = 10 –4 и отрезок [–,]. Определим n =2k из неравенства:

|R2k(x)| =  < < < e = 10-4.

Таким образом, на отрезке  –,  функция y = sinx с точностью до e = 10-4 равна многочлену 5-ой степени:

sinx = x – +  = x – 0.1667x3 + 0.0083x5.

Пример 4.2.

Найдем приближение функции y = ex многочленом Тейлора на отрезке [0, 1] с точностью e = 10 –5.

Выберем a = ½, т. е в середине отрезка. При этом величина погрешности в левой части (4.2) принимает минимальное значение. Из математического анализа известно, что для k-ой производной от ex справедливо равенство:

(ex)(k) = ex.

Поэтому

(ea)(k)  = ea = e1/2,

Следовательно, многочлен Тейлора для функции y = ex имеет вид:

ex  = e1/2 + e1/2(x – ½) +  (x – ½)2 + … +  (x – ½)n+ Rn(x),

При этом, учитывая, что xÎ [0, 1], получим оценку погрешности:

|Rn(x)|  < .                                                                     (4.4)

Составим таблицу погрешностей, вычисленных по формуле (4.4):

n	2	3	4	5	6
Rn	0.057	0.0071	0.00071	0.000059	0.0000043

Таким образом, следует взять n = 6.

4.3 Интерполяция  функции многочленами Лагранжа

Рассмотрим другой подход к приближению функции многочленами. Пусть функция y = f(x) определена на отрезке [a, b] и известны значения этой функции в некоторой системе узлов xi Î [a, b], i = 0, 1, … , n. Например, эти значения получены в  эксперименте при наблюдении некоторой величины в определенных точках или в определенные моменты времени x0, x1, … , xn. Обозначим эти значения следующим образом: yi = f(xi), i = 0, 1, … , n. Требуется найти такой многочлен P(x)  степени m,

P(x) =  a0 + a1x + a2x2 + … + amxm,                                                (4.5)

который бы в узлах xi, i = 0, 1, … , n принимал те же значения, что и исходная функция y = f(x), т. е.

P(xi) = yi,  i = 0, 1, … , n.                                                                   (4.6)

Многочлен (4.5), удовлетворяющий условию (4.6), называется интерполяционным многочленом.

Другими словами, ставится задача построения функции y = P(x), график которой проходит через заданные точки (xi, yi), i = 0, 1, … , n (рис. 4.1).

Рис. 4.1

Объединяя (4.5) и (4.6), получим:

a0 + a1xi + a2x + … + amx = yi,      i = 0, 1, … , n.                         (4.7)

В искомом многочлене  P(x) неизвестными являются m +1 коэффициент a0 , a1, a2, …, am. Поэтому систему (4.7) можно рассматривать как систему из n +1 уравнений с m +1 неизвестными. Известно, что  для существования единственного решения такой системы необходимо , чтобы выполнялось условие: m = n. Таким образом, систему (4.7) можно переписать в развернутом виде:

a0 + a1 x0 + a2x + … + anx = y0

a0 + a1 x1 + a2x + … + anx = y1

a0 + a1 x2 + a2x + … + anx = y2                                                       (4.8)

.

a0 + a1 xn + a2x + … + anx = yn

Вопрос о существовании и единственности интерполяционного многочлена решает следующая теорема:

 

Теорема 4.1. Существует единственный интерполяционный многочлен степени n, удовлетворяющий условиям (4.6).

Имеются различные формы записи интерполяционного многочлена. Широко распространенной формой записи является многочлен Лагранжа

Ln(x) =  = .   (4.9)

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

L1(x) = y0+ y1,

L2(x) = y0+ y1+ y2.

Пример 4.3.

Построим интерполяционный многочлен Лагранжа по следующим данным:

x	0	2	3	5
y	1	3	2	5

Степень многочлена Лагранжа для n +1 узла равна n. Для нашего примера многочлен Лагранжа имеет третью степень. В соответствии с (4.9)

L3(x) = 1+3 + 2 + 5 = 1 + x –  x2 + x3.

Пример 4.4.

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

Для функции y = sinx известны следующие данные.

x	0	p/6	p/3	p/2
y	0	½		1

Вычислим y(0.25).

Найдем многочлен Лагранжа третьей степени:

L3(x) = 0 + +

 + 1.

При x = 0.25 получим y(0.25) = sin 0.25 » 0.249.

Погрешность интерполяции. Пусть интерполяционный многочлен Лагранжа построен для известной функции f(x). Необходимо выяснить, насколько этот многочлен близок к функции в точках отрезка [a, b], отличных от узлов. Погрешность интерполяции равна |f(x) – Pn(x)|. Оценку погрешности можно получить на основании следующей теоремы.

Теорема 4.2. Пусть функция f(x) дифференцируема n +1 раз на отрезке [a, b], содержащем узлы интерполяции xi Î [a, b], i = 0, 1, … , n. Тогда для погрешности интерполяции в точке x Î [a, b] справедлива оценка:

|f(x) – Ln(x)|£ |wn+1(x)|,                                                         (4.10)

где

Mn+1 = |f(n+1)(x)|,

wn+1(x) = (x – x0)(x – x1)…. (x – xn).

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

|f(x) – Ln(x)| £  |wn(x)|                                                 (4.11)

Пример 4.5.

Оценим погрешность приближения функции f(x) =  в точке x = 116  и на всем отрезке [a, b], где a = 100, b = 144, с помощью интерполяционного много члена Лагранжа L2(x) второй степени, построенного с узлами x0 = 100, x2 = 144.

Найдем первую, вторую и третью производные функции f(x):

f '(x)= x – 1/2,   f "(x)= – x –3/2,    f'''(x)= x –5/2.

M3 = | f'''(x)| =  100 –5/2 =  10 –5.

В соответствии с (4.9) получим оценку погрешности в точке x = 116:

| –  L2(116)| £ |(116 – 100)(116 – 121)(116 – 144)| = 10 –5×16×5×28 = 1.4×10 – 3.

Оценим погрешность приближения функции f(x) =  на всем отрезке в соответствии с (4.11):

| – L2(x)| £ |(x – 100)(x – 121)(x –144)| » 2.5×10–3.

4.4 Аппроксимация функций. Метод наименьших квадратов

В инженерной деятельности часто возникает необходимость описать в виде функциональной зависимости связь между величинами, заданными таблично или в виде набора точек с координатами (xi, yi), i = 0, 1, 2,... , n, где n –  общее количество точек. Как правило, эти табличные данные получены экспериментально и имеют погрешности (рис. 2.5)

Рис.4.2

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

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

S =,                                                              (4.12)

обращается в минимум.

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

D = .                                                                              (4.13)

В качестве функциональной зависимости рассмотрим многочлен

Pm(x)=a0 + a1x + a2x2+...+amxm.                                                  (4.14)

Формула (4.12) примет вид

S =

Условия минимума S можно записать, приравнивая нулю частные производные S по  всем переменным a0,  a1,  a2, … , am. Получим систему уравнений

 =  –= 0, или

= 0,    k = 0, 1, … , m.                 (4.15) 

Систему уравнений (4.15) перепишем в следующем виде:

a0+ a1+ … +am= , k = 0, 1, … , m     (4.16)

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

ck = , bk = .

Система (4.16) может быть записана так:

a0ck + a1ck+1 + … + ck+mam = bk, k = 0, 1, … , m.                                   (4.17)

Перепишем систему (4.17) в развернутом виде:

  c0a0 + c1a1 + c2a2… + cmam = b0

  c1a0 + c2a1 + c3a2… + cm+1am  = b1

(4.18)

cma0 + cm+1a1 + cm+2a2… + c2mam  = bm

Матричная запись системы (4.18) имеет следующий вид:

Ca = b.                                                                                          (4.19)

Для определения коэффициентов ak, k = 0, 1, … , m, и, следовательно, искомого многочлена (4.14) необходимо вычислить суммы ck, bk и решить систему уравнений (4.18). Матрица  C системы (4.19) называется матрицей Грама и является симметричной и положительно определенной. Эти полезные свойства используются при решении.

Погрешность приближения в соответствии с формулой (4.13) составит

D = .                                                                (4.20)

Рассмотрим частные случаи m =1 и m = 2.

1. Линейная аппроксимация (m = 1).

P1(x) = a0 + a1x.

c0 = = n + 1; c1 = = ; c2 = ;                         (4.21)

b0 = = ; b1 = = .                                  (4.22)

           c0           c1                         n+1        

C =                             =                                      ,

           c1               c2                                   

b = (b0, b1)T = (,)T.

Решение системы уравнений  Ca = b найдем по правилу Крамера:

a0 = , a1 = ,

где úCú  –  определитель матрицы C, аúCiú – определитель матрицы Ci, полученной из матрицы C заменой i-го столбца столбцом свободных членов b, i = 1, 2.

Таким образом,

a0 = ,  a1 = .                                                    (4.23)

Алгоритм 4.1 (Алгоритм метода наименьших квадратов. Линейная аппроксимация).

Шаг 1. Ввести исходные данные: xi, yi, i=0, 1, 2, ... , n.

Шаг 2. Вычислить коэффициенты c0,  c1, b0, b1 по формулам (4.21), (4.22).

Шаг 3. Вычислить a0, a1 по формулам (4.23).

Шаг 4. Вычислить величину погрешности

D1 = .                                                          (4.24)

Шаг 5. Вывести на экран результаты: аппроксимирующую линейную функцию P1(x) = a0 + a1x и величину погрешности D1.

2. Квадратичная аппроксимация (m = 2).

P2(x) = a0 + a1x + a2x2.

c0 == n+1; c1 ==; c2 =; c3 =; c4 =.  (4.25)

b0 ==; b1 ==; b2 = .                      (4.26)

           c0    c1    c2                 

C =     c1     c2    c3    .

           c2    c3    c4

b = (b0, b1, b2)T .

Решение системы уравнений  Ca = b найдем по правилу Крамера:

ai = , i = 0, 1,

где úCú  –  определитель матрицы C, аúCiú – определитель матрицы Ci, полученной из матрицы C заменой i-го столбца столбцом свободных членов b.

úCú = c0c2c4 + 2c1c2c3  –  c –  сc4  – cc0.                                       (4.27)

               b0    c1  c2 

úC1ú =   b1    c2   c3  = b0c2c4 + b2c1c3 + b1c2c3 – b2c– b1c1c4 – b0c.    (4.28)

              b2    c3   c4

             c0    b0   c2   

úC2ú =  c1    b1    c3 = b1c0c4 + b0c2c3 + b2c1c2 – b1c– b0c1c4 – b2c0c3.   (4.29)

             c2    b2   c4

  c0    c1    b0 

úC3ú =  c1    c2    b1  = b2c0c2 + b1c1c2 + b0c1c3 – b0c– b2c – b1c0c3.    (4.30)

c2    c3    b2

a0 = , a1 = , a2 = .                                                                 (4.31)

 

Алгоритм 4.2 (Алгоритм метода наименьших квадратов. Квадратичная аппроксимация).

Шаг 1. Ввести исходные данные: xi, yi, i=0, 1, 2, ... , n.

Шаг 2. Вычислить коэффициенты c0,  c1, c2, c3, c4, b0, b1, b2, по формулам (4.25), (4.26).

Шаг 3. Вычислить úCú, úC1ú, úC2ú, úC3ú по формулам (4.27) – (4.30).

Шаг 4. Вычислить a0, a1, a2 по формулам (4.31).

Шаг 5. Вычислить величину погрешности

D2 = .                                              (4.32)

Шаг 5. Вывести на экран результаты : аппроксимирующую квадратичную функцию P2(x) = a0 + a1x + a2x2  и величину погрешности D2.

Пример 4.6.

Построим по методу наименьших квадратов многочлены первой и второй степени и оценим степень приближения. Значения yi в точках xi , i =0, 1, 2, 3, 4 приведены в таблице 2.3.

Таблица 4.1

i	0	1	2	3	4
xi	1	2	3	4	5
yi	–1	1	2	4	6

Вычислим коэффициенты c0,  c1, c2, c3, c4, b0, b1, b2, по формулам (4.25), (4.26):

c0 = 5; c1 = 15; c2 = 55; c3 = 225; c4 = 979;

b0 = 12;  b1 = 53;  b2 = 235.

1. Линейная аппроксимация (m =1).

Система уравнений для определения коэффициентов a0 и a1 многочлена первой степени P2(x) = a0 + a1x + a2x2  имеет вид

5a0 + 15a1 = 12

15a0 + 55a1 = 53

По формулам (4.23) найдем  коэффициенты a0 и a1:

a0 =  » –2.7,  a1 = » 1.7.    

P1(x) = a0 + a1x =  –2.7 + 1.7x.

2. Квадратичная аппроксимация (m =2).

Система уравнений для определения коэффициентов a0, a1 и a2 многочлена второй степени P2(x) = a0 + a1x + a2x2  имеет вид

5a0 +   15a1 +   55a2 =   12

15a0 +   55a1 + 225a2 =   53

55a0 + 225a1 + 979a2 = 235

По формулам (4.31) найдем  коэффициенты a0, a1 и a2:

a0 » –2.20, a1 » 1.27, a2 » 0.07.

P2(x) = a0 + a1x + a2x2 = –2.20 + 1.27x + 0.07x2.

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

Таблица 4.2

i	0	1	2	3	4
xi	1	2	3	4	5
yi	–1	1	2	4	6
P1(xi)	–1	0.7	2.4	4.1	5.8
P2(xi)	–1	0.62	2.24	4	6.9

Погрешность приближения в соответствии с формулами (4.24) и (4.32) составит

D1 =  = 0.245.

D2 =  = 0.084.

Тема 5. Численное интегрирование функций одной переменной

5.1 Постановка задачи численного интегрирования

Далеко не все интегралы можно вычислить по известной из математического анализа формуле Ньютона – Лейбница:

I == F(b) – F(a),                                (5.1)

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

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

» ,                                                                                  (5.2)

где xi – некоторые точки на отрезке [a, b],называемые узлами квадратурной формулы, Ai – числовые коэффициенты, называемые весами квадратурной формулы, n ³ 0 – целое число.

5.2 Метод прямоугольников

Формулу  прямоугольников можно получить из геометрической интерпретации интеграла. Будем интерпретировать интеграл  как площадь криволинейной трапеции, ограниченной графиком функции y = f(x), осью абсцисс и прямыми x = a и x = b (рис. 5.1).

Рис. 5.1

Разобьем отрезок [a, b] на n равных частей длиной h, так, что h = . При этом получим точки a = x0 < x1< x2 < … < xn = b и  xi+1 = xi  + h, i = 0, 1, … , n – 1 (рис. 5.2)

Рис. 5.2

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

Рис. 5.3

Эта фигура состоит из n прямоугольников. Основание i-го прямоугольника образует отрезок [xi, xi+1] длины h, а высота основания равна значению функции в середине отрезка [xi, xi+1],  т е. f(рис. 5.4).

Рис. 5.4

Тогда получим квадратурную формулу средних прямоугольников:

I =» Iпр =                                                    (5.3)

Формулу (5.3) называют также формулой средних прямоугольников. Иногда используют формулы

I  » I =  ,                                                                            (5.4)

I  » I =  ,                                                                            (5.5)

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

Геометрические иллюстрации этих формул приведены на рис. 5.5 и 5.6.

Рис. 5.5

Рис. 5. 6

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

Теорема 5.1. Пусть функция f дважды непрерывно дифференцируема на отрезке [a, b]. Тогда для формулы прямоугольников справедлива следующая оценка погрешности:

| I  –  Iпр | £  h2,                                                                          (5.6)

где M2  = |f "(x)|

Пример 5.1.

Вычислим значение интеграла  по формуле средних прямоугольников (5.3) с шагом h = 0.1.

Составим таблицу значений функции e(табл. 5.1):

Таблица 5.1

x	e	x	e
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50	1.0000000 0.9975031 0.9900498 0.9777512 0.9607894 0.9394131 0.9139312 0.8847059 0.8521438 0.8166865 0.7788008	0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00	0.7389685 0.6976763 0.6554063 0.6126264 0.5697828 0.5272924 0.4855369 0.4448581 0.4055545 0.3678794

Производя вычисления по формуле (5.3), получим:

Iпр  = 0.74713088.

Оценим погрешность полученного значения. Имеем:

f "(x) = (e)" = (4x2  –  2) e.

Нетрудно убедиться, что | f "(x)| £ M2 = 2. Поэтому по формуле(5.4)

| I  –  Iпр | £ (0.1)2  » 0.84× 10-3.

5.3 Метод трапеций

Выведем формулу трапеций так же, как и формулу прямоугольников, из геометрических соображений. Заменим график функции y = f(x) (рис.5.1) ломаной линией (рис.5.7), полученной следующим образом. Из точек a = x0, x1, x2,…, xn = b проведем ординаты до пересечения с кривой y = f(x). Концы ординат соединим прямолинейными отрезками.

Рис. 5.7

Тогда площадь криволинейной трапеции приближенно можно считать равной площади фигуры, составленной из трапеций. Так как площадь трапеции, построенной на отрезке [xi, xi+1] длины h =  , равна h , то, пользуясь этой формулой для i = 0, 2, … , n – 1, получим квадратурную формулу трапеций:

I=»Iтр =h= (5.7)

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

Теорема 5.2. Пусть функция f дважды непрерывно дифференцируема на отрезке [a, b]. Тогда для формулы трапеций справедлива следующая оценка погрешности:

| I  –  Iтр | £  h2,                                                                         (5.8)

где M2  = |f "(x)|.

Пример 5.2.

Вычислим значение интеграла по формуле трапеций (5.7) и сравним полученный результат с результатом примера 5.1.

Используя таблицу значений функции eиз примера 5.1 и производя вычисления по формуле трапеций (5.7), получим: Iтр = 0.74621079.

Оценим погрешность полученного значения. В примере (5.1) получили оценку:  | f "(x)| £ M2 = 2. Поэтому по формуле (5.8)

 I  –  Iтр | £ (0.1)2  » 1.7× 10-3.

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

5.4 Метод Симпсона (метод парабол)

Заменим график функции y = f(x) на отрезке  [xi, xi+1], i = 0, 2, … , n – 1, параболой, проведенной через точки (xi, f(xi)), (x,f(x)), (xi+1, f(xi+1)), где x - середина отрезка [xi, xi+1]. Эта парабола есть интерполяционный многочлен второй степени L2(x) с узлами xi, x, xi+1. Нетрудно убедиться, что уравнение этой параболы имеет вид:

y = L2(x) =

f(x) + (x – x) + (x - x)2,    (5.9)

где h = .

Проинтегрировав функцию (5.9) на отрезке  [xi, xi+1], получим

Ii = »  = ( f(xi) + 4f(x) + f(xi+1)).                     (5.10)

Суммируя  выражение (5.10) по i = 0, 1, 2, … , n – 1, получим квадратурную формулу Симпсона (или формулу парабол):

 I =» IС = ( f(x0) + f(xn) + 4 + 2).           (5.11)

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

 

Теорема 5.2. Пусть функция f имеет на отрезке [a, b] непрерывную производную четвертого порядка f (4)(x). Тогда для формулы Симпсона (5.9) справедлива следующая оценка погрешности:

| I  –  IС | £ h4,                                                                      (5.12)

где M4 = | f (4)(x)|.

Замечание. Если число элементарных отрезков, на которые делится отрезок  [a, b], четно , т.е. n = 2m, то параболы можно проводить через узлы с целыми индексами, и вместо элементарного отрезка  [xi, xi+1] длины h рассматривать отрезок [x2i, x2i+2] длины 2h. Тогда формула Симпсона примет вид:

I » (f(x0) + f(x2m) + 4 + 2),                                   (5.13)

а вместо оценки (5.10) будет справедлива следующая оценка погрешности:

| I  –  IС | £  h4,                                                                    (5.14)

Пример 5.3.

Вычислим значение интеграла  по формуле Симпсона (5.11) и сравним полученный результат с результатами примеров 5.1 и 5.2.

Используя таблицу значений функции eиз примера 5.1 и производя вычисления по формуле Симпсона (5.11) , получим:

IС = 0.74682418.

Оценим погрешность полученного значения. Вычислим четвертую производную f (4)(x).

f (4)(x) = (16x4 – 48x2 + 12) e, | f (4)(x)| £ 12.

Поэтому

| I  –  IС | £ (0.1)4 » 0.42 × 10-6.

Сравнивая результаты примеров 5.1, 5.2 и 5.3, видим , что метод Симпсона имеет меньшую погрешность, чем метод средних прямоугольников и метод трапеций.

5.5 Правило Рунге практической оценки погрешности

Оценки погрешности по формулам (5.4), (5.8) и (5.12) являются априорными. Они зависят от длины элементарного отрезка h, и при достаточно малом h справедливо приближенное равенство:

 I  –  Ih  » Chk,                                                                                     (5.15)

где Ih приближенное значение интеграла, вычисленное по одной из формул (5.3), (5.5), (5.9), C ¹ 0 и k > 0 – величины, не зависящие от h.

Если уменьшить шаг h в два раза, то, в соответствии с (5.15) получим:

 I  –  Ih/2 » Chk » ( I  –  Ih).                                                              (5.16)

Непосредственное использование оценок  погрешности (5.4), (5.8) и (5.12) неудобно, так как при этом требуется вычисление производных функции f (x). В вычислительной практике используются другие оценки.

Вычтем из равенства (5.15) равенство (5.16):

Ih/2  –  Ih  » Chk(2k – 1).                                                                      (5.17)

Учитывая приближенное равенство (5.16), получим следующее приближенное равенство:

I  –  Ih/2 »  .                                                                                (5.18)

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

Для формул прямоугольников и трапеций k = 2, а для формулы Симпсона k = 4. Поэтому для этих формул приближенное равенство (5.18) принимает вид:

I  –  Iпр  » ,                                                                           (5.19)

I  –  Iтр  » ,                                                                            (5.20)

I  –  IС  » .                                                                             (5.21)

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

Пример 5.4.

Найдем значение интеграла   с точностью e = 10-4, используя формулу трапеций и применяя вышеизложенную процедуру дробления шага. В примере 5.2 было получено значение I при h1 = 0.1, Ih =0.74621079. Уменьшим шаг вдвое: h2 = 0.05 и вычислим I= 0.74667084, e2  = ( I- I) = (0.74667084 – 0.74621079) » 1.5×10-4. Так как |e2| > e, то снова дробим шаг: h3 = 0.025, вычисляем I= 0.74678581, e2  = ( I- I) = (0.74678581 – 0.74667084) » 4×10-5. Поскольку  |e3| < e,  требуемая точность достигнута и I » 0.7468 ± 0.0001.

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

6.1 Постановка задачи Коши

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

y' (t) = f(t, y(t)).                                                                                      (6.1)

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

Рис. 6.1

Производную y'(t) в каждой точке (t, y) можно геометрически интерпретировать как тангенс угла a наклона касательной к графику решения, проходящего через эту точку, т е.: k = tga = f(t, y).

Уравнение (6.1) определяет целое семейство решений. Чтобы выделить одно решение, задают начальное условие:

y(t0 ) = y0,                                                                                                  (6.2)

где t0 – некоторое заданное значение аргумента t, а y0 – начальное значение функции.

Задача Коши заключается в отыскании функции y = y(t), удовлетворяющей уравнению (6.1) и начальному условию (6.2). Обычно определяют решение задачи Коши на отрезке, расположенном справа от начального значения t0, т. е.  для t Î [t0, T].

Разрешимость задачи Коши определяет следующая теорема.

 

Теорема 6.1.  Пусть функция f(t, y) определена и непрерывна при t0 £  t £  T, -¥ < y < ¥  и удовлетворяет условию Липшица:

| f(t, y1) – f(t, y2)| £ L| y1 – y2|,

где L некоторая постоянная, а  y1 , y2  –  произвольные значения.

Тогда для каждого начального значения y0 существует единственное решение y(t) задачи Коши для t Î [t0, T].

Даже для простых дифференциальных уравнений первого порядка не всегда удается получить аналитическое решение. Поэтому большое значение имеют численные методы решения. Численные методы позволяют определить приближенные значения искомого решения y(t) на некоторой выбранной сетке значений аргумента ti, (i = 0, 1, …). Точки ti называются узлами сетки, а величина hi = ti+1 –  ti –  шагом сетки. Часто рассматривают равномерные сетки, для которых шаг hi постоянен, hi = h = . При этом решение получается в виде таблицы, в которой каждому узлу сетки ti соответствуют приближенные значения функции y(t) в узлах сетки yi » y(ti).

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

Сходимость численных методов решения задачи Коши.  Пусть y(t) – решение задачи Коши. Назовем глобальной погрешностью (или просто погрешностью) численного метода функцию ei = y(ti) –  yi , заданную в узлах сетки ti. В качестве абсолютной погрешности примем величину R = | y(ti) –  yi|

Численный метод решения задачи Коши называется сходящимся, если для него R ® 0 при h ® 0. Говорят, что метод имеет p-ый порядок точности, если для погрешности справедлива оценка R £ Chp, p > 0, C – константа, C ¹ 0.

6.2 Метод Эйлера

Простейшим методом решения задачи Коши является метод Эйлера.

Будем решать задачу Коши

y' (t) = f(t, y(t)).

y(t0 ) = y0,

на отрезке [t0, T]. Выберем шаг h = , и построим  сетку с системой узлов

ti = t0 + ih, i = 0, 1, …, n.

В методе Эйлера вычисляются приближенные значения функции y(t) в узлах сетки :yi » y(ti).

Заменив производную y' (t) конечными разностями на отрезках [ti, ti+1],  i = 0, 1, …, n – 1, получим приближенное равенство:

 = f(ti, yi),          i = 0, 1, …, n – 1,

которое можно переписать  так:

yi+1 = yi + h f(ti, yi),          i = 0, 1, …, n – 1.                                     (6.3)

Формулы (6.3)  и начальное условие (6.2) являются расчетными формулами метода Эйлера.

Геометрическая интерпретация одного шага метода Эйлера заключается в том, что решение на отрезке [ti, ti+1] заменяется касательной y = y' (ti)( t - ti), проведенной в точке (ti, y(ti)) к интегральной кривой, проходящей через эту точку. После выполнения n шагов неизвестная интегральная кривая заменяется ломаной линией (ломаной Эйлера).

Оценка погрешности. Для оценки погрешности метода Эйлера воспользуемся следующей теоремой.

 

Теорема 6.2. Пусть функция f удовлетворяет условиям:

 £ K,            =  £ L.                                                (6.4)

Тогда для метода Эйлера справедлива следующая оценка погрешности:

R = | y(ti) –  yi| £  = ,

где l – длина отрезка  [t0, T]. Мы видим , что метод Эйлера имеет первый порядок точности.

Оценка погрешности метода Эйлера часто бывает затруднительна, так как требует  вычисления производных функции f(t, y(t)). Грубую оценку погрешности дает правило Рунге (правило двойного пересчета), которое используется для различных одношаговых методов, имеющих p-ый порядок точности. Правило Рунге заключается в следующем. Пусть y – приближения, полученные  с шагом , а y – приближения, полученные  с шагом h. Тогда справедливо приближенное равенство:

|y- y(ti)|  » |y- y| .                                                                   (6.5)

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

R  » |y- y|                                                                                    (6.6)

Так как метод Эйлера имеет первый порядок точности, т. е. p = 1, то приближенное равенство (6.6) примет вид

R  » |y- y|                                                                                             (6.7)

Используя правило Рунге, можно построить процедуру приближенного вычисления решения задачи Коши с заданной точностью e. Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y, i = 0, 1, …, n. Вычисления прекращаются тогда, когда  будет выполнено условие:

R  » |y- y| < e.                                                                           (6.8)

Для метода Эйлера условие (6.8) примет вид

R » |y- y| < e                                                                                       (6.9)

Приближенным решением будут значения y, i = 0, 1, …, n.

Пример 6.1.

Найдем решение на отрезке [0, 1] следующей задачи Коши:

y' (t) = y – ,                                                                                          (6.10)

y(0) = 1.

Возьмем шаг h = 0.2. Тогда   n =  = 5.

В соответствии с (6.3) получим расчетную формулу метода Эйлера:

yi+1 = yi + 0.2  ,  y0 = 1,        i = 0, 1, 2, 3, 4, 5.

Решение представим в виде таблицы 6.1:

Таблица 6.1

i	0	1	2	3	4	5
ti	0	0.2	0.4	0.6	0.8	1.0
yi	1.0000	1.2000	1.3733	1.5294	1. 6786	1.8237

Уравнение (6.10) есть уравнение Бернулли. Его решение можно найти в явном виде:

y = .                                                                                        (6.11)

Для сравнения точного и приближенного решений представим точное решение (6.11) в виде таблицы 6.2:

Таблица 6.2

i	0	1	2	3	4	5
ti	0	0.2	0.4	0.6	0.8	1.0
y(ti)	1.0000	1.1832	1.3416	1.4832	1. 6124	1.7320

Из таблицы видно, что погрешность составляет R = | y(ti) –  yi| = 0.0917.

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

Первый модифицированный метод Эйлера. Суть этого метода состоит в следующем. Сначала вычисляются вспомогательные значения искомой функции y  в точках t = ti +  с помощью формулы:

y = yi +  fi  = yi +f(ti, yi).

Затем находится значение правой части уравнения (6.1) в средней точке

f = f(t, y)

и затем полагается

yi+1 = yi + h f,          i = 0, 1, …, n – 1.                                     (6.12)

Формулы (6.12) являются расчетными формулами первого модифицированного метода Эйлера.

Первый модифицированный метод Эйлера является одношаговым методом со вторым порядком точности

Второй модифицированный метод Эйлера – Коши. Суть этого метода состоит в следующем. Сначала вычисляются вспомогательные значения

 = yi + h f(ti, yi).                                                                               (6.13)

Затем приближения искомого решения находятся по формуле:

yi+1 = yi + [f(ti, yi) + f(ti+1, )],   i = 0, 1, …, n – 1.                           (6.14)

Формулы (6.14) являются расчетными формулами второго модифицированного метода Эйлера – Коши.

Второй модифицированный метод Эйлера – Коши, так же, как и первый,  является одношаговым методом со вторым порядком точности.

Оценка погрешности. Приближенная оценка погрешности модифицированных методов Эйлера осуществляется как и для простого метода Эйлера с использованием правила Рунге (см. предыдущий раздел 6.2). Так как оба модифицированных метода Эйлера имеют второй порядок точности, т. е. p = 2, то оценка погрешности (6.6) примет вид

R  » |y- y|.                                                                                       (6.15)

Используя правило Рунге, можно построить процедуру приближенного вычисления решения задачи Коши модифицированными методами Эйлера с заданной точностью e.  Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y, i = 0, 1, …, n. Вычисления прекращаются тогда, когда  будет выполнено условие:

R  » |y- y| < e.                                                                           (6.16)

Приближенным решением будут значения y, i = 0, 1, …, n.

Пример 6.2.

Применим первый  модифицированный метод Эйлера для решения  задачи Коши

y' (t) = y – ,  y(0) = 1,

рассмотренной ранее в примере 6.1.

Возьмем шаг h = 0.2. Тогда   n =  = 5.

В соответствии с (6.3) получим расчетную формулу первого модифицированного метода Эйлера:

yi+1 = yi + h f =  yi + 0.2 f, где

f = f(t, y) = y – ,

t = ti +  =  ti + 0.1,

y =  yi +f(ti, yi) = yi +0.1,

t0  = 0, y0 = 1, i = 0, 1, …, 4.

Решение представим в виде таблицы 6.3:

Таблица 6.3

i	ti	yi	f(ti, yi)	t	y	h f
0 1 2 3 4 5	0 0.2 0.4 0.6 0.8 1.0	1 1.1836 1.3426 1.4850 1.6152 1.7362	0.1 0.0850 0.0747 0.0677 0.0625	0.1 0.3 0.5 0.7 0.9	1.1 1.2682 1.4173 1.5527 1.6777	0.1836 0.1590 0.1424 0.1302 0.1210

Третий столбец таблицы 6.3 содержит приближенное решение yi, i = 0, 1, …, 5.

Сравним полученное приближенное решение с точным решением (6.11), представленном в таблице 6.2. Виднм, что погрешность составляет R = | y(ti) –  yi| = 0.0042.

Пример 6.3.

Применим второй  модифицированный метод Эйлера – Коши  для решения  задачи Коши

y' (t) = y – , y(0) = 1,

рассмотренной ранее в примерах 6.1 и 6.2. Так же, как и ранее, зададим шаг h = 0.2. Тогда   n =  = 5.

В соответствии с (6.14) получим расчетную формулу метода Эйлера – Коши:

yi+1 = yi + [f(ti, yi) + f(ti+1, )] = yi + 0.1[f(ti, yi) + f(ti+1, )],

где

f(ti, yi) = yi –

 = yi + h f(ti, yi) = yi + 0.1

t0  = 0, y0 = 1, i = 0, 1, …, 4.

Решение представим в виде таблицы 6.4:

Таблица 6.4

i	ti	yi	f(ti, yi)	ti+1		f(ti+1,)
0 1 2 3 4 5	0 0.2 0.4 0.6 0.8 1.0	1 1.1867 1.3484 1.4938 1.6272 1.7542	0.1 0.0850 0.0755 0.0690 0.0645	0.2 0.4 0.6 0.8 1.0	1.2 1.3566 1.4993 1.6180 1.7569	0.867 0.767 0.699 0.651 0.618

Таблица 6.4 заполняется последовательно по строкам, сначала первая строка, затем вторая и т. д. Третий столбец таблицы 6.4 содержит приближенное решение yi, i = 0, 1, …, 5.

Сравним полученное приближенное решение с точным решением (6.11), представленном в таблице 6.2. Видим, что погрешность составляет R = | y(ti) –  yi| = 0.0222.

6.4 Метод Рунге – Кутта

Метод Рунге – Кутта является одним из наиболее употребительных методов высокой точности. Метод Эйлера можно рассматривать как простейший вариант метода Рунге – Кутта.

Рассмотрим задачу Коши для дифференциального уравнения

y' (t) = f(t, y(t))

с начальным условием y(t0 ) = y0.

Как и в методе Эйлера, выберем шаг h =    и построим  сетку с системой узлов ti = t0 + ih, i = 0, 1, …, n.

Обозначим через yi приближенное значение искомого решения в точке ti.

Приведем  расчетные формулы метода Рунге – Кутта четвертого порядка точности:

yi+1 = yi + h(k+ 2k+ 2k + k),

k = f(ti, yi),

k = f(ti + , yi + k),                                                                          (6.17)

k = f(ti + , yi + k),

k = f(ti +h, yi + hk),

i = 0, 1, …, n.

Оценка погрешности. Оценка погрешности метода Рунге –  Кутта затруднительна. Грубую оценку погрешности дает правило Рунге (см. раздел 6.2). Так как метод Рунге - Кутта имеет четвертый порядок точности, т. е. p = 4, то оценка погрешности (6.6) примет вид

R  » |y- y|.                                                                                   (6.18)

Используя правило Рунге, можно построить процедуру приближенного вычисления решения задачи Коши методом Рунге – Кутта четвертого порядка точности с заданной точностью e. Нужно, начав вычисления с некоторого значения шага h, последовательно уменьшать это значение в два раза, каждый раз вычисляя приближенное значение y, i = 0, 1, …, n. Вычисления прекращаются тогда, когда  будет выполнено условие:

R  » |y- y| < e.                                                                           (6.19)

Приближенным решением будут значения y, i = 0, 1, …, n.

Пример 6.4.

Методом Рунге – Кутта четвертого порядка точности найдем решение на отрезке [0, 1] следующей задачи Коши.

y' (t) = 2ty,  y(0) = 1.                                                                               (6.20)

Возьмем шаг h = 0.1. Тогда   n =  = 10.

В соответствии с (6.17) расчетные формулы примут вид:

yi+1 = yi + h(k+ 2k+ 2k + k),

k = 2tiyi,

k = 2(ti + )(yi + k),                                                                          (6.21)

k = 2(ti + )(yi + k),

k = 2(ti +h)(yi + hk),

i = 0, 1, …, 10.

Задача (6.20) имеет точное решение: y(t) = e, поэтому погрешность определяется как абсолютная величина разности между точными и приближенными значениями ei = | y(ti) –  yi|.

Найденные по формулам (6.21) приближенные значения решения yi  и их погрешности ei представлены  в таблице 6.5:

Таблица 6.5

ti	yi	ei	ti	yi	ei
0.1 0.2 0.3 0.4 0.5	1.01005 1.04081 1.09417 1.17351 1.28403	10-9 4×10-9 2×10-8 6×10-8 2×10-7	0.6 0.7 0.8 0.9 1.0	1.43333 1.63232 1.89648 2.24790 2.71827	5×10-7 2×10-6 3×10-6 6×10-6 2×10-5

Задачи к зачету по курсу “Вычислительные методы”

Указание. Каждый студент вначале должен определить параметр своего контрольного задания, s = log10(1 +), где k - номер студента в списке группы, k = 1, 2, …  Решение задач должно быть оформлено аккуратно и содержать все промежуточные расчеты. В качестве образца можно взять примеры, рассмотренные в соответствующих разделах методических указаний.

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

4(1 –  x2)  –  ex = s  с точностью e  = 10-3.

2. Методом Зейделя решить систему уравнений с точностью e =10-3.

           6.2+s     2.2+s     1.2+s               16.55+s

A =     2.2+s     5.5+s    -1.5+s   ,  b =    10.55+s     .

          1.2+s    -1.5+s     7.2 +s                16.80+s

3. Найти приближение функции f(x) = esx  на отрезке [0, 1] многочленом Тейлора с точностью e = 10-3 . Вычислить es.

4. Вычислить приближенно по формуле средних прямоугольников интеграл  при n = 4 и оценить погрешность результата.

5. Методом Эйлера найти численное решение задачи Коши

y' = 2sy; y(0) = 1,  на отрезке [0, 1] с шагом h = 0.2.

Сравнить с точным решением.

Указания к выполнению лабораторных работ

Программой курса предусмотрено проведение четырех лабораторных работ. Лабораторные работы ориентированы на использование системы Maple.

Система Maple V была создана группой символьных вычислений в 1980 году в университете Waterloo, Канада. В конце 1997 года вышла реализация Maple V R5.

Maple V принадлежит к классу прикладных программных пакетов, объединенных под общим названием Computer Algebra Systems (CAS) - системы компьютерной алгебры. Самым важным отличием Maple от таких пакетов как MathCad, MatLAB, Mathematica, является то, что она была изначально задумана как символьный пакет. Как и любой представитель данного семейства продуктов, Maple ориентирована на решение широкого ряда математических проблем. Она включает в себя большое количество специальных пакетов для решения задач линейной и тензорной алгебры, евклидовой и аналитической геометрии, теории чисел, теории графов, теории вероятностей, математической статистики, комбинаторики, теории групп, численной аппроксимации и линейной оптимизации, задач финансовой математики и многих других.

В основу Maple положен алгоритмический язык высокого уровня, предназначенный для реализации обычного процедурного программирования. Maple-язык "понимает" все стандартные объекты типа циклов (while, for), операторов условного перехода (if-then-else), массивов (array), списков (list), наборов (set), таблиц и т.д. Есть также возможность работы с файлами, что позволяет строить системы, состоящие из множества модулей, подгружая необходимые процедуры в процессе выполнения программы, а также реализовывать ввод и вывод больших объемов данных. Реализованы также все стандартные процедуры обработки строковой информации.

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

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

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

Лабораторная работа №1.

Решение нелинейных уравнений и систем линейных уравнений.

Используемые функции: solve, fsolve,  plot.

1. Найти точное решение уравнения:5x2+2x – n = 0.

2. Найти приближенное решение этого же уравнения.

3. Построить график левой части уравнения.

4. Найти приближенное решение уравнения x2ex – n = 0.

5. Построить график левой части уравнения.

6. Найти точное решение системы уравнений.

2x1 + 6x2   – x3 = –12 + n

5x1 –  x2 + 2x3 =   29 + n

–3x1 – 4x2 +  x3 =     5 + n

7. Найти приближенное решение этой же системы уравнений.

Лабораторная работа №2.

Построение интерполяционных многочленов.

Используемые функции: interp, plot, subs.        

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

x      1          3         5         7         9    

y    0+n     4+n     2+n     6+n     8+n

2. Построить график полученного интерполяционного многочлена .

3. Найти значение функции в точке x = 6.

Лабораторная работа №3

Вычисление определенных интегралов.

Используемые функции: int, plot, evalf.

1. Найти аналитическое выражение для неопределенного интеграла .

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

3. Вычислить значение этого интеграла в пределах от 2 до n + 2:

4. Вычислить приближенное значение интеграла .

Лабораторная работа №4

Решение обыкновенных дифференциальных уравнений.

Используемые функции: dsolve, plot, odeplot, op, with.

1. Найти аналитическое решение задачи Коши: y'(t) = (1/n)(t + y), y(0) = n.

2. Построить график найденного решения на отрезке [0, n].

3. Найти численное решение задачи Коши y'(t) = sin(ny(t))+t2),  y(0) = n  в точках t = 1 и t = 2.

4. Построить график найденного решений на отрезке [0, 5].

Указания к выполнению курсовых работ

Цель курсовой работы – приобретение студентами практического опыта реализации на ЭВМ алгоритмов численных методов  для конкретных задач. Язык программирования выбирает студент.

Требования к выполнению курсовой работы

Результаты курсовой работы оформляются в виде отчета. Отчет по курсовой работе должен содержать следующие разделы:

1. Постановка задачи.

2. Описание математического метода.

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

4. Листинг программы.

5. Контрольный пример. Анализ полученных результатов.

Темы курсовых работ

Решение нелинейных уравнений

Указание. В курсовых работах 1 – 10 необходимо проанализировать два предложенных метода решения нелинейных уравнений, написать алгоритмы и программы этих методов. С помощью этих программ решить контрольный пример, предварительно локализовав корни уравнения (п. 2.2). Дать сравнительный анализ полученных результатов.

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

Контрольный пример. Найти один действительный корень уравнения x5 – x –  1 = 0 с точностью e  = 10-5.

Указание. При применении метода простых итераций преобразовать исходное уравнение так, чтобы итерационный процесс сходился (п. 2.4).

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

Контрольный пример. Найти три корня уравнения x3 – 4x2 + 2 = 0  с точностью e  = 10-5.

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

Контрольный пример. Найти три корня уравнения x3 + 3x2 – 1 = 0  с точностью e  = 10-5.

4. Решение нелинейных уравнений методом деления отрезка пополам и методом ложного положения.

Контрольный пример. Найти три корня уравнения x3 + 3x2 – 1 = 0  с точностью e  = 10-5.

5. Решение нелинейных уравнений методом  простых итераций и методом Ньютона.

Контрольный пример. Найти один действительный корень уравнения  x = 0.5 с точностью e  = 10-5.

6. Решение нелинейных уравнений методом  простых итераций и методом секущих.

Контрольный пример. Найти один действительный корень уравнения  x = 0.5 с точностью e  = 10-5.

7. Решение нелинейных уравнений методом  простых итераций и методом ложного положения.

Контрольный пример. Найти один действительный корень уравнения  x = 0.5 с точностью e  = 10-5.

8. Решение нелинейных уравнений методом секущих и методом Ньютона.

Контрольный пример. Найти три корня уравнения x3 + 3x2 – 3 = 0 с точностью e  = 10-5.

9. Решение нелинейных уравнений методом Ньютона и методом ложного положения.

Контрольный пример. Найти три корня уравнения x3 + x2 – 10x +8 = 0 с точностью e  = 10-5.

10. Решение нелинейных уравнений методом секущих и методом ложного положения.

Контрольный пример. Найти три корня уравнения x3 – x2 – 4x +4 = 0 с точностью e  = 10-5.

Решение систем линейных алгебраических уравнений

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

Контрольный пример. Решить систему уравнений

2.1x1  – 4.5x2  –  2.0x3   =  19.07

3.0x1 + 2.5x2 +  4.3x3  =  3.21

–6.0x1  + 3.5x2 + 2.5x3  =  –18.25

12. Решение системы линейных алгебраических уравнений  методом исключения Гаусса с выбором главного элемента по столбцу

Контрольный пример. Решить систему уравнений

1.00x1 + 0.42x2  +  0.54x3 +  0.66x4  =  0.3

0.42x1 + 1.00x2 +  0.32x3  +  0.44x4 =  0.5

0.54x1  + 0.32x2 + 1.00x3 +  0.22x4  =  0.7

0.66x1 + 0.22x2 + 1.00x3  –  1.0x4  =  0.9

13. Решение системы линейных алгебраических уравнений методом простых итераций Якоби.

Контрольный пример. Решить систему уравнений с точностью e  = 10-5.

–3.0x1  + 0.5x2  + 0.5x3 =  –56.65

0.5x1  – 6.0x2 +  0.5x3 = –160

0.5x1  + 0.5x2  – 3.0x3  =  –210

14. Решение системы линейных алгебраических уравнений методом Зейделя.

Контрольный пример. Решить систему уравнений с точностью e  = 10-5.

10x1  + 2x2  + x3 =  10

x1 + 10x2 +  2x3 = 12

x1  + x2  + 10x3  =  8

15. Вычисление определителя методом исключения Гаусса.

Контрольный пример. Вычислить определитель

det A = 3.0    1.5    0.1    1.0

0.4    0.5    4.0    6.5

0.3    1.2    3.0    0.7

1.8    2.2    2.5    1.4

16. Вычисление обратной матрицы методом исключения Гаусса.

Контрольный пример. Вычислить обратную матрицу A-1 для матрицы

A =    6.4375    2.1849   –3.7474    1.8822

 2.1356    5.2101    1.5220   –1.1234

 –3.7362    1.4998    7.6421    1.2324

 1.8666   –1.1004    1.2460    8.3312

17. Интерполяция  функции многочленами Лагранжа.

Контрольный пример. Построить интерполяционный многочлен Лагранжа для функции y = eпо точкам, заданным таблицей

x	0.00	0.25	0.50	0.75	1.00
e	1.0000000	0.9394131	0.7788008	0.7389685	0.3678794

Оценить погрешность интерполяции на отрезке [0, 1]. Вычислить y(0.4) и y(0.8).

18. Метод наименьших квадратов. Линейная и квадратичная аппроксимация

Численное интегрирование функций одной переменной

Указание. В курсовых работах 19 – 22 необходимо проанализировать  предложенные методы численного интегрирования функций одной переменной, написать алгоритмы и программы этих методов. С помощью этих программ решить контрольный пример. Проконтролировать погрешность, использовав правило Рунге (п. 5.5). Если можно, вычислить точное значение интеграла. Дать сравнительный анализ полученных результатов.

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

Контрольный пример. Вычислить ,  n = 10.

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

Контрольный пример. Вычислить ,  n = 10.

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

Контрольный пример. Вычислить   , n = 10.

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

Контрольный пример. Вычислить ,  n = 10.

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

Указание. В курсовых работах 23 – 26 необходимо проанализировать  предложенные методы численного решения задачи Коши, написать алгоритмы и программы этих методов. С помощью этих программ решить контрольный пример. Проконтролировать погрешность, использовав правило Рунге (пп. 6.2, 6.3, 6.4). Найти точное решение. Дать сравнительный анализ полученных результатов.

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

Контрольный пример. Найти численное решение задачи Коши

 

y' = y3,  y(0) = 0.5

на отрезке [0, 2] с шагом h = 0.2.

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

Контрольный пример. Найти численное решение задачи Коши

 

y' = t2,  y(0) = 1

на отрезке [0, 2] с шагом h = 0.2.

25. Решение задачи Коши для обыкновенных дифференциальных уравнений первым модифицированным  методом Эйлера и вторым  модифицированный метод Эйлера – Коши.

Контрольный пример. Найти численное решение задачи Коши

 

y' = sint,  y(0) = 1

на отрезке [0, 2] с шагом h = 0.2.

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

Контрольный пример. Найти численное решение задачи Коши

 

y' = 2cost,  y(0) = 0.

на отрезке [0, 2] с шагом h = 0.2.

Краткие сведения о математиках

Гаусс Карл Фридрих (1777 – 1855) – немецкий математик и физик, работы которого оказали большое влияние на развитие высшей алгебры, геометрии, теории чисел, теории электричества и магнетизма.

Зейдель Людвиг (1821 – 1896) – немецкий астроном и математик.

Коши Огюстен Луи (1789 – 1857) – французский математик, один из создателей современного математического анализа, теории дифференциальных уравнений и др.

Крамер Габриэль (1704 – 1752) – швейцарский математик.

Кутта В. М. (1867 – 1944) – немецкий математик.

Лагранж Жозеф Луи (1736 – 1813) – французский математик, механик и астроном. Один из создателей математического анализа, вариационного исчисления, классической аналитической механики.

Липшиц Рудольф (1832 – 1903) – немецкий математик.

Лейбниц Готфрид Вильгельм (1646 – 1716) – немецкий математик, физик и философ. Один из создателей дифференциального и интегрального исчислений.

Ньютон Исаак (1643 – 1727) – английский физик, механик, астроном, заложивший основы современного естествознания.

Рунге Карл Давид Тольме (1856 – 1927) – немецкий физик и математик.

Симпсон Томас (1710 – 1761) – английский математик.

Тейлор Брук (1685 – 1731) – английский математик и философ. Широко известная формула разложения функции в степенной ряд была получена им в 1712 г.

Эйлер Леонард (1707 – 1783) – математик, физик, механик, астроном. Родился в Швейцарии, с 1726 по 1741 г. и с 1776 по 1783 г. работал в России.

Якоби Карл Густав Якоб (1804 – 1851) – немецкий математик.

Список литературы

1. Амосов А., Дубинский Ю. А., Копченова Н.В. Вычислительные методы для инженеров: Учеб. пособие. – М.: Высш. шк., 1994.

2. Бахвалов Н. С. Численные методы. –  М.: Наука, 1973.

3. Волков Е. А. Численные методы. – М.: Наука, 1987.

4. Дьяконов В. П. Математическая система Maple V R3/R4/R5. – М.: Изд-во "СОЛОН", 1998.

5. Калиткин Н.Н. Численные методы. – М.: Наука, 1978.

6. Копченова Н.В., Марон И. А. Вычислительная математика в примерах и задачах. – М.: Наука, 1972.

7. Пирумов У.Г. Численные методы.: Учебное пособие. – М.: Изд-во МАИ, 1998.