Лабораторная работа № 2
МЕТОДЫЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ ФУНКЦИЙ
Цель работы
Ознакомление с методами численного интегрирования, с понятием порядка точности численного метода, а также со способами контроля численных результатов.
Описание метода
Пусть необходимо вычислить интеграл
. (2.1)
Разобьем отрезок интегрирования на n частей. Введем в рассмотрение последовательность узловых точек xi Î[ a, b ]: xi=ih+a, i= 0 ,...,n. Величина 
называется шагом разбиения.
Все основные способы численного интегрирования сводятся к интерполяции функции по ее значениям в узловых точках f (xi) и интегрированию интерполяционного многочлена. При этом значение интеграла получается приближенно равным сумме
. (2.2)
При различном выборе Ai и xi получаются различные квадратурные формулы. Каждая из них обладает некоторой погрешностью e m, которую можно оценить следующим образом:
çe m ç£ chk, (2.3)
где c >0 - некоторая постоянная, не зависящая от h (зависящая от a, b, вида f (x) и метода интегрирования),
k -некоторое целое число, называемое порядком точности метода. Чем больше k, тем быстрее убывает погрешность при уменьшении h.
Предлагается рассмотреть квадратурные формулы Ньютона-Котеса, к которым, в частности, относятся формулы прямоугольников (левых, правых и симметричных), трапеций, парабол. В таблице 2.1 представлены эти формулы и значения констант для оценки погрешности по формуле (2.3). В таблице обозначено
, f ( j )(x) - j -ая производная f (x).
Таблица 2.1
| № | Название метода | Квадратурная формула | c | k |
| Левых прямоугольников | 
| 
| ||
| Правых прямоугольников | 
|
| ||
| Симметричных прямоугольников | 
| 
| ||
| Трапеций | 
| 
| ||
| Парабол | 
| 
|
Одним из способов практической оценки погрешности дискретизации, которая возникает при применении численных методов, (в том числе численного интегрирования, дифференцирования и т.п.) является правило Рунге, заключающееся в последовательном увеличении (например, удвоении) числа узловых точек n и соответствующем уменьшении шага дискретизации h. Оценка по правилу Рунге основана на предположении, что искомую величину J можно представить в виде
, (2.4)
где J – точное значение,
Jh – приближенный результат, полученный при шаге дискретизации равном h,
c – коэффициент, который предполагается независящим от h,
k – порядок точности метода,
d (h) - составляющая погрешности, которая считается пренебрежимо малой по сравнению с chk. В этом случае, уменьшив шаг дискретизации в два раза и отбросив d(h), нетрудно найти c и оценку погрешности
, (2.5)
где J - точное, а Jh, Jh/ 2 - приближенные значения интегралов, полученные, соответственно, с шагом h и h/ 2, k - порядок точности метода.
Тогда при заданной точности e величина h должна выбираться так, чтобы выполнялось условие
. (2.6)
Критерием допустимости отбрасывания малых величин можно считать стабильность величины K D
, (2.7)
полученной при уменьшении h в 4, 8 и т.д. раз:
Оценка погрешностей, связанных с машинным представлением чисел
Вычислительные ошибки этого типа порождаются ограниченной разрядностью представления чисел в ЭВМ. Эти ошибки резко возрастают в ситуациях, близких к математическим неопределенностям типа 0/0, ¥–¥, 0·¥.
Рассмотрим некоторое число A= 0.2354868971×10 p в машинном представлении с плавающей точкой
| Знак числа | Мантисса (M разрядов) | Знак порядка | Порядок | ||||||||
| ± | ± | p | |||||||||
Последние цифры (9,7,1), помещенные в нижней строке не умещаются в m разрядов и теряются. В худшем случае все потерянные цифры равны 9. Следовательно, предельная погрешность равна единице последнего разряда
.
Относительная предельная погрешность
. (2.8)
Здесь использовано правило записи числа в нормализованном виде: среди множества способов записи числа с плавающей точкой выбирается тот, при котором старшая значащая цифра располагается непосредственно за точкой (этим минимизируется объем памяти, необходимый для записи числа, так как не нужно хранить незначащие нули и позицию точки).
Отметим, что в машинном представлении используется двоичная система счисления, поэтому на самом деле
,
где M 2 - количество двоичных разрядов в мантиссе. Здесь мы используем десятичную систему только для удобства восприятия.
При сложении и вычитании двух чисел A ± B производится выравнивание порядков операндов по большему:
| A | + | + |
±
| B | + | - |
¯
| B¢ | + | + | ||||||||||
При этом последние разряды меньшего по порядку числа теряются и возникает погрешность, которая оценивается аналогично (2.8). Только необходимо помнить, что при оценке абсолютной погрешности число (2.8) нужно умножить на старшее по порядку число, участвующее в операции
. (2.9)
После выравнивания порядков производится операция, результат которой при сложении чисел одинакового знака может иметь мантиссу, превышающую единицу. При приведении числа к нормализованной форме производится сдвиг разрядов вправо. В результате возникает погрешность, которая может превысить (2.9). Поэтому общая погрешность операции оценивается следующим образом:
. (2.10)
Рассмотрим пример вычитания двух близких чисел:
| A | + | + |
-
| B | + | + |
=
| A - B | + | + |
Результат операции, преобразованный в нормализованную форму:
| A - B | + | - |
Пять нулей, записанные после цифр результата операции введены произвольно. Поскольку каждое число A и B могло быть усечено, то вместо нулей на самом деле могли бы стоять любые цифры, в том числе и девятки. Поэтому формула (2.9) дает реальную оценку и в этом случае.
Учитывая (2.10) найдем оценку относительной погрешности результата операции сложения и вычитания:
. (2.11)
Теперь рассмотрим квадратурные формулы типа (2.2):
(2.12)
(последнее равенство следует из того, что интеграл от функции f (x)=const должен вычисляться точно). Пусть
.
Тогда ошибка исходных данных (усечения значений функции) 
. Ошибка суммы приближенных значений
. (2.13)
При вычислении суммы накоплением возникает ситуация, описанная выше, когда складываются слагаемые разного порядка.
Оценка одного слагаемого суммы
.
Поэтому в соответствии с (2.9) ошибка округления при очередной операции сложения
.
а таких операций необходимо совершить n. Кроме того, в соответствии с (2.12) для приближенного вычисления интеграла сумму надо умножить на h. В связи с этим оценка погрешности округления
. (2.14)
Тогда, с учетом ошибок округления равенство (2.4) может принять вид
, (2.15)
причем последнее слагаемое обусловлено ограниченной разрядностью. Можно приближенно указать значения h и n, при которых оценка суммарной погрешности имеет минимальное значение. Для этого запишем общую оценку погрешности квадратурной формулы
,
и найдем минимум D(h):
,
,
,
.
Таким образом, можно считать, что
, (2.16)
где M - эквивалентное количество десятичных знаков мантиссы (при расчетах с обычной точностью M »7-8, с двойной точностью M »16).
Поскольку наличие значительной погрешности округления мешает использованию оценки (1.4.8), то при расчетах приходится ограничиваться меньшими n и большими h, чем это следует из (2.16). Кроме того, существуют различные способы, чтобы ограничить возрастание погрешности, связанное с математическими неопределенностями.
Возможность контроля погрешности округления несколько облегчает то обстоятельство, что эта погрешность, в отличие от остальных типов погрешностей, как правило, ведет себя достаточно хаотично, и по уровню этой хаотической составляющей можно судить, хотя и очень приближенно, о ее величине.
Пример
В приведенных ниже таблицах показаны результаты численного интегрирования функции f(x)= 6 x 5 на интервале [0,1] методом парабол (точное значение интеграла равно 1). Величины K D и DРунге получены по формулам (2.7) и (2.6), Dтеор – по (2.3) с учетом данных табл. 2.1, Dточное равно разности между точным и приближенным значением. Результаты, приведенные в таблице 2.2, получены путем вычисления с двойной точностью (мантисса 16 десятичных знаков), в таблице 2.3 – с одинарной точностью (мантисса 7-8 знаков). Из таблиц видно, что в данном случае коэффициент уменьшения погрешности K D весьма стабилен до значений n примерно равных n 0 (2.16). Кроме того, видно, что при этих n оценка по Рунге DРунге практически совпадает с Dточное, в то время как оценка через производную D теор превышает их в два раза.
Таблица 2.2
| N | K D | Dточное | DРунге | Dтеор |
| – | -1.2500×10-1 | – | 2.5000×10-1 | |
| – | -7.8125×10-3 | -7.8125×10-3 | 1.5625×10-2 | |
| 16.0 | -4.8828×10-4 | -4.8828×10-4 | 9.7656×10-4 | |
| 16.0 | -3.0518×10-5 | -3.0518×10-5 | 6.1035×10-5 | |
| 16.0 | -1.9073×10-6 | -1.9073×10-6 | 3.8147×10-6 | |
| 16.0 | -1.1921×10-7 | -1.1921×10-7 | 2.3842×10-7 | |
| 16.0 | -7.4506×10-9 | -7.4506×10-9 | 1.4901×10-8 | |
| 16.0 | -4.6566×10-10 | -4.6566×10-10 | 9.3132×10-10 | |
| 16.0 | -2.9104×10-11 | -2.9104×10-11 | 5.8208×10-11 | |
| 16.0 | -1.8190×10-12 | -1.8190×10-12 | 3.6380×10-12 | |
| 16.0 | -1.1366×10-13 | -1.1369×10-13 | 2.2737×10-13 | |
| 16.1 | -7.2997×10-15 | -7.0906×10-15 | 1.4211×10-14 | |
| 13.9 | 3.6082×10-16 | -5.0823×10-16 | 8.8818×10-16 | |
| -68.7 | 2.4980×10-16 | 7.4015×10-18 | 5.5511×10-17 | |
| 0.2 | -1.9429×10-16 | 3.2078×10-17 | 3.4694×10-18 | |
| 2.6 | -4.1633×10-16 | 1.2338×10-17 | 2.1684×10-19 | |
| 0.1 | -1.7486×10-15 | 8.6353×10-17 | 1.3553×10-20 |
Таблица 2.3
| N | K D | Dточное | DРунге | Dтеор |
| – | -1.2500×10-1 | – | 2.5000×10-1 | |
| – | -7.8125×10-3 | -7.8125×10-3 | 1.5625×10-2 | |
| 16.0 | -4.8828×10-4 | -4.8828×10-4 | 9.7656×10-4 | |
| 16.0 | -3.0518×10-5 | -3.0518×10-5 | 6.1035×10-5 | |
| 16.0 | -1.9073×10-6 | -1.9073×10-6 | 3.8147×10-6 | |
| 16.0 | -1.1921×10-7 | -1.1921×10-7 | 2.3842×10-7 | |
| 45.0 | -1.1921×10-7 | -2.6491×10-9 | 1.4901×10-8 | |
| 0.2 | 1.1921×10-7 | -1.4570×10-8 | 9.3132×10-10 | |
| 2.2 | 2.3842×10-7 | -6.6227×10-9 | 5.8208×10-11 |