# 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