7.4 Приближенное дифференцирование, основанное на

ин терполяционной формуле Ньютона

Идея численного дифференцирования заключается в том, что функцию y(x), заданную в

равноотстоящих точках xi (i = 0, 1, ..., n) отрезка [a, b] с помощью значений yi=f(xi)

приближенно заменяют интерполяционным полиномом Ньютона, построенном для системы

узлов x0, x1, ..., xk (kn), и вычисляют производные y'=f'(x), y''=f''(x) и т.д.

y '  x0=1

h y0−

2 y0

2 

3 y0

3 −

4 y0

4 

5 y0

5 −... , (7.2)

На практике приближенное дифференцирование применяют в основном для функций

заданных в виде таблицы.

В Scilab численное дифференцирование реализовано командой dy=diff(y [,n ]),

где y _ значения функции y(x) в виде вектора вещественных чисел, n _ порядок

дифференцирования. Результат работы функции - вектор вещественных чисел dy,

представляющий собой разности порядка n интерполяционного полином Ньютона y, 2y, ...,

ky. Рассмотрим работу функции на примере.

ЗАДАЧА 7.5

Найти y'(50) функции y=lg(x), заданной в виде таблицы.

Решение данной задачи с комментариями представлено в листинге 7.7.

-->h=5;x=50:5:65;

-->y=log10(x)

y = 1.69897 1.7403627 1.7781513 1.8129134

-->dy=diff(y)

dy = 0.0413927 0.0377886 0.0347621

-->dy2=diff(y,2)

dy2 = - 0.0036041 - 0.0030265

-->dy3=diff(y,3)

dy3 = 0.0005777

-->//Приближенное значение y'(50) по формуле (7.2)

-->Y=(dy(1)-dy2(1)/2+dy3(1)/3)/h

Y = 0.0086775

-->//Значение y'(50) для lg'(x)=1/ln(10)/x

-->1/log(10)/x(1)

ans = 0.0086859

-->//Приближенное значение y'(x), x=50,55,60 (7.2)

-->Y=(dy-dy2(1:$-1)/2+dy3(1:$-2)/3)/h

Y = 0.0086389 0.0079181 0.0073128

-->//Значение y'(x), x=50,55,60, для lg'(x)=1/ln(10)/x

-->(1/log(10))./x(1:$-1)

ans = 0.0086859 0.0078963 0.0072382

Листинг 7.7

7.5 Вычисление производной функции в т очке.

Приближенное вычисление час тных производных.

Более универсальной командой дифференцирования является команда

g=numdiff(fun,x)

здесь fun - имя функции, задающей выражение для дифференцирования. Функция должна

быть задана в виде y=fun(x [, p1, p2,..., pn]), где x - переменнаяпо которой

будет проводится дифференцирование. Если параметры p1, p2,..., pn присутствуют в

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

g=numdiff(list(fun,p1,p2,..pn),x).

Результат работы функции _ матрица gij=

df i

dx j

.

Рассмотрим несколько примеров.

ЗАДАЧА 7.6

Вычислить f'(1), если f(x)=(x+2)3+5x.

Листинг 7.8 содержит решение данной задачи.

-->function f=my(x), f=(x+2)^3+5*x, endfunction;

-->numdiff(my,1)

ans = 32.

-->x=1;3*(x+2)^2+5

ans = 32.

Листинг 7.8

5

ЗАДАЧА 7.7

Вычислить f'(x) в точках 0, 1, 2, 3 для f(x)=(x+2)3+5x.

Решение:

-->v=0:3;

-->numdiff(my,v)

ans =

17. 0. 0. 0.

0. 32. 0. 0.

0. 0. 52.999999 0.

0. 0. 0. 80.000002

-->function f1=my1(x), f1=3*(x+2)^2+5, endfunction;

-->my1(v)

ans = 17. 32. 53. 80.

Листинг 7.9

ЗАДАЧА 7.8

Задана функция многих переменных y  x1, x2, x3=x1 x2

x3x1 2 x3 . Вычислить

dy

dx1

, dy

dx2

, dy

dx3

в точке (1, 2, 3).

Решение задачи представлено в листинге 7.10.

-->function [Y]=f(X), Y=X(1)*X(2)^X(3)+X(1)^2*X(3),endfunction

-->X=[1 2 3];

-->numdiff(f,X)

ans = 14. 12. 6.5451775

-->//--------------------------------

-->function [Y]=f1(X),

Y(1)=X(2)^X(3)+2*X(1)*X(3),

Y(2)=X(1)*X(3)*X(2)^(X(3)-1),

Y(3)=x(1)*X(2)^X(3)*log(X(2))+X(1)^2,

endfunction

-->f1(X)

ans =

14.

12.

6.5451774

Листинг-->function f=my(x), f=(x+2)^3+5*x, endfunction;

-->numdiff(my,1)

ans = 32.

-->x=1;3*(x+2)^2+5

ans = 32.

Листинг 7.8

5

ЗАДАЧА 7.7

Вычислить f'(x) в точках 0, 1, 2, 3 для f(x)=(x+2)3+5x.

Решение:

-->v=0:3;

-->numdiff(my,v)

ans =

17. 0. 0. 0.

0. 32. 0. 0.

0. 0. 52.999999 0.

0. 0. 0. 80.000002

-->function f1=my1(x), f1=3*(x+2)^2+5, endfunction;

-->my1(v)

ans = 17. 32. 53. 80.

Листинг 7.9

ЗАДАЧА 7.8

Задана функция многих переменных y  x1, x2, x3=x1 x2

x3x1 2 x3 . Вычислить

dy

dx1

, dy

dx2

, dy

dx3

в точке (1, 2, 3).

Решение задачи представлено в листинге 7.10.

-->function [Y]=f(X), Y=X(1)*X(2)^X(3)+X(1)^2*X(3),endfunction

-->X=[1 2 3];

-->numdiff(f,X)

ans = 14. 12. 6.5451775

-->//--------------------------------

-->function [Y]=f1(X),

Y(1)=X(2)^X(3)+2*X(1)*X(3),

Y(2)=X(1)*X(3)*X(2)^(X(3)-1),

Y(3)=x(1)*X(2)^X(3)*log(X(2))+X(1)^2,

endfunction

-->f1(X)

ans =

14.

12.

6.5451774

Листинг