In my "Numerical Analysis" course, I learned the Romberg Algorithm to numerically calculate the integral.
The Romberg Algorithm as shown below:
$$T_{2n}(f)+\frac{1}{4^1-1}[T_{2n}(f)-T_{n}(f)]=S_n(f) \\ S_{2n}(f)+\frac{1}{4^2-1}[S_{2n}(f)-S_{n}(f)]=C_n(f) \\ C_{2n}(f)+\frac{1}{4^3-1}[C_{2n}(f)-C_{n}(f)]=R_n(f)$$
where $$T_n=\frac{h}{2}\left[ f(a)+ 2\sum_{i=1}^{n-1}f(x_i) +f(b)\right]$$ and $h=\frac{b-a}{n}$.
My code:
trapezium[func_, n_, {a_, b_}] :=
With[{h = (b - a)/n},
1/2 h (func@a + 2 Sum[func[a + i h], {i, 1, n - 1}] + func@b)
]
rombergCalc[func_,iter_, {a_, b_}] :=
Module[{m = 1},
Nest[
MovingAverage[#, {-1,4^(m++)}] &,
Table[trapezium[func, 2^i, {a, b}], {i, 0., iter}], 3]
]
The calculation process comes from my textbook
Fixed one bug 1
Update
Fixed bug 2
Test
rombergCalc[Exp, 5, {0, 1}]//InputForm
{1.7182818287945305, 1.7182818284603885, 1.7182818284590502}
My Question update
- In function
rombergCalc
, I utilized the usagem++
that I believe is not suitable in Mathematica Programming. Is there any other method to replacem++
or implement Romberg Algorithm?
- Why Block[{$MinPrecision = precision, $MaxPrecision = precision}..]
cannot give the result that contain significance digit
that I gave(seeing the graphic of textbook)?
(Thanks for @xzczd's solution for dealing with precision problem)
N[{a, b}, precision]
and replace
trapezium[func, 2^i, {a, b}], {i, 0., iter}]
with
trapezium[func, 2^i, {a, b}], {i, 0, iter}]
- Except for SetOptions[SelectedNotebook[], PrintPrecision -> 16]
or InputForm
, is there other solutions to set precision conveniently?
Block
and make$MinPrecision
equal to$MaxPrecision
. $\endgroup$