เป็นตัวอย่างที่ผมเขียนขึ้นตอนใช้สอนการเขียน function ในภาษา R ครับ
library(magick)
# the formula for the changing of angle over time
# https://en.wikipedia.org/wiki/Pendulum
theta <- function(t,period,theta0){
theta0*cos(2*pi*t/period)
}
# for drawing the pendulum at a given point (x,y), i.e. drawing a line between (0,0) and (x,y)
draw.pendulum <- function(pen.x, pen.y){
pivot.x <- 0
pivot.y <- 0
plot(x=pivot.x,y=pivot.y,type = 'l',ylim = c(-1.5,0.1),asp = 1)
lines(x=c(pivot.x,pen.x),y=c(pivot.y,pen.y))
points(x = pen.x, y=pen.y,col='red')
}
# calculate the position of the pendulum over time
# l-length of the pendulum
# theta0 - the starting angle (in radian)
# runtime - the simulation time (in seconds)
# step - time step
pendulum <- function(theta0,l,runtime,step = 0.1){
g <- 9.8
period <- 2*pi*sqrt(l/g)
times <- seq(0.0,runtime, by = step)
pen.x <- l*sin(theta(times, period = period,theta0 = theta0))
pen.y <- -l*cos(theta(times,period = period,theta0 = theta0))
output <- matrix(c(times,pen.x,pen.y),ncol = 3)
return(output)
}
pendulum.sim <- function(theta0, l, runtime, step = 0.1){
sim.data <- pendulum(theta0 = theta0, l=l, runtime = runtime)
## export each frame as png
for(i in 1:nrow(sim.data)){
png(filename = paste0(tempdir(),"/",i,".png"))
draw.pendulum(sim.data[i,2],sim.data[i,3])
dev.off()
}
frames <- image_read(paste0(tempdir(),"/",1:nrow(sim.data),".png"))
animate <- image_animate(image = frames, fps = 1/step)
return(animate)
}
animate <- pendulum.sim(theta0 = pi/4, l = 1, runtime = 10)
print(animate)
มีคนถามมาว่าช่วยทำตัวอย่างให้ดูหน่อยในการใช้ OpenCL ใน Mathematica ซึ่งผมก็บอกไปแล้วว่ามันตัวอย่างเยอะเลยในhelp Mathematica เอง 🙂 แต่ก็ยังอยากให้ผมทำให้ดู ตัวอย่างที่ทำให้ดูก็เป็นการประมาณค่า Pi โดยอาศัยที่เรียกว่า Monte carlo โดยประมาณจากการสุ่มจุดลงในสี่เหลี่ยมจัตุรัสที่มีวงกลมรัศมีครึ่งหนึ่งของความยาวด้านสี่เหลี่ยมจัตุรัส แล้วดูว่ามีกี่จุดที่อยู่ในวงกลม ซึ่งค่า Pi จะสามรถประมาณได้จาก
Pi = 4 * จำนวนจุดในวงกลม / จำนวนจุดที่สุ่มมาทั้งหมด
