maemod package

จากที่เคยเขียนไว้เกี่ยวกับ r package อันหนึ่งที่ผมเขียนชื่อ maemod สำหรับช่วยให้คนที่สนใจอยากคำนวณพวก ode ได้ง่ายขึ้น (ดูข้างล่าง)

maemod (แม่มด)

มีคนสนใจว่าถ้าพวกตัวแปร state ต่างๆนั้นเป็นแบบ array จะทำอย่างไร ผมเลยเขียนตัวอย่างพร้อมกับเพิ่มความสามารถด้านarrayนี้เข้าไป พอใช้ได้ไปคร่าวๆก่อน ดูตัวอย่างข้างล่างนี้ครับ

# Example from Berkeley Madonna
 # for using Array
 # METHOD RK4
 #
 # STARTTIME = 0
 # STOPTIME = 20
 # DT = 0.02
 #
 # d/dt (A[1]) = -k[1]*A[1]
 # d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i]
 # d/dt (A[n]) = k[n-1]*A[n-1]
 #
 # init A[1] = 100
 # init A[2..n]=0
 #
 # k[1..n] =1
 # n = 14
 #
 #

arrayEx <-'
 !ExtraFunctions
 n<-14

!Equations
 dA <-rep(0,n)
 dA[1] <- -k[1]*A[1]
 #ArrayB
 "dA[i.indx] <- k[i.indx-1] *
 A[i.indx-1]-k[i.indx]*A[i.indx]"%@@%list("i.indx"%=>%2:(n-1))
 #ArrayE
 dA[n] <- k[n-1]*A[n-1]

!Parameters
 1,1,1,1,1,1,1,1,1,1,1

!Inits
 100,0,0,0,0,0,0,0,0,0,0,0,0,0

!Outputs
 c(dA)

!Plots

!MAEMOD_End
 '
 out<-maemod.ode(input.text = arrayEx ,timegrid = seq(0,20,0.02),sys.template = Maemod_Array)
 plot(out[,c(1,2)],type = 'l')
 cols<-rainbow(14)
 for(i in 3:n) lines(out[,c(1,i)],col=cols[i])

ในเบื้องต้นนี้ ตัวstate กับ parameter จะต้องเป็น A กับ k เท่านั้น จากตัวอย่างเราสามารถที่จะdefine ตัวแปรที่ index ต่างได้ เช่น ต้องการ

d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i]

โดยที่ i มีค่าตั้งแต่ 2 ถึง  n

เราก็สามารถพิมพ์ตามนี้ได้

#ArrayB

"dA[i] <- k[i-1]*A[i-1]-k[i]*A[i] "%@@%list('i'%=>%2:n)

#ArrayE

เจ้าเครื่องหมาย %@@% เป็นตัวบอกว่าเราจะใช้ array กับ สมการด้านซ้ายมือของเคื่องหมายและindexจะอยู่ด้านขวาของเครื่องหมาย ซึ่งจะต้องเป็นตัวแปรแบบ list

ส่วนเครื่องหมาย %=>% บอกว่าสัญลักษณ์หรือรูปแแบบindexด้านขวาจะมีค่าเปลี่ยนไปตาม vector ด้านซ้าย เช่น ‘i’%=>%2:n จะหมายถึง i มีค่าตั้งแต่ 2 ถึง n

ส่วน #ArrayB และ #ArrayE เป็น block ที่บอกว่ามีการใช้ ตัวแปรแบบ array ในข้อความที่อยู่ใน block นี้ครับ

codeที่มีarray blockนี้เวลาเอาไปใช้กับ maemod.ode จะต้องใช้กับtemplate เฉพาะที่ชื่อ Maemod_Array นั่นก็คือต้องเซ็ต option ชื่อ sys.template = Maemod_Array ครับ

สนใจอยากลองใช้งานก็ดูได้ที่

https://github.com/slphyx/maemod

 

สมการดอกกุหลาบ

Rose[x_,theta_]:=Module[{phi=(Pi/2)Exp[-theta/(8 Pi)],X=1-(1/2)((5/4)(1-Mod[3.6 theta,2 Pi]/Pi)^2-1/4)^2},
y=1.95653 x^2 (1.27689 x-1)^2 Sin[phi];
r=X(x Sin[phi]+y Cos[phi]);
{r Sin[theta],r Cos[theta],X(x Cos[phi]-y Sin[phi])}
];

ParametricPlot3D[Rose[x,theta],{x,0,1},{theta,-2 Pi,15 Pi},PlotPoints->{25,576},PlotStyle->Red,Mesh->None,Axes->False,Boxed->False]

rose

ลอกมาจาก www.bugman123.com

HappyValentinesDay

R code สำหรับวันแห่งความรัก 🙂

t<-seq(-pi, pi,0.05)

colors<-rainbow(length(t))

x<-16*sin(t)^3

y<-13*cos(t)-5*cos(2*t)-2*cos(3*t)-1*cos(4*t)

p<-c(“H”, “a”, “p”,”p”,”y”,”V”,”a”,”l”,”e”,”n”,”t”,”i”,”n”,”e”,”s”,”D”,”a”,”y”)

plot(x,y, type=”p”,pch=p,col=colors,xlab=”X”,ylab=”Y”)

heart R