Solving Equation on the Computer

• Solving Equation on the Computer と こちら のコメントで教えていただいた数値計算について
• 数値計算
  • Euler's method
  • Improved Euler method
  • Runge-Kutta method
  • その他（Wiki）
• について
• Euler's method
  • 
• Improved Euler method
  • 2点の傾きの平均をとる
  • 
  • 
• Runge-Kutta method
  • 
  • 

• ソースで用いた条件は以下
  • 
  • 初期条件：
  • の範囲の計算

f<-function(x){exp(-x)}
T0<-0
T<-10
x0<-0

# 解析解
F<-function(t){log(t+1)}
t<-seq(from=T0,to=T,by=0.01)
plot(t,F(t),type="l",col=1,ylim=c(0,F(T)))

# オイラー法
dt<-0.5
Nt<-T/dt
xold<-x0
xn<-c(x0)
for(nt in 1:Nt){
	xnew<-xold+f(xold)*dt
	xn<-c(xn,xnew)
	xold<-xnew
	}

t<-seq(from=T0,to=T,by=dt)
par(new=T)
plot(t,xn,col=2,ylim=c(0,F(T)))

# improved-オイラー法
dt<-0.5
Nt<-T/dt
xold<-x0
xn<-c(x0)
for(nt in 1:Nt){
	x1<-xold+f(xold)*dt
	xnew<-xold+(f(xold)+f(x1))/2*dt
	xn<-c(xn,xnew)
	xold<-xnew
	}

t<-seq(from=T0,to=T,by=dt)
par(new=T)
plot(t,xn,col=3,ylim=c(0,F(T)))

# ルンゲ-クッタ法：fourth-order Runge-Kutta
dt<-1
Nt<-T/dt
xold<-x0
xn<-c(x0)
for(nt in 1:Nt){
	k1<-f(xold)*dt
	k2<-f(xold+k1/2)*dt
	k3<-f(xold+k2/2)*dt
	k4<-f(xold+k3)*dt
	xnew<-xold+(k1+k2*2+k3*2+k4)/6
	xn<-c(xn,xnew)
	xold<-xnew
	}

t<-seq(from=T0,to=T,by=dt)
par(new=T)
plot(t,xn,col=4,ylim=c(0,F(T)))

legend(x=c(4,10),y=c(0.5,1),c("Analytical solution","Euler's method","improved Euler method","fourth-order Runge-Kutta method"),col=c(1,2,3,4),lty=c(1,0,0,0),pch=c(-1,1,1,1))