Попытка реализовать формулу разности в MATLAB - PullRequest
Новая
       39

Попытка реализовать формулу разности в MATLAB

0 голосов
/ 05 января 2019

Я пытаюсь реализовать формулу разности 

f'(x) ≈ [ f(x+h) - f(x) ] / h

с использованием MATLAB для x=1 и h=10^-k, где k=0,...,16. Кроме того, я хочу подготовить ошибку.

Ниже мой код. Я вижу, что ошибка составляет около 3, что я считаю слишком большой. Должно быть близко к 0. 

syms f(x)
f(x) = tan(x);
df = diff(f,x);
x = 1;
for k = 0:16
    h = 10^-k;
    finitediff = double((f(x+h)-f(x))/h);
    err = double(abs(finitediff-df(x)));
end

Ответы [ 2 ]

0 голосов
/ 05 января 2019

Вы вычисляете x+h численно. с x=1 наименьшее число, большее x, которое вы можете сделать, составляет x+eps(x) = 1+eps(1). eps(1) - это 2.2204e-16. Таким образом, добавление h=1e-16 не меняет x. Теперь вы символически вычисляете tan(1)-tan(1), что равно 0. Следовательно, ваше конечное разностное приближение к производной равно 0.

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

x = 1;
h = 1e-14;
f(x+h)

возвращает tan(1). То есть символический набор инструментов считает 1 + 1e-14 в контексте функции tan достаточно близким к 1, чтобы гарантировать округление. Как ни странно, 

f(x+h)-f(x)

возвращает 0, тогда как 

tan(x+h)-tan(x)

возвращает 3.4195e-14 (что близко к ожидаемому 3.4255e-14).

Обратите внимание, что, как видно из графика , Хантер строит графики , наилучшее приближение для h=10^-8, так как при уменьшении h ошибки округления в x+h начинают увеличиваться (используйте set(gca,'yscale','log') чтобы увидеть это).

0 голосов
/ 05 января 2019

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

  • ошибки, сгенерированные вычислением a/b, когда оба значения a и b очень и очень малы.

  • вычисление a-b, когда a и b очень очень близки к тому, что MATLAB даст 0.

Вот результат, когда k изменяется от 1 до 15:

enter image description here

Спасибо за проницательный комментарий @CrisLuengo!

, который показывает, что err мгновенно падает почти до нуля, и снова повышается, когда h становится 1e-9 (ситуация 1 происходит после этого). Наконец df становится 0, когда h становится 1e-14 (здесь происходит ситуация 2).

Я добавил несколько строк кода к вашей, чтобы показать это: 

clc;
clear;
format long
syms f(x)
f(x) = tan(x);
h=1;
df = diff(f,x);
double(df(1));
x=1;


range=1:15;
[finitediff,err]=deal(zeros(size(range)));
for k=range
    h=10^-k;
    finitediff(k)=double((f(x+h)-f(x))/h);
    err(k)=double(abs(finitediff(k)-df(1)));
end

figure(1)

subplot(1,2,1)
hold on
plot(err)
plot(err,'bx')
set(gca,'yscale','log')
title('err')

subplot(1,2,2)
hold on
ezplot(df)
axis([0.5 1.5 0 5])
plot(ones(size(range)),finitediff,'rx','MarkerSize',7)
for ii=range
  text(1,finitediff(ii),['h=1e-' num2str(ii)])
end
...