Nejlepší odpověď
Máte na mysli funkci, která hodnotí \ frac {\ sin x} {x} pro x \ ne 0 a 1 pro x = 0?
Je zřejmé, že by bylo snadné jen něco naprogramovat pomocí příkazu jako
function result = mySinc(x)
\% My implementation of the sinc function.
if x == 0
result = 1;
else
result = sin(x)/x;
end
end
Ale to by ve skutečnosti mohlo být méně než ideální. Jistě, vyhýbá se singularitě a má správnou matematickou strukturu, ale samozřejmě si musíme pamatovat, že hodnoty v x
nejsou čísla: jsou to vzory bitů reprezentované přepínači v hardwaru a pouze přibližné počty. Většinou fungují docela dobře, ale existují okolnosti, kdy věci prostě nefungují správně. Například při hodnocení a + log(1+x)
pro malé x
, nebo při hodnocení a-b
když a
a b
mají velmi podobné hodnoty. může přijít o hodně přesnosti, pokud ji naivně vyhodnotíte (i když někdy možná nemáte jinou možnost, bohužel). Dělení malým počtem může také způsobit problémy, i když naštěstí jsou \ sin x a x srovnatelné, když | x | \ ll 1. Možná by to nebyl problém, ale při práci s aritmetikou s plovoucí desetinnou čárkou je zřídka dobrým zvykem mít kód jako if x == foo
.
Co dělat v tomto případě? Vezměme si Taylorovu řadu \ frac {\ sin x} {x} se středem na x = 0. Nebudu to psát, ale určitě víte, jak to vypočítat.
Poté nahraďte test if x == 0
znakem if abs(x) < sqrt(eps)
, kde v Matlabu eps
vrací stroj epsilon (srovnatelný s 10 ^ {- 16} na většině systémů v těchto dnech) a nahradí řádek result = 1;
by result = ...;
, kde ...
bude váš výpočet Taylorovy řady zkrácené po příslušném pořadí (obvykle asi 2 ale měli byste sami zkontrolovat měřítko termínu nejvyššího řádu, který vyberete pro x \ sim 10 ^ {- 8} a zda bude výsledek tak přesný, jak potřebujete, nebo ne).
The pěkná věc je, že to, co ztratíte na rychlosti běhu tohoto kódu, pravděpodobně získáte přesnost, zvláště při hodnocení sinusového integrálu \ operatorname {Si} (x) = \ int\_0 ^ xf (t) dt, kde f (t ) je kontinuální rozšíření funkce \ operatorname {sinc} do všech realů jako její doména.