読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ggplot2 で接線場や方向場を描く

微分方程式 R グラフ

接線場

微分方程式の解の振る舞いをみるために接線場を描いてみる。

以下で与えられるロジスティック方程式を考える。

x'(t)=x(t)(1-x(t))

接線場は t-x 平面の各点に傾き x(t)(1-x(t)) の小さな線分を描いたものである。

logis <- function(t,x,a=1){
  a*x*(1-x)
}
#点(t1,x1)を通り, 傾きが slope の直線
TheLine <- function(t1,x1,slope,t){
  slope(t1,x1)*(t-t1)+x1
}
tv <- xv <- seq(-1.5,1.5,by=0.1)
slopefield <- function(tv,xv,f){
  tx <-expand.grid(x=tv,y=xv)
  data.frame(t=tx[,1],
             x=tx[,2],
             t2=tx[,1]+0.05,
             x2=TheLine(tx[,1],tx[,2],f,tx[,1]+0.05))
}
out <- slopefield(tv,xv,logis)
library(ggplot2)
library(cowplot)
ggplot(out,aes(x=t,y=x))+
  geom_segment(aes(xend=t2,yend=x2))

f:id:abrahamcow:20160607214449p:plain

このグラフから x(0) > 0 なる解は 1 に近づくこと、x(0) < 0 なる解は t が大きくなるとき無限大に発散することがわかる。

また x=0 と x=1 で傾きが 0、解は定数関数になることがわかる。

とはいえロジスティック方程式は解析的に解けるので接線場のありがたみが少ない。

今度は、

 x'(t)= x(1-x)-(1+\sin(2\pi t))

の接線場を描いてみる。

こんな風だ。

f <- function(t,x,a=1,h=1){
  a*x*(1-x)+h*(1+sin(2*pi*t))
}
out <- slopefield(tv,xv,f)
library(ggplot2)
library(cowplot)
ggplot(out,aes(x=t,y=x))+
  geom_segment(aes(xend=t2,yend=x2))

f:id:abrahamcow:20160607215952p:plain

方向場

微分方程式

 x'=y\\y'=-x

を考える。

接線場と同じように各点 (x,y) において x 成分が yy 成分が -x のベクトルを考えることで、解曲線の動きを把握できる。

ただし、ベクトルをそのまま図示すると互いに重なって見づらくなるので、ベクトルの向きを保ったまま適当に縮小して、方向場を描く。

f <- function(x,y){
  c(y,-x)
}
xv <- yv <- seq(-1,1,by=0.1)
vectorfield <- function(f,xv,yv){
  xy <-expand.grid(x=xv,y=yv)
  xy2 <-mapply(f,xy[,1],xy[,2])
  data.frame(x=xy[,1],y=xy[,2],x2=xy2[1,],y2=xy2[2,])
}
out <-vectorfield(f,xv,yv)
library(ggplot2)
ggplot(out,aes(x=x,y=y))+
  geom_segment(aes(xend=x+x2/10,yend=y+y2/10),arrow = arrow(length = unit(0.3,"cm")))

f:id:abrahamcow:20160607221228p:plain

参考文献

Hirsch・Smale・Devaney 力学系入門―微分方程式からカオスまで

Hirsch・Smale・Devaney 力学系入門―微分方程式からカオスまで

expand.grid関数・outer関数を使って2重forループを美しく書く - My Life as a Mock Quant

そういえば接線場という言葉はあまり見かけない。slope field に他の訳語があるのだろうか。