Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н

Из (VI.24) видно, что основной член погрешности на каждом

a)/N , то полная погрешность метода левых прямоугольников

пропорциональна 0( N 2 ) .

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

вычисляется по формуле f (xi ) f (xi 1) x / 2, то погрешность

Заменив в (VI.25) первый член выражением (VI.23), значение функции в точке xi 1 будет разложением в ряд Тэйлора:

f (xi 1) f (xi ) f '(xi ) x 12 f ''(xi )( x)2

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

погрешность формулы трапеций на отрезке [а,b] по порядку величины равна 0(N 2 ).

148

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

пропорциональный f IV (x)( x)4 . Следовательно, погрешность формулы Симпсона на отрезке xi , xi 1 пропорциональна f IV (xi )( x)5 . , а полная погрешность на отрезке [а,b] по

порядку величины составляет 0(N 4 ).

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

Подставив (VI.26) в (VI.27), выполнив интегрирование и приведя подобные члены, получаем, что член, пропорциональный f (xi , yi ), сокращается, а интеграл от

(x xi )dx дает ( х)2/2. Интеграл от данного выражения по dy дает еще один множитель y. Аналогичный вклад дает

интеграл от члена, пропорционального ( y yi ) . Так как

Из (VI.28) видно, что погрешность интегрирования по одному параллелепипеду составляет 0(( х)3). Так как имеется N параллелепипедов, полная погрешность по порядку величины

прямоугольников i ~ 0(N-1).

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

равны 0(N-1) и 0(N-2 ) . Вообще можно показать, что если для одномерного случая погрешность составляет 0(N- ) , то в d- мерном случае она равна 0(N- /d ) .

6.5. Вычисление интегралов методом Монте-Карло

Проиллюстрируем идеи метода Монте-Карло на примере вычисления определенного интеграла от функции, зависящей от одной переменной. Пусть нам необходимо вычислить интеграл (VI.11) от некоторой заданной функции f(x) на интервале [а,b]. В предыдущем разделе мы рассмотрели несколько различных формул интегрирования, в которых использовались значения функции f(x), вычисляемые в

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

Рис. VI.7. К объяснению метода Монте-Карло

Представим себе прямоугольник высотой Н и длиной (b - а), такой, что функция f(x) целиком лежит внутри данного прямоугольника (рис. VI.7).

Сгенерируем N пар случайных чисел, равномерно распределенных в данном прямоугольнике:

yi f (xi ), является оценкой отношения интеграла от функции

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

N — полное количество точек, А — площадь прямоугольника. 151

Можно предложить и другой путь вычисления определенного интеграла, рассматривая его как среднее значение функции f(x) на отрезке [а,b]:

где xi — последовательность случайных чисел с равномерным

законом распределения на отрезке [а,b] .

Отметим, что в отличие от ранее упомянутых методов погрешность метода Монте - Карло не зависит от размерности

и меняется как 0(N-1/2 ) . Следовательно, для достаточно

больших d интегрирование по методу Монте - Карло будет приводить к меньшим погрешностям при тех же значениях N.

/ 2

Пример VI.6. Вычисление интеграла sin(x)dx методом

0

Монте - Карло в пакете MATLAB.

Решение:

%задание координат вершит прямоугольника

» X min = 0;

» X max = pi/2; » Y min = 0;

» Y max =1.5;

%генерация случайных координат

»N=2000;

»x =X min+ (X max - X min) *rand (N, 1) ;

»y =Y min+(Y max - Y min) *rand (N,l) ;

% подсчет числа точек, попавших под график функции

»s=0 ;

»for i =l:N

if y(i)<=feval(f,x(i)) s=s+l;

end;

end;

» s* (Xmax-Xmin) * (Ymax-Ymin) /N % вычисление значения

% интеграла ans =

1.0261 ;

% вычисление интеграла в соответствие c (VI.31)

»Fr=feval(f ,х) ;

»(Xmax-Xmin) /N*sum(Fr) ans =

1.0091

153

VII. МЕТОДЫ ОБРАБОТКИ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

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

7.1. Метод наименьших квадратов

Пусть в результате измерений в процессе опыта получена таблица некоторой зависимости f(x) (табл. VII.1).

Таблица VII.1

Исходные данные для метода наименьших квадратов

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

x1, x2,…, xn

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

154

y

yn .

. . . .

.

y1

x

y0

Рис. VII.1. К объяснению метода наименьших квадратов

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

(VII.1)

которая в точках x1, x2,…, xn

принимает значения как можно более близкие к табличным значениям y1, y2,…, yn.

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

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

Требование близости табличных значений y1, y2,…, yn и значений (VII.2) можно истолковать таким образом. Будем рассматривать совокупность значений функции f(x) из табл. VII.1 и совокупность значений (VII.2) как координаты двух точек n-мерного пространства. С учетом этого задача приближения функции может быть переформулирована следующим

метрикой Евклидова пространства, приходим к требованию, чтобы величина

была наименьшей. Это равносильно следующему: сумма квадратов

должна быть наименьшей.

Приведем окончательную формулировку задачи приближения функции f(x): для функции f(x), заданной табл. VII.1, найти функцию F(x) определенного вида так, чтобы сумма квадратов (VII.4) была наименьшей. Эта задача называется приближением функции методом наименьших квадратов. В качестве приближающих функций в зависимости от характера точечного графика f (x) часто используют функции, представленные далее. (Здесь а, b, т — неизвестные параметры.)

156

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

Рассмотрим метод нахождения параметров приближающей функции в общем виде на примере приближающей функции, зависящей от трех параметров:

Сумма квадратов разностей соответствующих значений функций f(x) и F(x) имеет вид:

Сумма является функцией Ф(а,b,с) трех переменных. Используя необходимое условие экстремума:

Ф 0, Ф 0, Ф 0 ,

а b с

получаем систему уравнений:

n

yi F xi , а,b,c Fn xi , а,b,с 0,

i 1

157