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

6. 2. Особенностизадачичисленногодифференцирования функций, заданныхтаблицей

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

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

Полагая, что погрешность интерполирования определяется формулой

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

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

6.3.Интегрирование функций, заданных аналитически

Сгеометрической точки зрения определенный интеграл

— есть площадь фигуры, ограниченная графиком функции f(x) и прямыми х = а , х = Ь (рис. VI.3).

138

y

f(x)

a b x

Рис. VI.3 К объяснению геометрического смысла определенного интеграла

Разделим отрезок [а,b] на N равных отрезков длиной x, где

Тогда координата правого конца i-го отрезка определяется по формуле

Рис. VI.4. Геометрическая интерпретация метода левых прямоугольников

В этом случае фигура, ограниченная графиком функции и прямымиx xi , x xi 1 , является трапецией. Искомый определенный интеграл определяется как сумма площадей всех трапеций:

Более высокая точность вычисления интегралов обеспечивается при использовании параболической интерполяции (полиномом второй степени) по трем соседним точкам:

Для нахождения коэффициентов а,b,с полинома, проходящего через точки (x0 , y0 ), (x1, y1), (x2 , y2 ), нужно найти решение

следующей системы линейных уравнений:

относительно неизвестных а ,b , с.

Решив систему (VI.18) относительно неизвестных а, b, с любым известным методом (например, Крамера), подставив найденные выражения в (VI.17) и выполнив элементарные преобразования, получаем

(VI.19)

Отметим, что формула (VI.19) совпадает с интерполяционным полиномом Лагранжа, подробно рассмотренным в главе V, при интерполяции по трем точкам.

141

Площадь под параболой у = у(х) на интервале [ x0 , x2 ]

где x x1 x0 x2 x1 .

Искомый определенный интеграл находится как площадь всех параболических сегментов (формула Симпсона):

Обратите внимание на то обстоятельство, что в формуле Симпсона N должно быть четным числом.

/ 2

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

0

прямоугольников в пакете MATLAB.

Решение:

»f=inline('sin(x)'); % задание подынтегральной функции

»Xmin=0;

»Xmax=pi/2;

»N=2001;

»i=l:N;

»dx=(Xmax-Xmin)/(N-l); % шаг интегрирования

»x=Xmin:dx:Xmax; % вычисление координат узлов сетки

»y=feval(f,x); % вычисление значений функции в узлах сетки % вычисление интеграла по формуле правых прямоугольников

»m=2:N;

»yl(m-l)=y(m) ;

»Fr=sum(yl) *dx

3.9284e-004

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

» m=l:N-l;

» yl (m) =у (m) ; » Fl=sum(yl)*dx Fl =

0.9996 » Fl-1 ans =

-3.9295e-004

%вычисление интеграла методом трапеций

» s=0;

for i=2:N-l s=s+y[i];

end;

Ft=(0.5*y(l)+s+0.5*y(N))*dx Ft =

1.0000

»Ft-1 ans =

-5.1456e-008

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

»s=0;

for i=2:N-l

if i-2*ceil(i/2)= =0 k=4;

else k=2;

end;

s=s+k*y(i);

end;

Fs=(y(l)+s+y(N))*dx/3 Fs =

1.0000

143

» Fs-1 ans =

1.5543e-015

Для вычисления значений определенных интегралов в пакете

MATLAB есть функции quad ( ), guadl( ), trapz( ), cumtrapz( ),

которые имеют следующий синтаксис:

q = quad (fun, a, b) — возвращает значение интеграла от функции fun на интервале [а, b], при вычислении используется адаптивный метод Симпсона;

q = quad(fun,a,b,tol) — возвращает значение интеграла от функции fun с заданной относительной погрешностью tol (по умолчанию tol=10-3);

q = quad (fun, a, b,tol, trace) — возвращает значение интеграла от функции fun на интервале [а, b] на каждом шаге итерационного процесса;

q = quad(fun,a,b, tol, trace, pi,p2,...) — возвращает значение интеграла от функции fun на на интервале [а, b] на каждом шаге итерационного процесса, pi, р2 — параметры, передаваемые в функцию fun;

[q,fcnt] = quadi(fun,a,b,...) — возвращает в переменную fcnt до-

полнительно к значению интеграла число выполненных итераций.

Функция guad1 ( ) возвращает значения интеграла, используя для вычислений метод Лоббато (Lobbato). Синтаксис этой функции аналогичен синтаксису функции quad ().

/ 2

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

0

использованием функций quard ().

Решение:

»q=quad('sin',0,pi/2,10^-4) q =

1.0000

»q-1

ans =

144

-3.7216e-008

»q-1 ans =

-2.1269e-009

»[q, fnct]=quad( 'sin' ,0, pi/2 , 10^-6, 'trace' ) ;

» 17

Функция trapz( ) вычисляет интеграл, используя метод трапеций. Синтаксис этой функции следующий:

z = trapz(Y) — возвращает значение определенного интеграла

интервале [х(1), x(N)];

z = trapz(Х,Y,dim) — интегрирует вектор Y, формируемый из чисел, расположенных в размерности dim многомерного массива.

Функция cumtrapz( ), вычисляет интеграл, как функцию с переменным верхним пределом. Синтаксис функции cumtrapz() аналогичен синтаксису функции trapz( ).

145

/ 2

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

0

использованием встроенной функции пакета MATLAB.

Решение:

»x =0:0.01:pi/2; % задание координат узловых точек

»y= sin(x); % вычисление значений подынтегральной функции

%в узловых точках

»trapz (у) % вычисление значения интеграла в предположении

%о том, что шаг интегрирования равен единице

ans = 99.9195

»trapz (х, у) % вычисление значения интеграла на отрезке

ans = 0.9992

Пример VI.5. Вычислить интеграл с переменным верхним

x

пределом sin( )d на интервале 0, 3 / 2 .

0

Решение:

»х=0:0.01:3*pi/2; % задание координат узловых точек

»y=sin(x); % вычисление значений подынтегральной функции

%в узловых точках

»z=cumtrapz(x,y); % вычисление значений интеграла с

%переменным верхним пределом в узловых точках

» plot(x,y,x,z) % построение графиков подынтегральной

%функции и интеграла с переменным

%верхним пределом (рис. VI.6)

146

x

Рис. VI.6. Графики функций f(x)=sin (x) (1), (x) sin( )d

0

6.4. Погрешность численного интегрирования

Для нахождения зависимостей погрешности вычисления определенного интеграла на отрезке [а, b] от числа отрезков

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

разности между точным значением интеграла и его оценкой f (xi ) x : 147