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")

อีกตัวอย่างครับ สำหรับแก้สมการลอเร้นซ์   แต่แทนที่จะเซฟสมการและพารามิเตอร์ใน 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')

ถ้าสนใจอยากติดตั้ง 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 ได้เลยครับ