常微分方程式を勉強する

講座 数学の考え方〈7〉常微分方程式論

• 勉強会の目次
• 常微分方程式。正規系で表せる常微分方程式最高次の導関数はそれより次数の低い導関数と元の関数と変数そのもので決まる。今が変化の仕方を決めるというルール。
• 参考→こちら
• 常微分方程式などの既知解の検索サイト→EqWorldについて
• Rの常微分方程式パッケージ deSolove(こちら)
• 常微分方程式の現実世界のために「発散しない条件」を検討しよう→こちらとこちらとこちらと
• 常微分方程式　ordinary differential equation
• 一般解と特殊解と特異解
• 線型のそれ、非線形のそれ
• 初期値問題
• 境界値問題(２階以上だと、両端２か所に初期値を与えられる)
  • 境界条件
    • ディリクレ型
    • ノイマン型
    • 周期型
• 『初期値問題と境界値問題の違いってなんですか？』：→こちら
• 正規形の微分方程式
• リプシッツ連続(こちら)
  • 連続より強い平滑条件
  • コーシー・リプシッツの定理(初期値問題の解が一意に決まるときにリプシッツ条件が出てくる)→こちら
• EqWorldを使って、常微分方程式の厳密解を調べ、Rでプロットしてみよう
  • はaの正負0の３通りに書き分けられる
    • 
    • 
    • 
      • EqWorldのこの記事ではとなっているが、誤植だろう。aを変えてプロットしたときの移行のよさを考慮すれば

• 導関数の動きを見よう

C1<-runif(1)
C2<-runif(1)
as<-seq(from=-1,to=1,length=11)

xs<-seq(from=-10,to=10,length=101)

# ys1,ys2は第１、第２導関数値
ys<-ys1<-ys2<-matrix(0,length(as),length(xs))

for(i in 1:length(ys[,1])){
	a<-as[i]
	absa<-abs(a)
	if(a<0){
		ys[i,]<-C1*sinh(xs*sqrt(absa))+C2*cosh(xs*sqrt(absa))
		# 微分しよう
		dx<-deriv(~C1*sinh(xs*sqrt(absa))+C2*cosh(xs*sqrt(absa)),"xs")
		# 値の評価
		ys1[i,]<-eval(dx)
	}else if(a>0){
		ys[i,]<-C1*sin(xs*sqrt(absa))+C2*cos(xs*sqrt(absa))
		dx<-deriv(~C1*sin(xs*sqrt(absa))+C2*cos(xs*sqrt(absa)),"xs")
		ys1[i,]<-eval(dx)
	}else{
		#ys[i,]<-C1+C2*xs
		ys[i,]<-C1*xs+C2
		dx<-deriv(~C1*xs+C2,"xs")
		ys1[i,]<-eval(dx)
	}
	ys2[i,]<-(-a)*ys[i,]
}

matplot(t(ys),type="l",ylim=range(ys[11,])*5)

library(rgl)
plot3d(xs,ys[1,],ys1[1,])

# a>0 の範囲では
# ysとその第１導関数は周期が半周ずれていて、第２導関数の周期はもとに戻っている

xys<-cbind(xs,ys[8,],ys1[8,],ys2[8,])
plot(as.data.frame(xys))
matplot(xys[,2:4],type="l")

# a<0 の範囲では
# 発散する
xys<-cbind(xs,ys[3,],ys1[3,],ys2[3,])
plot(as.data.frame(xys))
matplot(xys[,2:4],type="l")