Comparison of 1st and 2nd order¶
using FundamentalsNumericalComputation
f = x -> sin( exp(x+1) );
exact_value = cos(exp(1))*exp(1);
We’ll run both formulas in parallel for a sequence of \(h\) values.
h = @. 4. ^(-1:-1:-8)
FD1 = []; FD2 = [];
for h in h
push!(FD1, (f(h)-f(0)) / h )
push!(FD2, (f(h)-f(-h)) / 2h )
end
pretty_table([h FD1 FD2],["h","FD1","FD2"],backend=:html)
| h | FD1 | FD2 |
|---|---|---|
| 0.25 | -3.0100196002533077 | -2.3924525328620603 |
| 0.0625 | -2.6451014322652484 | -2.473905902362808 |
| 0.015625 | -2.521133723616341 | -2.478075703599398 |
| 0.00390625 | -2.4891011707672988 | -2.478332620635271 |
| 0.0009765625 | -2.4810408642263155 | -2.478348663491829 |
| 0.000244140625 | -2.4790227172861705 | -2.478349666114127 |
| 6.103515625e-5 | -2.478517991587978 | -2.4783497287785394 |
| 1.52587890625e-5 | -2.4783917983913852 | -2.4783497326825454 |
All that’s easy to see from this table is that FD2 appears to converge to the same result as FD1, but more rapidly. In each case \(h\) is decreased by a factor of 4, so that the error is reduced by a factor of 4 in the first-order method and 16 in the second-order method.
error_FD1 = @. exact_value - FD1
error_FD2 = @. exact_value - FD2
table = (n=Int.(log2.(h)),fd1=error_FD1,fd2=error_FD2)
pretty_table(table,["log_2(h)","error in FD1","error in FD2"],backend=:html)
| log_2(h) | error in FD1 | error in FD2 |
|---|---|---|
| -2 | 0.5316698672980729 | -0.08589720009317459 |
| -4 | 0.16675169931001355 | -0.004443830592427034 |
| -6 | 0.04278399066110605 | -0.00027402935583697996 |
| -8 | 0.010751437812063891 | -1.711231996370799e-5 |
| -10 | 0.002691131271080671 | -1.0694634058339147e-6 |
| -12 | 0.0006729843309356554 | -6.684110775978525e-8 |
| -14 | 0.00016825863274316788 | -4.176695433955047e-9 |
| -16 | 4.2065436150373614e-5 | -2.7268942659475215e-10 |
A graphical comparison can be clearer. On a log-log scale, the error should (roughly) be a straight line whose slope is the order of accuracy. However, it’s conventional in convergence plots, to show \(h\) decreasing from left to right, which negates the slopes.
plot(h,abs.([error_FD1 error_FD2]),m=:o,label=["error FD1" "error FD2"],
xflip=true, xaxis=(:log10,"\$h\$"), yaxis=(:log10,"error in \$f'(0)\$"),
title="Convergence of finite differences", leg=:bottomleft)
plot!(h,[h h.^2],l=:dash,label=["order1" "order2"]) # perfect 1st and 2nd order