3. Численное дифференцирование с заданной точностью

Алгоритм численного дифференцирования с заданной точностью построим на последовательном делении в точке х текущего приращения аргумента =hпополам, и вычислении отношенияD(x,h)=:

1. зададим приращение аргумента h, его значение х, точность вычисления eps и вычислим начальное значение функции y=D(x,h);

2. сохраним старое значение функции ys, уменьшими приращение аргумента вдвое и вычислим новое значение y=D(x,h/2);

3. если abs(y-ys) > eps , то повторим процесс, возвратившись к пункту 2, а иначе – выведем результат.

Составим фрагмент программы:

consteps=1.0e-6;

type fun: function(x:real):real; {$F+}

function f(x:real):real;

begin f:=exp(x) end;

function Y1( x: real; ff: fun):real;

function D(x,h:real):real;

begin D:=0.5*(ff(x+h)-ff(x-h))/h end;

var y,ys:real; k,kmax: integer;

begin

h:=0.1; y:=D(x,h);

k:=1; kmax:=100;

yy:=ff(x);

repeat

Inc(k);

if ( k >kmax) then begin writeln(‘Error!’); Exit end

ys:=y; h:=0.5*h; y:=D(x,h);

until abs(y-ys) < eps;

Y1:=y

end;

Приведем результаты расчетов:

k h Y1 abs(Y1 – )

1 0.0500000000 2.71941458750 0.0011327590

2 0.0250000000 2.71856499170 0.0002831633

3 0.0125000000 2.71835261770 0.0000707892

4 0.0062500000 2.71829952570 0.0000176973

5 0.0031250000 2.71828625290 0.0000044245

6 0.0015625000 2.71828293450 0.0000011061

7 0.0007812500 2.71828293450 0.0000002760

8 0.0003906250 2.71828210450 0.0000000711

9 0.0001953125 2.71828189960 0.0000000199

10 0.0000976563 2.71828184840 0.0000000106

11 0.0000488281 2.71828185770 0.0000000292

12 0.0000244141 2.71828182040 0.0000000080

13 0.0000122070 2.71828174590 0.0000000825

14 0.0000061035 2.71828144790 0.0000003806

15 0.0000030518 2.71828234200 0.0000005135

16 0.0000015259 2.71828293800 0.0000011095

17 0.0000007629 2.71827936170 0.0000024667

18 0.0000003815 2.71827697750 0.0000048509

4. Модификация алгоритма численного дифференцирования Использование центральной разности (6.3) для приближения производной позволяет проводить вычисления с точность порядка :

(6.4)

где , а С – произвольная константа. Уменьшая шаг в два раза, получим

(6.5)

где . Из соотношений (6.4) и (6.5) нетрудно получить, что

(6.6)

Таким образом, не проводя лишних вычисления, можно значительно улучшить результат, используя комбинацию величин и. Формула (6.6) имеет порядок точности. Убедитесь в этом самостоятельно, вводя небольшую поправку в только что разработанную подпрограмму и проведя расчеты.

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

где . Разделив текущий интервал пополам, получим выражение для:

Подставляя эти соотношения в (6.6) нетрудно вывести следующую формулу

(6.7)

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

Описанную выше идею комбинирования можно обобщить, записав выражение:

, (6.8)

где коэффициент .

Составим подпрограмму, заменив переменную массивом q[i,j] :

function DifR( x: real; f:fun):real;

function D(x,h:real):real;

begin D:=0.5*(f(x+h)-f(x-h))/h end;

const n=10;

var q:array{1..n,1..n] of real; hm2:real; I,j:integer;

begin

h:=0.1; q[1,1]:=D(x,h); i=1;

repeat

h:=0.5*h; m2:=1;

q[i,1] :=D(x,h); Inc(i);

if (i >=n) then

begin writeln(‘Error!’); Exit end;

for j:=2 to i do

begin m2:=m2*4; q[i,j]:=(m2*q[i.j-1] - q[i-1.j-1])/(m2-1) end;

until abs(q[i.j] - q[i-1.j-1]) < eps;

DifR:= q[i.j]

end;