library(Ryacas)
In R
, the package Rmpfr
package provide
some functions for arbitrary-precision arithmetic. The way it is done is
to allow for using more memory for storing numbers. This is very good
for some use-cases, e.g. the example in Rmpfr
’s a vignette
“Accurately Computing \(\log(1 -
\exp(-|a|))\)”.
Here, we demonstrate how to do investigate other aspects of
arbitrary-precision arithmetic using Ryacas
, e.g. solving
systems of linear equations and matrix inverses.
First see the problem when using floating-point arithmetic:
0.3 - 0.2) - 0.1 (
## [1] -2.775558e-17
yac_str("(0.3 - 0.2) - 0.1") # decimal number does not
## [1] "0"
# always work well; often
# it is better to represent as
# rational number if possible
# (1/3 - 1/5) = 5/15 - 3/15 = 2/15
1/3 - 1/5) - 2/15 (
## [1] -2.775558e-17
yac_str("(1/3 - 1/5) - 2/15")
## [1] "0"
The yacas
function N()
gives the numeric (floating-point) value for a precision.
yac_str("1/3")
## [1] "1/3"
yac_str("N(1/3)")
## [1] "0.3333333333"
yac_str("N(1/3, 1)")
## [1] "0.3"
yac_str("N(1/3, 200)")
## [1] "0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333"
Consider the polynomial \((x-1)^7\) with root 1 (with multiplicity 7). We can consider it both in factorised form and in expanded form:
<- "(x-1)^7"
pol <- pol %>% yac_expr()
p1 p1
## expression((x - 1)^7)
<- pol %>% y_fn("Expand") %>% yac_expr()
p2 p2
## expression(x^7 - 7 * x^6 + 21 * x^5 - 35 * x^4 + 35 * x^3 - 21 *
## x^2 + 7 * x - 1)
Mathematically, these are the same objects, but when it comes to numerical computations, the picture is different:
eval(p1, list(x = 1.001))
## [1] 1e-21
eval(p2, list(x = 1.001))
## [1] 8.881784e-15
We can of course also evaluate the polynomials in the different forms
in yacas
using WithValue()
:
# First try with 1.001:
<- paste0("WithValue(x, 1.001, ", pol, ")")
pol_val pol_val
## [1] "WithValue(x, 1.001, (x-1)^7)"
yac_str(pol_val)
## [1] "0"
#... to get result symbolically, use instead the number as a fraction
<- paste0("WithValue(x, 1001/1000, ", pol, ")")
pol_val pol_val
## [1] "WithValue(x, 1001/1000, (x-1)^7)"
yac_str(pol_val)
## [1] "1/1000000000000000000000"
%>% y_fn("Denom") %>% yac_str() pol_val
## [1] "1000000000000000000000"
%>% y_fn("Denom") %>% y_fn("IntLog", "10") %>% yac_str() pol_val
## [1] "21"
<- paste0("WithValue(x, 1001/1000, Expand(", pol, "))")
pol_val pol_val
## [1] "WithValue(x, 1001/1000, Expand((x-1)^7))"
yac_str(pol_val)
## [1] "1/1000000000000000000000"
The reason for the difference is the near cancellation problem when using floating-point representation: Subtraction of nearly identical quantities. This phenomenon is illustrated in the plot below.
This example could of course have been illustrated directly in
R
but it is convenient that Ryacas
provides
expansion of polynomials directly so that students can experiment with
other polynomials themselves.
In R
’s solve()
help file the Hilbert matrix
is introduced:
<- function(n) {
hilbert <- 1:n
i <- 1 / outer(i - 1, i, "+")
H return(H)
}
We will now demonstrate the known fact that it is ill-conditioned (meaning e.g. that it is difficult to determine its inverse numerically).
First we make the Hilbert matrix symbolically:
<- function(n) {
hilbert_sym <- matrix("", nrow = n, ncol = n)
mat
for (i in 1:n) {
for (j in 1:n) {
<- paste0("1 / (", (i-1), " + ", j, ")")
mat[i, j]
}
}
return(mat)
}
<- hilbert(8)
A A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1.0000000 0.5000000 0.3333333 0.25000000 0.20000000 0.16666667 0.14285714
## [2,] 0.5000000 0.3333333 0.2500000 0.20000000 0.16666667 0.14285714 0.12500000
## [3,] 0.3333333 0.2500000 0.2000000 0.16666667 0.14285714 0.12500000 0.11111111
## [4,] 0.2500000 0.2000000 0.1666667 0.14285714 0.12500000 0.11111111 0.10000000
## [5,] 0.2000000 0.1666667 0.1428571 0.12500000 0.11111111 0.10000000 0.09090909
## [6,] 0.1666667 0.1428571 0.1250000 0.11111111 0.10000000 0.09090909 0.08333333
## [7,] 0.1428571 0.1250000 0.1111111 0.10000000 0.09090909 0.08333333 0.07692308
## [8,] 0.1250000 0.1111111 0.1000000 0.09090909 0.08333333 0.07692308 0.07142857
## [,8]
## [1,] 0.12500000
## [2,] 0.11111111
## [3,] 0.10000000
## [4,] 0.09090909
## [5,] 0.08333333
## [6,] 0.07692308
## [7,] 0.07142857
## [8,] 0.06666667
<- hilbert_sym(8)
B1 B1
## [,1] [,2] [,3] [,4] [,5]
## [1,] "1 / (0 + 1)" "1 / (0 + 2)" "1 / (0 + 3)" "1 / (0 + 4)" "1 / (0 + 5)"
## [2,] "1 / (1 + 1)" "1 / (1 + 2)" "1 / (1 + 3)" "1 / (1 + 4)" "1 / (1 + 5)"
## [3,] "1 / (2 + 1)" "1 / (2 + 2)" "1 / (2 + 3)" "1 / (2 + 4)" "1 / (2 + 5)"
## [4,] "1 / (3 + 1)" "1 / (3 + 2)" "1 / (3 + 3)" "1 / (3 + 4)" "1 / (3 + 5)"
## [5,] "1 / (4 + 1)" "1 / (4 + 2)" "1 / (4 + 3)" "1 / (4 + 4)" "1 / (4 + 5)"
## [6,] "1 / (5 + 1)" "1 / (5 + 2)" "1 / (5 + 3)" "1 / (5 + 4)" "1 / (5 + 5)"
## [7,] "1 / (6 + 1)" "1 / (6 + 2)" "1 / (6 + 3)" "1 / (6 + 4)" "1 / (6 + 5)"
## [8,] "1 / (7 + 1)" "1 / (7 + 2)" "1 / (7 + 3)" "1 / (7 + 4)" "1 / (7 + 5)"
## [,6] [,7] [,8]
## [1,] "1 / (0 + 6)" "1 / (0 + 7)" "1 / (0 + 8)"
## [2,] "1 / (1 + 6)" "1 / (1 + 7)" "1 / (1 + 8)"
## [3,] "1 / (2 + 6)" "1 / (2 + 7)" "1 / (2 + 8)"
## [4,] "1 / (3 + 6)" "1 / (3 + 7)" "1 / (3 + 8)"
## [5,] "1 / (4 + 6)" "1 / (4 + 7)" "1 / (4 + 8)"
## [6,] "1 / (5 + 6)" "1 / (5 + 7)" "1 / (5 + 8)"
## [7,] "1 / (6 + 6)" "1 / (6 + 7)" "1 / (6 + 8)"
## [8,] "1 / (7 + 6)" "1 / (7 + 7)" "1 / (7 + 8)"
<- ysym(B1) # convert to yacas symbol
B B
## {{ 1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8},
## { 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9},
## { 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10},
## { 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11},
## { 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12},
## { 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13},
## { 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14},
## { 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15}}
<- solve(A)
Ainv Ainv
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 64 -2016 20160 -92400 221760 -288288
## [2,] -2016 84672 -952560 4656960 -11642400 15567552
## [3,] 20160 -952560 11430720 -58212000 149688000 -204324119
## [4,] -92400 4656960 -58212000 304919999 -800414996 1109908794
## [5,] 221760 -11642400 149688000 -800414996 2134439987 -2996753738
## [6,] -288288 15567552 -204324119 1109908793 -2996753738 4249941661
## [7,] 192192 -10594584 141261119 -776936154 2118916782 -3030050996
## [8,] -51480 2882880 -38918880 216215998 -594593995 856215351
## [,7] [,8]
## [1,] 192192 -51480
## [2,] -10594584 2882880
## [3,] 141261119 -38918880
## [4,] -776936155 216215998
## [5,] 2118916783 -594593995
## [6,] -3030050996 856215352
## [7,] 2175421226 -618377753
## [8,] -618377753 176679358
<- solve(B) # result is still yacas symbol
Binv1 Binv1
## {{ 64, -2016, 20160, -92400, 221760, -288288, 192192, -51480},
## { -2016, 84672, -952560, 4656960, -11642400, 15567552, -10594584, 2882880},
## { 20160, -952560, 11430720, -58212000, 149688000, -204324120, 141261120, -38918880},
## { -92400, 4656960, -58212000, 304920000, -800415000, 1109908800, -776936160, 216216000},
## { 221760, -11642400, 149688000, -800415000, 2134440000, -2996753760, 2118916800, -594594000},
## { -288288, 15567552, -204324120, 1109908800, -2996753760, 4249941696, -3030051024, 856215360},
## { 192192, -10594584, 141261120, -776936160, 2118916800, -3030051024, 2175421248, -618377760},
## { -51480, 2882880, -38918880, 216216000, -594594000, 856215360, -618377760, 176679360}}
<- as_r(Binv1) # convert to R numeric matrix
Binv Binv
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 64 -2016 20160 -92400 221760 -288288
## [2,] -2016 84672 -952560 4656960 -11642400 15567552
## [3,] 20160 -952560 11430720 -58212000 149688000 -204324120
## [4,] -92400 4656960 -58212000 304920000 -800415000 1109908800
## [5,] 221760 -11642400 149688000 -800415000 2134440000 -2996753760
## [6,] -288288 15567552 -204324120 1109908800 -2996753760 4249941696
## [7,] 192192 -10594584 141261120 -776936160 2118916800 -3030051024
## [8,] -51480 2882880 -38918880 216216000 -594594000 856215360
## [,7] [,8]
## [1,] 192192 -51480
## [2,] -10594584 2882880
## [3,] 141261120 -38918880
## [4,] -776936160 216216000
## [5,] 2118916800 -594594000
## [6,] -3030051024 856215360
## [7,] 2175421248 -618377760
## [8,] -618377760 176679360
- Binv Ainv
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.022656e-07 -7.897346e-06 0.0000742203 -0.0002794876 4.892563e-04
## [2,] -7.516163e-06 2.443442e-04 -0.0016348482 0.0022857217 8.007960e-03
## [3,] 6.596705e-05 -1.459678e-03 -0.0013010856 0.0928798467 -4.353278e-01
## [4,] -2.222202e-04 6.508809e-04 0.1018696651 -1.0164485574 3.740381e+00
## [5,] 3.084455e-04 1.383389e-02 -0.4781522453 3.8412284851 -1.304656e+01
## [6,] -1.008651e-04 -3.762038e-02 0.9126323462 -6.7141628265 2.190641e+01
## [7,] -1.274265e-04 3.748520e-02 -0.7896100283 5.5373635292 -1.763533e+01
## [8,] 8.357911e-05 -1.313744e-02 0.2562877908 -1.7438534796 5.464779e+00
## [,6] [,7] [,8]
## [1,] -3.859913e-04 9.308662e-05 1.726509e-05
## [2,] -2.777486e-02 2.954033e-02 -1.067452e-02
## [3,] 8.310004e-01 -7.192242e-01 2.335241e-01
## [4,] -6.457403e+00 5.287978e+00 -1.657728e+00
## [5,] 2.159914e+01 -1.723701e+01 5.309898e+00
## [6,] -3.544080e+01 2.786332e+01 -8.493719e+00
## [7,] 2.811761e+01 -2.189091e+01 6.626538e+00
## [8,] -8.625638e+00 6.669370e+00 -2.008779e+00
max(abs(Ainv - Binv))
## [1] 35.4408
max(abs((Ainv - Binv) / Binv))
## [1] 1.136963e-08
As seen, already for \(n=8\) is the numeric solution inaccurate.