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

 

maemod (แม่มด)

ผมเขียน R package ชื่อ maemod สำหรับช่วยให้คนเพิ่งเริ่มเรียนรู้ใช้งาน R และอยากจะแก้สมการเชิงอนุพันธ์เบื้องต้น (ordinary differential equation: ode) ใน R ให้ทำได้ง่ายขึ้นครับ โดยที่ผู้ใช้งานก็เพียงเขียนสมการ ODE (ในรูปแบบคำสั่งของ R), ค่าของตัวแปรพารามิเตอร์ และค่าเริ่มต้นต่าง ๆ ของสมการ ตามตัวอย่างนี้ครับ
1. สร้างไฟล์สมมุติชื่อ test.txt ที่มีสมการตามนี้ครับ

!Equations
 dY1<-Y2
 dY2<-a*sin(Y1)+sin(t)
 
!Parameters
 a= -1.0
 
!Inits
 Y1=0,Y2=0
 
!Outputs
 c(dY1,dY2)
 
!ExtraFunctions

!Plots
!MAEMOD_End

โดยที่บรรทัดที่อยู่ถัดจาก !Equations คือสมการ ode, !Parameters คือค่าของตัวแปรต่าง ๆ, !Inits ค่าเริ่มต้นของแต่ล่ะสมการ ode, !Outputs คือlistของoutputที่ต้องการโดยที่คอลัมน์แรกซึ่งคือค่าของเวลาจะถูกเพิ่มเข้าไปอัตโนมัติเสมอ, !ExtraFunctions ในกรณีที่มี function หรือตัวแปรพิเศษก็สามารถใส่ในส่วนนี้ได้ครับ และสุดท้าย !MAEMOD_End เป็นตัวบอกว่าสมการหรือค่าต่างๆสิ้นสุดที่บรรทัดนี้ครับ พูดง่ายๆก็คือมีคีย์เวิร์ด อยู่ 6 คำคือ !Equations, !Parameters, !Inits, !Outputs, !ExtraFunctions และ !MAEMOD_End ครับ ไม่จำเป็นต้องเรียงตามลำดับนี้นะครับ อันจะขึ้นก่อนหลังก็ได้
2. วิธีรันก็เพียงใช้คำสั่ง maemod.ode เช่น
ต้องการจะคำนวณสมการข้างบนโดยที่ ค่า t เริ่มต้นที่ 0 เพิ่มขึ้นที่ล่ะ 0.1 จนถึง 20*pi

out<-maemod.ode(input.filename = "test.txt",timegrid = seq(0,20*pi,0.1))
head(out)  time Y1 Y2 
[1,] 0.0 0.0000000000 0.000000000 
[2,] 0.1 0.0001663662 0.004991383 
[3,] 0.2 0.0013277223 0.019866019 
[4,] 0.3 0.0044590899 0.044327146 
[5,] 0.4 0.0104962890 0.077883467 
[6,] 0.5 0.0203162441 0.119856470
plot(out[,c(2,3)],col="red")

heart

อีกตัวอย่างครับ สำหรับแก้สมการลอเร้นซ์   แต่แทนที่จะเซฟสมการและพารามิเตอร์ใน text file เราสามารถสร้างเป็นตัวแปรเก็บข้อความของระบบที่มีตัวแปรต่างๆได้ตามนี้เลยครับ (สังเกตุตัวแปรที่ชื่อ Lorenz ครับ)

library(scatterplot3d)
library(maemod)
lorenz<-"
 !Equations
 dx<-sigma*(y-x)
 dy<-x*(rho-z)-y
 dz<-x*y-beta*z
 
 !Parameters
 sigma=10, beta=2.66666, rho=28
 
 !Inits
 x=-10, y=-12, z=30.05
 
 !Outputs
 c(dx,dy,dz)

 !Plots
 !ExtraFunctions
 
 !MAEMOD_End
 "
lorenzout<-maemod.ode(input.text = lorenz, timegrid = seq(0,30,0.01))

scatterplot3d(lorenzout[,c(2,3,4)],type = 'l')

lorenzถ้าสนใจอยากติดตั้ง package นี้ก็สามารถทำได้ตามนี้เลยครับ

 library(devtools)
 install_github("slphyx/maemod")

code ทั้งหมดอยู่ที่ https://github.com/slphyx/maemod ครับ

จริงๆแล้วใน maemod นี้เราสามารถเลือกใช้ methods ต่างๆ เช่น lsoda, rk4 และอื่นๆอีกมากมาย สำหรับแก้ระบบสมการอนุพันธ์นี้ได้ครับ โดยวิธีการเลือกใช้ methods เหล่านี้สามารถดูได้ที่ package deSolve (http://desolve.r-forge.r-project.org/index.html) เลยครับเพราะค่าพารามิเตอร์สำหรับmethods ต่างๆนี้จะถูกส่งไปเรียกจาก deSolve อีกทีครับ

หากอยากช่วยเขียนปรับปรุงก็ Pull requests หรือ Fork ได้เลยครับ