จากที่เคยเขียนไว้เกี่ยวกับ r package อันหนึ่งที่ผมเขียนชื่อ maemod สำหรับช่วยให้คนที่สนใจอยากคำนวณพวก ode ได้ง่ายขึ้น (ดูข้างล่าง)
มีคนสนใจว่าถ้าพวกตัวแปร 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