iCAL Modern Calculus
with R and Python


Version 1.0 developed from May 2021 for Cal III by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
https://shen.sdsu.edu
Email:


Chapter 9: Integral in n-dimensional Space

setwd('/Users/sshen/CalculusR')

R plot Fig. 9.1: Example 9.1

library(plotly)
par(mar = c(0,0,0,0.0))
x = y = seq(0, 1, len = 101)
z = outer(x, y, function(x, y){ x^2 + y^2 })
p <- plot_ly(x = ~x, y = ~y, z = ~z,
             type = 'surface')
p
hide_colorbar(p)

Monte Carlo Method for a Double Integral on a Rectangle

N = 10000
x = runif(N, 0, 1)
#z = rnorm(N, 0, 1)
head(x) #Using head() because full results are too long
## [1] 0.19844039 0.62727818 0.82616087 0.08911126 0.26067699 0.76178226
f = function(x){sin(x^2)}
head(f(x))
## [1] 0.039368413 0.383402866 0.630767418 0.007940734 0.067900211 0.548285065
fbar = mean(f(x))
b=1
a =0
intf = fbar*(b-a)
intf
## [1] 0.3141545
N = 100000
x = runif(N, 0, 1)
#z = rnorm(N, 0, 1)
head(x)
## [1] 0.7881409 0.8223316 0.6847545 0.7735544 0.9456567 0.2481613
f = function(x){cos(x^2) + exp(-x^3)}
head(f(x))
## [1] 1.426092 1.353388 1.617441 1.455711 1.055361 1.982938
fbar = mean(f(x))
b=1
a =0
intf = fbar*(b-a)
intf
## [1] 1.712543

Example 9.2: Monte Carlo method for a double integral on a rectangle

N = 10000
x = runif(N, 0, 1)
head(x)
## [1] 0.5349670 0.3934666 0.8857839 0.7491726 0.4789188 0.3395228
par(mar=c(4,4,2,1))
plot(x, pch = 16, lwd = 0.01, cex = 0.09)

hist(x)

x = y = runif(N, 0, 1)
#z = rnorm(N, 0, 1)
head(x)
## [1] 0.0375758521 0.0434734607 0.6545006768 0.0496587814 0.0005776512
## [6] 0.4748596514
f = function(x,y){x^2 + y^2}
head(f(x,y))
## [1] 2.823889e-03 3.779884e-03 8.567423e-01 4.931989e-03 6.673618e-07
## [6] 4.509834e-01
fbar = mean(f(x,y))
fbar
## [1] 0.6680232
V = 1
intf = fbar*V
intf
## [1] 0.6680232

Example 9.3: Monte Carlo method for the volume of a 3D ball

N = 100000
R = 2 #Assume radius R = 2
#This can be adjusted according to needs
p = matrix(runif(3*N, -R, R), ncol =3)
head(p)
##            [,1]       [,2]       [,3]
## [1,] -1.8877396  1.6803422  0.5728856
## [2,] -0.2454757  0.2177943 -1.8880653
## [3,] -1.7991989  0.3564320  0.4331365
## [4,]  1.2550775  1.3862862 -0.8058860
## [5,] -1.2671633 -0.7430721 -1.8507988
## [6,]  0.4098398  0.3982701  1.2065036
k = 0
for(i in 1:N){if((p[i,]%*%p[i,])[1,1] < R^2)
  k = k+1}
k
## [1] 52512
k/N
## [1] 0.52512
V = (k/N)*(2*R)^3
V
## [1] 33.60768
#Formula value
(4/3)*pi*R^3
## [1] 33.51032

Example 9.3-4D case: Monte Carlo method for the volume of a 4D ball

N = 100000
R = 2 #Assume radius R = 2
#This can be adjusted according to needs
p = matrix(runif(9*N, -R, R), ncol =4)
head(p)
##            [,1]       [,2]       [,3]       [,4]
## [1,] -0.0642889  1.0552567  1.3807659  0.3222448
## [2,] -1.2328616  1.7363995  1.3913113  1.9030023
## [3,] -1.6268638  1.7246121 -0.9783401 -1.8250596
## [4,] -1.1357267  1.8153529  1.6019101 -0.2581115
## [5,] -1.2396843 -0.3784469  0.3141563  0.5300585
## [6,] -1.7615425  0.5211506  1.7935314 -1.4376762
k = 0
for(i in 1:N){if((p[i,]%*%p[i,])[1,1] < R^2)
  k = k+1}
k
## [1] 30932
k/N
## [1] 0.30932
V = (k/N)*(2*R)^4
V
## [1] 79.18592

Example 9.3-2D case: Monte Carlo method for the volume of a 2D round disc

N = 100000
R = 2 #Assume radius R = 2
#This can be adjusted according to needs
p = matrix(runif(9*N, -R, R), ncol =2)
head(p)
##            [,1]       [,2]
## [1,] -1.2515638  1.8113389
## [2,]  0.1051441 -1.8148275
## [3,] -0.8670114  0.0860607
## [4,] -0.4505007 -0.7803015
## [5,]  0.1314025  0.6172270
## [6,]  1.2991535  0.7852268
k = 0
for(i in 1:N){if((p[i,]%*%p[i,])[1,1] < R^2)
  k = k+1}
k
## [1] 78646
k/N
## [1] 0.78646
V = (k/N)*(2*R)^2
V  
## [1] 12.58336
#Formula value
pi*R^2
## [1] 12.56637

Example 9.4: Volume of a spherical cap Monte Carlo

N = 100000
R = 1 #Work on the unit ball
phi = 66.5*pi/180
p = matrix(runif(3*N, -R, R), ncol =3)
head(p)
##              [,1]        [,2]       [,3]
## [1,]  0.433250422 -0.19414015 -0.8805101
## [2,]  0.543172437  0.08476663 -0.3777364
## [3,]  0.006962278 -0.65376788  0.4040676
## [4,] -0.271858338  0.12372996  0.8882701
## [5,] -0.158401400  0.43626352 -0.2462571
## [6,] -0.772845326 -0.96960723  0.5902053
p[1:5,]
##              [,1]        [,2]       [,3]
## [1,]  0.433250422 -0.19414015 -0.8805101
## [2,]  0.543172437  0.08476663 -0.3777364
## [3,]  0.006962278 -0.65376788  0.4040676
## [4,] -0.271858338  0.12372996  0.8882701
## [5,] -0.158401400  0.43626352 -0.2462571
r = s = matrix(rep(5,3*N), ncol =3)
head(r)
##      [,1] [,2] [,3]
## [1,]    5    5    5
## [2,]    5    5    5
## [3,]    5    5    5
## [4,]    5    5    5
## [5,]    5    5    5
## [6,]    5    5    5
rep(3,5)
## [1] 3 3 3 3 3
k = 0
for(i in 1:N){if(t(p[i,])%*%p[i,] < R^2 ){
  k = k+1
  r[k,] = p[i,]}
}

k #points inside the ball
## [1] 52429
head(r[1:k,])
##              [,1]        [,2]       [,3]
## [1,]  0.543172437  0.08476663 -0.3777364
## [2,]  0.006962278 -0.65376788  0.4040676
## [3,] -0.271858338  0.12372996  0.8882701
## [4,] -0.158401400  0.43626352 -0.2462571
## [5,] -0.519812271  0.13922050  0.6319331
## [6,] -0.312380176 -0.33131871  0.2026529
k/N
## [1] 0.52429
m=0
for(i in 1:k){if(r[i,3] >= R*sin(phi) ){
  m = m + 1
  s[m,] = r[i,]
  print(paste('m=', m, 's[m,3]=', s[m,3]))}
}
## [1] "m= 1 s[m,3]= 0.921816658694297"
## [1] "m= 2 s[m,3]= 0.940909567289054"
## [1] "m= 3 s[m,3]= 0.93080237088725"
## [1] "m= 4 s[m,3]= 0.975483434274793"
## [1] "m= 5 s[m,3]= 0.981032277457416"
## [1] "m= 6 s[m,3]= 0.961193792056292"
## [1] "m= 7 s[m,3]= 0.922882820479572"
## [1] "m= 8 s[m,3]= 0.921208629850298"
## [1] "m= 9 s[m,3]= 0.92384293442592"
## [1] "m= 10 s[m,3]= 0.965838648378849"
## [1] "m= 11 s[m,3]= 0.975044265855104"
## [1] "m= 12 s[m,3]= 0.921715532895178"
## [1] "m= 13 s[m,3]= 0.94948077108711"
## [1] "m= 14 s[m,3]= 0.979863563552499"
## [1] "m= 15 s[m,3]= 0.987580810207874"
## [1] "m= 16 s[m,3]= 0.966287970077246"
## [1] "m= 17 s[m,3]= 0.97387523856014"
## [1] "m= 18 s[m,3]= 0.950209831353277"
## [1] "m= 19 s[m,3]= 0.952395372558385"
## [1] "m= 20 s[m,3]= 0.934897273778915"
## [1] "m= 21 s[m,3]= 0.917918430641294"
## [1] "m= 22 s[m,3]= 0.956200008280575"
## [1] "m= 23 s[m,3]= 0.926511798519641"
## [1] "m= 24 s[m,3]= 0.928110209293664"
## [1] "m= 25 s[m,3]= 0.932830648962408"
## [1] "m= 26 s[m,3]= 0.92246739147231"
## [1] "m= 27 s[m,3]= 0.924133832100779"
## [1] "m= 28 s[m,3]= 0.925658266060054"
## [1] "m= 29 s[m,3]= 0.938493492547423"
## [1] "m= 30 s[m,3]= 0.968822161201388"
## [1] "m= 31 s[m,3]= 0.931828898377717"
## [1] "m= 32 s[m,3]= 0.92828167276457"
## [1] "m= 33 s[m,3]= 0.947845994960517"
## [1] "m= 34 s[m,3]= 0.947608351241797"
## [1] "m= 35 s[m,3]= 0.949914176017046"
## [1] "m= 36 s[m,3]= 0.965681762900203"
## [1] "m= 37 s[m,3]= 0.932134729810059"
## [1] "m= 38 s[m,3]= 0.943194188643247"
## [1] "m= 39 s[m,3]= 0.928180811461061"
## [1] "m= 40 s[m,3]= 0.927245305385441"
## [1] "m= 41 s[m,3]= 0.963859253562987"
## [1] "m= 42 s[m,3]= 0.92101534595713"
## [1] "m= 43 s[m,3]= 0.925220956560224"
## [1] "m= 44 s[m,3]= 0.942395468242466"
## [1] "m= 45 s[m,3]= 0.920638816896826"
## [1] "m= 46 s[m,3]= 0.947308211121708"
## [1] "m= 47 s[m,3]= 0.927311429288238"
## [1] "m= 48 s[m,3]= 0.99594323663041"
## [1] "m= 49 s[m,3]= 0.91770010208711"
## [1] "m= 50 s[m,3]= 0.973927923012525"
## [1] "m= 51 s[m,3]= 0.982136631850153"
## [1] "m= 52 s[m,3]= 0.961269753519446"
## [1] "m= 53 s[m,3]= 0.917861492838711"
## [1] "m= 54 s[m,3]= 0.963325895369053"
## [1] "m= 55 s[m,3]= 0.94056200934574"
## [1] "m= 56 s[m,3]= 0.925087577663362"
## [1] "m= 57 s[m,3]= 0.941063900012523"
## [1] "m= 58 s[m,3]= 0.942983948625624"
## [1] "m= 59 s[m,3]= 0.93048397032544"
## [1] "m= 60 s[m,3]= 0.921281399670988"
## [1] "m= 61 s[m,3]= 0.918697745539248"
## [1] "m= 62 s[m,3]= 0.963446413632482"
## [1] "m= 63 s[m,3]= 0.951483107637614"
## [1] "m= 64 s[m,3]= 0.955124681815505"
## [1] "m= 65 s[m,3]= 0.922175273299217"
## [1] "m= 66 s[m,3]= 0.9478819668293"
## [1] "m= 67 s[m,3]= 0.959511217661202"
## [1] "m= 68 s[m,3]= 0.925715006422251"
## [1] "m= 69 s[m,3]= 0.929974530357867"
## [1] "m= 70 s[m,3]= 0.935583623126149"
## [1] "m= 71 s[m,3]= 0.934425842016935"
## [1] "m= 72 s[m,3]= 0.974141843151301"
## [1] "m= 73 s[m,3]= 0.919062351342291"
## [1] "m= 74 s[m,3]= 0.947194661013782"
## [1] "m= 75 s[m,3]= 0.965570126194507"
## [1] "m= 76 s[m,3]= 0.931413921527565"
## [1] "m= 77 s[m,3]= 0.964717651251704"
## [1] "m= 78 s[m,3]= 0.952360346913338"
## [1] "m= 79 s[m,3]= 0.957933951169252"
## [1] "m= 80 s[m,3]= 0.975665588397533"
## [1] "m= 81 s[m,3]= 0.917518503032625"
## [1] "m= 82 s[m,3]= 0.976079885847867"
## [1] "m= 83 s[m,3]= 0.949505184777081"
## [1] "m= 84 s[m,3]= 0.9202173454687"
## [1] "m= 85 s[m,3]= 0.930778985377401"
## [1] "m= 86 s[m,3]= 0.93433305574581"
## [1] "m= 87 s[m,3]= 0.979264239314944"
## [1] "m= 88 s[m,3]= 0.941309930756688"
## [1] "m= 89 s[m,3]= 0.940386050380766"
## [1] "m= 90 s[m,3]= 0.940049518831074"
## [1] "m= 91 s[m,3]= 0.940309017896652"
## [1] "m= 92 s[m,3]= 0.945718036033213"
## [1] "m= 93 s[m,3]= 0.921568205114454"
## [1] "m= 94 s[m,3]= 0.925554279703647"
## [1] "m= 95 s[m,3]= 0.921827775426209"
## [1] "m= 96 s[m,3]= 0.970778708811849"
## [1] "m= 97 s[m,3]= 0.952456841710955"
## [1] "m= 98 s[m,3]= 0.955255742650479"
## [1] "m= 99 s[m,3]= 0.937821406405419"
## [1] "m= 100 s[m,3]= 0.966749305836856"
## [1] "m= 101 s[m,3]= 0.934906275942922"
## [1] "m= 102 s[m,3]= 0.95888876914978"
## [1] "m= 103 s[m,3]= 0.983210764359683"
## [1] "m= 104 s[m,3]= 0.953572403173894"
## [1] "m= 105 s[m,3]= 0.931217984762043"
## [1] "m= 106 s[m,3]= 0.994126837234944"
## [1] "m= 107 s[m,3]= 0.961666144896299"
## [1] "m= 108 s[m,3]= 0.98052762215957"
## [1] "m= 109 s[m,3]= 0.92921404261142"
## [1] "m= 110 s[m,3]= 0.93053272459656"
## [1] "m= 111 s[m,3]= 0.963832209818065"
## [1] "m= 112 s[m,3]= 0.957733529154211"
## [1] "m= 113 s[m,3]= 0.941085759084672"
## [1] "m= 114 s[m,3]= 0.93679807940498"
## [1] "m= 115 s[m,3]= 0.936752073932439"
## [1] "m= 116 s[m,3]= 0.91915347892791"
## [1] "m= 117 s[m,3]= 0.928131247404963"
## [1] "m= 118 s[m,3]= 0.935898820869625"
## [1] "m= 119 s[m,3]= 0.928044005297124"
## [1] "m= 120 s[m,3]= 0.957772428169847"
## [1] "m= 121 s[m,3]= 0.926257667131722"
## [1] "m= 122 s[m,3]= 0.935048893559724"
## [1] "m= 123 s[m,3]= 0.967209161259234"
## [1] "m= 124 s[m,3]= 0.938055330887437"
## [1] "m= 125 s[m,3]= 0.986553183291107"
## [1] "m= 126 s[m,3]= 0.923009362537414"
## [1] "m= 127 s[m,3]= 0.952888578176498"
## [1] "m= 128 s[m,3]= 0.923166673164815"
## [1] "m= 129 s[m,3]= 0.981190074700862"
## [1] "m= 130 s[m,3]= 0.936275183688849"
## [1] "m= 131 s[m,3]= 0.936032983008772"
## [1] "m= 132 s[m,3]= 0.919604346621782"
## [1] "m= 133 s[m,3]= 0.973697488661855"
## [1] "m= 134 s[m,3]= 0.943585784174502"
## [1] "m= 135 s[m,3]= 0.973814107012004"
## [1] "m= 136 s[m,3]= 0.97704277606681"
## [1] "m= 137 s[m,3]= 0.917201887350529"
## [1] "m= 138 s[m,3]= 0.928665128536522"
## [1] "m= 139 s[m,3]= 0.941827264148742"
## [1] "m= 140 s[m,3]= 0.954250680748373"
## [1] "m= 141 s[m,3]= 0.91842049639672"
## [1] "m= 142 s[m,3]= 0.945673807989806"
## [1] "m= 143 s[m,3]= 0.937210804782808"
## [1] "m= 144 s[m,3]= 0.924652635119855"
## [1] "m= 145 s[m,3]= 0.918374888598919"
## [1] "m= 146 s[m,3]= 0.970123317558318"
## [1] "m= 147 s[m,3]= 0.967223901301622"
## [1] "m= 148 s[m,3]= 0.919438926037401"
## [1] "m= 149 s[m,3]= 0.948293436784297"
## [1] "m= 150 s[m,3]= 0.921867011580616"
## [1] "m= 151 s[m,3]= 0.985108263324946"
## [1] "m= 152 s[m,3]= 0.926272001583129"
## [1] "m= 153 s[m,3]= 0.929686455987394"
## [1] "m= 154 s[m,3]= 0.951155030168593"
## [1] "m= 155 s[m,3]= 0.948720713611692"
## [1] "m= 156 s[m,3]= 0.95897483592853"
## [1] "m= 157 s[m,3]= 0.927373229991645"
## [1] "m= 158 s[m,3]= 0.952535929623991"
## [1] "m= 159 s[m,3]= 0.926533003337681"
## [1] "m= 160 s[m,3]= 0.918366595171392"
## [1] "m= 161 s[m,3]= 0.969678011257201"
## [1] "m= 162 s[m,3]= 0.933429474942386"
## [1] "m= 163 s[m,3]= 0.976080914027989"
## [1] "m= 164 s[m,3]= 0.929612132720649"
## [1] "m= 165 s[m,3]= 0.939752352423966"
## [1] "m= 166 s[m,3]= 0.964468342252076"
## [1] "m= 167 s[m,3]= 0.962882536929101"
## [1] "m= 168 s[m,3]= 0.928887522313744"
## [1] "m= 169 s[m,3]= 0.960963332559913"
## [1] "m= 170 s[m,3]= 0.978731917683035"
## [1] "m= 171 s[m,3]= 0.925939037464559"
## [1] "m= 172 s[m,3]= 0.966890634037554"
## [1] "m= 173 s[m,3]= 0.933122595772147"
## [1] "m= 174 s[m,3]= 0.945178967900574"
## [1] "m= 175 s[m,3]= 0.919037730433047"
## [1] "m= 176 s[m,3]= 0.918461074586958"
## [1] "m= 177 s[m,3]= 0.9330243174918"
## [1] "m= 178 s[m,3]= 0.946689327713102"
## [1] "m= 179 s[m,3]= 0.941670071799308"
## [1] "m= 180 s[m,3]= 0.947703080717474"
## [1] "m= 181 s[m,3]= 0.947378731798381"
## [1] "m= 182 s[m,3]= 0.94491387763992"
## [1] "m= 183 s[m,3]= 0.973638461437076"
## [1] "m= 184 s[m,3]= 0.919693177100271"
## [1] "m= 185 s[m,3]= 0.977659301366657"
## [1] "m= 186 s[m,3]= 0.943097482901067"
## [1] "m= 187 s[m,3]= 0.931962450966239"
## [1] "m= 188 s[m,3]= 0.923390087205917"
## [1] "m= 189 s[m,3]= 0.923031544778496"
## [1] "m= 190 s[m,3]= 0.947412902489305"
## [1] "m= 191 s[m,3]= 0.96322843618691"
## [1] "m= 192 s[m,3]= 0.928405406419188"
## [1] "m= 193 s[m,3]= 0.92774456506595"
## [1] "m= 194 s[m,3]= 0.932994205970317"
## [1] "m= 195 s[m,3]= 0.920540016144514"
## [1] "m= 196 s[m,3]= 0.929662850219756"
## [1] "m= 197 s[m,3]= 0.980793733615428"
## [1] "m= 198 s[m,3]= 0.935675075743347"
## [1] "m= 199 s[m,3]= 0.983597755432129"
## [1] "m= 200 s[m,3]= 0.950636671856046"
## [1] "m= 201 s[m,3]= 0.950666298624128"
## [1] "m= 202 s[m,3]= 0.966353800147772"
## [1] "m= 203 s[m,3]= 0.92826171265915"
## [1] "m= 204 s[m,3]= 0.955256260931492"
## [1] "m= 205 s[m,3]= 0.948634788859636"
## [1] "m= 206 s[m,3]= 0.938711593393236"
## [1] "m= 207 s[m,3]= 0.924468506593257"
## [1] "m= 208 s[m,3]= 0.932103272527456"
## [1] "m= 209 s[m,3]= 0.918220650870353"
## [1] "m= 210 s[m,3]= 0.958966819103807"
## [1] "m= 211 s[m,3]= 0.930381816346198"
## [1] "m= 212 s[m,3]= 0.954855390358716"
## [1] "m= 213 s[m,3]= 0.955350789707154"
## [1] "m= 214 s[m,3]= 0.928189533762634"
## [1] "m= 215 s[m,3]= 0.940488322172314"
## [1] "m= 216 s[m,3]= 0.958526365458965"
## [1] "m= 217 s[m,3]= 0.95574408210814"
## [1] "m= 218 s[m,3]= 0.929274273570627"
## [1] "m= 219 s[m,3]= 0.971799735445529"
## [1] "m= 220 s[m,3]= 0.927297601941973"
## [1] "m= 221 s[m,3]= 0.940124351996928"
## [1] "m= 222 s[m,3]= 0.9554911586456"
## [1] "m= 223 s[m,3]= 0.936964606866241"
## [1] "m= 224 s[m,3]= 0.966340733226389"
## [1] "m= 225 s[m,3]= 0.921287683304399"
## [1] "m= 226 s[m,3]= 0.959020675159991"
## [1] "m= 227 s[m,3]= 0.930571067612618"
## [1] "m= 228 s[m,3]= 0.954417453613132"
## [1] "m= 229 s[m,3]= 0.920202979817986"
## [1] "m= 230 s[m,3]= 0.922760763205588"
## [1] "m= 231 s[m,3]= 0.925123470369726"
## [1] "m= 232 s[m,3]= 0.96867320779711"
## [1] "m= 233 s[m,3]= 0.948681527283043"
## [1] "m= 234 s[m,3]= 0.992493551690131"
## [1] "m= 235 s[m,3]= 0.921145315747708"
## [1] "m= 236 s[m,3]= 0.969481565523893"
## [1] "m= 237 s[m,3]= 0.991569709498435"
## [1] "m= 238 s[m,3]= 0.984011267311871"
## [1] "m= 239 s[m,3]= 0.917112372349948"
## [1] "m= 240 s[m,3]= 0.934163674712181"
## [1] "m= 241 s[m,3]= 0.923978944774717"
## [1] "m= 242 s[m,3]= 0.989529800135642"
## [1] "m= 243 s[m,3]= 0.94763559801504"
## [1] "m= 244 s[m,3]= 0.92717882944271"
## [1] "m= 245 s[m,3]= 0.929191455710679"
## [1] "m= 246 s[m,3]= 0.994188048411161"
## [1] "m= 247 s[m,3]= 0.928681154269725"
## [1] "m= 248 s[m,3]= 0.91859626956284"
## [1] "m= 249 s[m,3]= 0.946858004201204"
## [1] "m= 250 s[m,3]= 0.917962252628058"
## [1] "m= 251 s[m,3]= 0.950334845110774"
## [1] "m= 252 s[m,3]= 0.974516636691988"
## [1] "m= 253 s[m,3]= 0.940762863960117"
## [1] "m= 254 s[m,3]= 0.925440881866962"
## [1] "m= 255 s[m,3]= 0.95187696814537"
## [1] "m= 256 s[m,3]= 0.924293751362711"
## [1] "m= 257 s[m,3]= 0.963314971886575"
## [1] "m= 258 s[m,3]= 0.954340252093971"
## [1] "m= 259 s[m,3]= 0.93427815893665"
## [1] "m= 260 s[m,3]= 0.972898615524173"
s[1:m]
##   [1]  0.2109093033  0.1829032903  0.3279029210  0.0566157419  0.1559664020
##   [6] -0.1420859145 -0.1629892876 -0.3034387538  0.0080781193 -0.1042655818
##  [11]  0.0946708843 -0.0634831167 -0.0862028264  0.0957515826  0.0713915620
##  [16]  0.1176622817 -0.0309391120 -0.1057061804 -0.1071871403 -0.1127826064
##  [21] -0.2559709558  0.0397065524 -0.2541542677 -0.2184718577 -0.3454627586
##  [26] -0.0527270185 -0.3519836436  0.0806141901 -0.0416481108  0.1988918586
##  [31]  0.1582821873 -0.0806001052  0.0875498494  0.1321594580 -0.0557956602
##  [36]  0.1344987154 -0.1616505198  0.2451121276  0.1508668028 -0.0880249552
##  [41] -0.1106677121 -0.3199094399 -0.2297610785  0.1294357781  0.0803842782
##  [46] -0.2195653571  0.2336007985 -0.0198512948 -0.0448082932 -0.0962785450
##  [51] -0.0697708395  0.0775085208  0.0729171457 -0.0500901076 -0.0709717087
##  [56] -0.0925654983 -0.0514518782  0.1725929789 -0.0453793444 -0.1548932847
##  [61] -0.2493298030 -0.0436264137 -0.2926847357 -0.0960618556  0.1583829247
##  [66] -0.0757010141 -0.0404875772 -0.0888525955  0.0514713735  0.2797548240
##  [71]  0.0531560876  0.1555669927  0.3647809913  0.2625277047 -0.0650965204
##  [76]  0.3078141236 -0.0235370602 -0.0769135384 -0.1615397627 -0.1480623465
##  [81] -0.3527155458  0.0284948256 -0.0568872453 -0.0870719948 -0.1164517617
##  [86]  0.2201016941  0.1203966276  0.0416955207  0.0152256778 -0.0276253079
##  [91] -0.2524985601  0.0167203569  0.3239782285 -0.1407553516 -0.0691206669
##  [96] -0.1282598833 -0.0239774301  0.2292598421  0.0676513780 -0.0963094961
## [101]  0.0850641271  0.1259961235  0.0824191477  0.2721560108 -0.1321267383
## [106] -0.0247510257 -0.1884576320  0.0808790941  0.0988698686  0.1477726288
## [111] -0.0606763572  0.0492022336  0.1133265533  0.3313317969 -0.0893175686
## [116] -0.3197054588  0.0215580883 -0.3359424882  0.2425898463 -0.0204441897
## [121]  0.0133355306  0.0192454858 -0.0755982199  0.0865013255  0.0150913591
## [126] -0.0666618594  0.1473219828  0.1741376496 -0.1574787125 -0.0393989356
## [131] -0.0648840209 -0.3472042894  0.0700079259  0.0178911295 -0.0928253052
## [136] -0.0400279239 -0.1266932343 -0.3507259493 -0.1251854491 -0.2345020659
## [141]  0.1348159821  0.1481963098 -0.1898893774 -0.0790258255  0.1961481660
## [146] -0.1306078029  0.0040883063 -0.0132956994 -0.2742838985  0.1503916890
## [151] -0.0922490996 -0.0198594499  0.0708531672  0.0752296075  0.0802421113
## [156]  0.0351957558 -0.0377694988  0.1760769254  0.1728293709  0.1948543307
## [161] -0.1506409803  0.3169869483 -0.1902791434  0.0887382450 -0.0531494692
## [166] -0.0951323635  0.0332103088  0.2393720681 -0.2332251580  0.0870058942
## [171]  0.2998819579  0.1250235187  0.1208798676  0.0026376732  0.2376591722
## [176] -0.0556477671 -0.0760345627 -0.2379241381  0.2569029080  0.1876717913
## [181] -0.2121319612  0.1818076246 -0.0775751970 -0.2167260502  0.0556163103
## [186] -0.1214036224 -0.1983099105  0.3259726861 -0.0447178883  0.0465284940
## [191] -0.0259109298 -0.1201040475 -0.0895912307 -0.2284821719 -0.2648461070
## [196] -0.2616426614  0.1207628050  0.2891282840  0.1276352149 -0.1794495280
## [201]  0.2973853094  0.0283347168  0.2349228095 -0.1044310248  0.0591846672
## [206] -0.1210267027 -0.2881859583 -0.1937646889  0.0940797310  0.1832588720
## [211] -0.0128406524  0.1248352001  0.1736709042  0.1903305016  0.1327950153
## [216] -0.1361503694  0.0532212490 -0.1040221290 -0.0514006144 -0.1412864863
## [221]  0.1161118424  0.0953468936  0.1112248641 -0.1523253527 -0.3351295898
## [226]  0.0822833376  0.0086781499 -0.2113626529  0.0172704570  0.3233398888
## [231]  0.0749199279  0.0716366554 -0.1107610390  0.0949401902  0.1961997775
## [236] -0.1693492946 -0.0335970633 -0.0396351377  0.0004740171 -0.0453896234
## [241]  0.1123070521 -0.0123655871 -0.0474910717  0.1152510527 -0.1805820777
## [246] -0.0251660352 -0.0089528584 -0.1403491204 -0.0790415420 -0.1314180973
## [251]  0.1775764469  0.1244415473 -0.1635731626  0.0367435827  0.1769915354
## [256]  0.3044887427 -0.0128437215  0.0637261108  0.2216758751  0.1366160922
s[1:5,]
##            [,1]        [,2]      [,3]
## [1,] 0.21090930 -0.25722607 0.9218167
## [2,] 0.18290329 -0.10140653 0.9409096
## [3,] 0.32790292  0.05979938 0.9308024
## [4,] 0.05661574 -0.07320641 0.9754834
## [5,] 0.15596640 -0.10847711 0.9810323
R*sin(phi) 
## [1] 0.9170601
plot(s[1:m,3], type = 'o')

m #points inside the cap
## [1] 260
m/N
## [1] 0.0026
V = (m/N)*(2*R)^3
V
## [1] 0.0208
#[1] 0.02024 i.e., the cap volume is 0.02024*R^3
0.02024/(4*pi/3) 
## [1] 0.004831944
#[1] 0.004831944 means about 0.5% of the ball's volume

#R plot the points in y-z plane when x=0
par(mar = c(4.2,4,0.5,0.5))
plot(p[,2], p[,3], pch=16, cex =0.1,
     xlab = 'y', ylab = 'z')
points(r[1:k,2], r[1:k,3], 
       pch=16, cex =0.2, col= 'red')
points(s[1:m,2], s[1:m,3], 
       pch =16, cex =0.4, col = 'blue')

#Plot the points in x-z plane when x=0
par(mar = c(4.2,4,0.5,0.5))
plot(p[,1], p[,3], pch=16, cex =0.1,
     xlab = 'x', ylab = 'z')
points(r[1:k,1], r[1:k,3], 
       pch=16, cex =0.2, col= 'red')
points(s[1:m,1], s[1:m,3], 
       pch =16, cex =0.4, col = 'blue')

Monte Carlo integral without plotting

N = 100000
R = 1 #Work on the unit ball
phi = 66.5*pi/180
p = matrix(runif(3*N, -R, R), ncol =3)
k = 0
for(i in 1:N){if(t(p[i,])%*%p[i,] <= R^2 && 
                 p[i,3] >= R*sin(phi)){k = k+1}}
(k/N)*(2*R)^3 #volume of the cap
## [1] 0.01984

Monte Carlo integral int_D exp(-x2-y2)

N = 100000
R = 2 #Work on the unit disk
p = q = matrix(runif(2*N, -R, R), ncol =2)
f = function(x,y){exp(-x^2 - y^2)}
k = 0
for(i in 1:N){if(t(p[i,])%*%p[i,] <= R^2){
  k = k+1
  q[k,] = p[i,]}
}
k
## [1] 78455
x = q[1:k,1]; y = q[1:k,2];
mean(f(x, y))*pi*R^2 #The integral
## [1] 3.083626

Monte Carlo method on a sphere

N = 10000
R = 1
phi = 75
#Find points randomly uniformly distributed on a sphere
#Algorithm ref: Marsaglia (1972) Ann. Math. Stat. 43,645-646 
p = q = matrix(runif(2*N, -R, R), ncol =2)
k = 0
for(i in 1:N){if(t(p[i,])%*%p[i,] <= R^2){
  k = k+1
  q[k,] = p[i,]}
}
M = k

Find the points on a polar cap

u = q[1:M, 1]
v = q[1:M, 2]
s = u^2 + v^2
x = 2*u*sqrt(1-s)
y = 2*v*sqrt(1-s)
z = 1-2*s
M
## [1] 7797
x2 = y2 = z2 = x
m=0
for(i in 1:M){if(abs(z[i]) >= R*sin(phi*pi/180)){
  m = m+1 
}
}
Mcap = m
Mcap
## [1] 252
#The area of the Artic polar cap
0.5*(Mcap / M)*4*pi*R^2
## [1] 0.2030733

Example 9.6: Parametric integration for a line integral

f = function(t){t + t^2 + 2*t^3}
integrate(f, 0, 2)
## 12.66667 with absolute error < 1.4e-13
#Example 9.6 from the original integral formation
n = 1000
a = 0
b = 2
m = n+1
x = seq(a, b, len = m)
y = x^2
dx = (b-a)/n
dy = 2*x*dx
I = sum((x-y)*dx + (x + y)*dy)
I #The line integral is equal to
## [1] 12.68868
f = function(x){x^2}
integrate(f, 0, 1)
## 0.3333333 with absolute error < 3.7e-15
f = function(t){sqrt(1 + 3*t^2)}
integrate(f, -9, 6)
## 103.483 with absolute error < 0.0054