a=0.5
 b=2
 cs=c("red","skyblue","navy","darkgrey")
 
 f<-function(x) {
   return( 0.9*sin(9*x^0.6+0.3)+x+0.9 )
 }
 
 mc_plot<-function(n,new=1) {
   set.seed(6911)
   x=seq(a*0.6,b*1.05,length=700)
   if (n>0) {
     xi=sort((b-a)*runif(n)+a)
   } else {
     xi=NA
   }
   Xi=sort((b-a)*runif(new)+a)
   pts=sort(c(a,xi,Xi,b))
   plot(x,f(x),type="n",ylim=range(f(x),0),xlab="",ylab="")
   polygon(c(pts,b,a),c(f(pts),0,0),col=cs[4], border=cs[2],lwd=2)
   for (i in xi) {
     segments(i,0,i,f(i),col=cs[2],lwd=3)
   }
   for (i in Xi) {
     segments(i,0,i,f(i),col=cs[3],lwd=3)
   }
   lines(x,f(x),col=cs[1],lwd=3)
   ##calcualte area (integral)
   area=0
   for (p in 1:(length(pts)-1) ) {
     area <- area + (pts[p+1]-pts[p]) * (f(pts[p+1])+f(pts[p])) / 2
   }
   text(x=0.45,y=3.2, labels=substitute(integral(f(x)*dx, A, B)%~~% AREA,list(A=a,B=b,AREA=format(area,digits=3))),cex=1.5,col=cs[4],pos=4)
 }
 
 
 png(filename = "NumInt-MC.png", width=1200, height=900, pointsize = 12)
 par(bg="grey90",mfrow=c(2,2), oma=c(0,0,3,0), cex.axis=0.85)
 mc_plot(n=0)
 mc_plot(n=1)
 mc_plot(n=2)
 mc_plot(n=3,new=5)
 title(main="Numerische Integration mit Monte Carlo (Trapezmethode)", outer=TRUE,cex.main=1.9)
 dev.off()