Модификация алгоритма численного дифференцирования

Использование центральной разности (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;

Результаты вычислений сведем в таблицу:

i h q[i,1] q[i,2] q[i,3] q[i,4]

1 0.10000 2.7228145640

2 0.05000 2.7194145875 2.7182812620

3 0.02500 2.7185649917 2.7182817932 2.7182818286

4 0.01250 2.7183526177 2.7182818263 2.7182818285 2.7182818285

Точность результатов, приведенных на диагонали таблицы, не может не радовать пользователя.

Не нашли, что искали? Воспользуйтесь поиском: