数値シミュレーション入門　サイエンス社　河村哲也先生

```{r}


#p.41

#d^2x/dt^2=-x, x(0)=0, dx(0)/dt=1 をオイラー法で解く

x=0
y=1
t=0

h=0.0001

ite=100000
x_vec=c()
y_vec=c()
t_vec=c()

for(l in 1:ite){

x_vec=c(x_vec,x)
y_vec=c(y_vec,y)
t_vec=c(t_vec,t)
x=x+h*y
y=y-h*x
t=t+h
#print(paste0(x,":",y))
  
}

plot(t_vec,x_vec)
plot(t_vec,y_vec)

```


```{r}

#p.45

#3階の常微分方程式d^3x/dt^3=-t+2*dx/dt+3*d^2x/dt^2

x=1 #x(0)=1
y=1 #y(0)=1
z=1 #z(0)=1
t=0


ite=10000
h=0.0001

x_vec=c()
y_vec=c()
z_vec=c()
t_vec=c()

for(l in 1:ite){

x_vec=c(x_vec,x)
y_vec=c(y_vec,y)
z_vec=c(z_vec,z)
t_vec=c(t_vec,t)

x=x+h*y
y=y+h*z
z=z+h*(2*y+3*z-t)
t=t+h
  
}

res=cbind(t_vec,x_vec,y_vec,z_vec)

plot(t_vec,x_vec,xlim=c(min(t_vec),max(t_vec)),ylim=c(min(c(x_vec,y_vec,z_vec)),max(c(x_vec,y_vec,z_vec))),xlab="t",ylab="values",col=2)

par(new=T)

plot(t_vec,y_vec,xlim=c(min(t_vec),max(t_vec)),ylim=c(min(c(x_vec,y_vec,z_vec)),max(c(x_vec,y_vec,z_vec))),xlab="t",ylab="values",col=3)

par(new=T)

plot(t_vec,z_vec,xlim=c(min(t_vec),max(t_vec)),ylim=c(min(c(x_vec,y_vec,z_vec)),max(c(x_vec,y_vec,z_vec))),xlab="t",ylab="values",col=4)




```

```{r}

#p.48

#4次のルンゲクッタ

#dx/dt=t^2+t+1-(2*t+1)*x+x^2、 x(0)=0.5、　を解く
#解析解x(t)=(t*exp(t)+t+1)/(exp(t)+1)

f=function(tt,xx){
  
return(tt^2+t+1-(2*tt+1)*xx+xx^2)  
  
}

t=0.5
x=0.5

x_vec=c()
t_vec=c()

ite=100
h=0.1

for(l in 1:ite){
  
x_vec=c(x_vec,x)
t_vec=c(t_vec,t)  
  
s1=f(t,x)
s2=f(t+h/2,x+h*s1/2)
s3=f(t+h/2,x+h*s2/2)
s4=f(t+h,x+h*s3)

x=x+h*(s1+2*s2+2*s3+s4)/6
t=t+h

}

plot(t_vec,x_vec,xlab="t",ylab="x",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(x_vec),max(x_vec)),col=2)

par(new=T)

plot(t_vec,(t_vec*exp(t_vec)+t_vec+1)/(exp(t_vec)+1),xlab="t",ylab="x",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(x_vec),max(x_vec)),col=3)

```

```{r}

#p.48 

#4次のルンゲクッタ　連立微分方程式

#dx/dt=-3*x-2*y+2*t

f=function(tt,xx,yy){
  
return(-3*xx-2*yy+2*tt)  
  
}

#dy/dt=2*x+y-sin(t)

g=function(tt,xx,yy){
  
return(2*xx+yy-sin(tt))  
  
}

t=0
x=4.5
y=-6.5

x_vec=c()
y_vec=c()
t_vec=c()

ite=2000
h=0.01

for(l in 1:ite){
  
x_vec=c(x_vec,x)
y_vec=c(y_vec,y)
t_vec=c(t_vec,t)  
  
s11=f(t,x,y)
s12=g(t,x,y)
s21=f(t+h/2,x+h*s11/2,y+h*s11/2)
s22=g(t+h/2,x+h*s12/2,y+h*s12/2)
s31=f(t+h/2,x+h*s21/2,y+h*s21/2)
s32=g(t+h/2,x+h*s22/2,y+h*s22/2)
s41=f(t+h,x+h*s31,y+h*s31)
s42=g(t+h,x+h*s32,y+h*s32)
x=x+h*(s11+2*s21+2*s31+s41)/6
y=y+h*(s12+2*s22+2*s32+s42)/6
t=t+h

}

plot(t_vec,x_vec,xlab="t",ylab="x",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(x_vec),max(x_vec)),col=2)

plot(t_vec,y_vec,xlab="t",ylab="y",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(y_vec),max(y_vec)),col=3)

```

```{r}

#p.48 

#オイラー法　連立微分方程式

#dx/dt=-3*x-2*y+2*t

f=function(tt,xx,yy){
  
return(-3*xx-2*yy+2*tt)  
  
}

#dy/dt=2*x+y-sin(t)

g=function(tt,xx,yy){
  
return(2*xx+yy-sin(tt))  
  
}

t=0
x=4.5
y=-6.5

x_vec=c()
y_vec=c()
t_vec=c()

ite=2000
h=0.01

for(l in 1:ite){
  
x_vec=c(x_vec,x)
y_vec=c(y_vec,y)
t_vec=c(t_vec,t)  
x=x+h*f(t,x,y)
y=y+h*g(t,x,y)
t=t+h

}

plot(t_vec,x_vec,xlab="t",ylab="x",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(x_vec),max(x_vec)),col=2)

plot(t_vec,y_vec,xlab="t",ylab="y",xlim=c(min(t_vec),max(t_vec)),ylim=c(min(y_vec),max(y_vec)),col=3)

Euler_x_vec=x_vec
Euler_y_vec=y_vec



```

```{r}


#p.122　渦の運動

#オイラー法　連立微分方程式

t=0
xA=1
yA=3
xB=2
yB=4

kB=-20
kA=10


xA_vec=c()
yA_vec=c()
xB_vec=c()
yB_vec=c()
t_vec=c()

ite=200
dt=0.1

for(l in 1:ite){

xA_vec=c(xA_vec,xA)
yA_vec=c(yA_vec,yA)
xB_vec=c(xB_vec,xB)
yB_vec=c(yB_vec,yB)
t_vec=c(t_vec,t)  


xA=xA-kB*(yA-yB)*dt/((xA-xB)^2+(yA-yB)^2)
yA=yA+kB*(xA-xB)*dt/((xA-xB)^2+(yA-yB)^2)
xB=xB+kA*(yA-yB)*dt/((xA-xB)^2+(yA-yB)^2)
yB=yB-kA*(xA-xB)*dt/((xA-xB)^2+(yA-yB)^2)
t=t+dt

}

#Aの表示
plot(xA_vec,yA_vec,xlab="x",ylab="y",xlim=c(min(c(xA_vec,xB_vec)),max(c(xA_vec,xB_vec))),ylim=c(min(c(yA_vec,yB_vec)),max(c(yA_vec,yB_vec))),col=2,main="渦の表示")

par(new=T)

#Bの表示
plot(xB_vec,yB_vec,xlab="x",ylab="y",xlim=c(min(c(xA_vec,xB_vec)),max(c(xA_vec,xB_vec))),ylim=c(min(c(yA_vec,yB_vec)),max(c(yA_vec,yB_vec))),col=3,main="渦の表示")



```


```{r}


#p.123　渦の運動

#オイラー法　連立微分方程式

t=0
xA=1
yA=3
xB=2
yB=4
xC=5
yC=6
xD=7
yD=8

kB=-20
kA=10
kC=5
kD=1


xA_vec=c()
yA_vec=c()
xB_vec=c()
yB_vec=c()
xC_vec=c()
yC_vec=c()
xD_vec=c()
yD_vec=c()
t_vec=c()

ite=200
dt=0.1

for(l in 1:ite){

xA_vec=c(xA_vec,xA)
yA_vec=c(yA_vec,yA)
xB_vec=c(xB_vec,xB)
yB_vec=c(yB_vec,yB)
xC_vec=c(xC_vec,xC)
yC_vec=c(yC_vec,yC)
xD_vec=c(xD_vec,xD)
yD_vec=c(yD_vec,yD)
t_vec=c(t_vec,t)  


xA=xA-kB*(yA-yB)*dt/((xA-xB)^2+(yA-yB)^2)-kC*(yA-yC)/((xA-xC)^2+(yA-yC)^2)-kD*(yA-yD)/((xA-xD)^2+(yA-yD)^2)
yA=yA-kB*(xA-xB)*dt/((xA-xB)^2+(yA-yB)^2)+kC*(xA-xC)/((xA-xC)^2+(yA-yC)^2)+kD*(xA-xD)/((xA-xD)^2+(yA-yD)^2)
xB=xB-kB*(yB-yA)*dt/((xB-xA)^2+(yB-yA)^2)-kC*(yB-yC)/((xB-xC)^2+(yB-yC)^2)-kD*(yB-yD)/((xB-xD)^2+(yB-yD)^2)
yB=yB-kB*(xB-xA)*dt/((xB-xA)^2+(yB-yA)^2)+kC*(xB-xC)/((xB-xC)^2+(yB-yC)^2)+kD*(xB-xD)/((xB-xD)^2+(yB-yD)^2)
xC=xC-kB*(yC-yB)*dt/((xC-xB)^2+(yC-yB)^2)-kC*(yC-yA)/((xC-xA)^2+(yC-yA)^2)-kD*(yC-yD)/((xC-xD)^2+(yC-yD)^2)
yC=yC-kB*(xC-xB)*dt/((xC-xB)^2+(yC-yB)^2)+kC*(xC-xA)/((xC-xA)^2+(yC-yA)^2)+kD*(xC-xD)/((xC-xD)^2+(yC-yD)^2)
xD=xD-kB*(yD-yB)*dt/((xD-xB)^2+(yD-yB)^2)-kC*(yD-yC)/((xD-xC)^2+(yD-yC)^2)-kD*(yD-yA)/((xD-xA)^2+(yD-yA)^2)
yD=yD-kB*(xD-xB)*dt/((xD-xB)^2+(yD-yB)^2)+kC*(xD-xC)/((xD-xC)^2+(yD-yC)^2)+kD*(xD-xA)/((xD-xA)^2+(yD-yA)^2)
t=t+dt

}

#Aの表示
plot(xA_vec,yA_vec,xlab="x",ylab="y",xlim=c(min(c(xA_vec,xB_vec,xC_vec,xD_vec)),max(c(xA_vec,xB_vec,xC_vec,xD_vec))),ylim=c(min(c(yA_vec,yB_vec,yC_vec,yD_vec)),max(c(yA_vec,yB_vec,yC_vec,yD_vec))),col=2,main="渦の表示")

par(new=T)

#Bの表示
plot(xB_vec,yB_vec,xlab="x",ylab="y",xlim=c(min(c(xA_vec,xB_vec,xC_vec,xD_vec)),max(c(xA_vec,xB_vec,xC_vec,xD_vec))),ylim=c(min(c(yA_vec,yB_vec,yC_vec,yD_vec)),max(c(yA_vec,yB_vec,yC_vec,yD_vec))),col=2,main="渦の表示")

par(new=T)

#Cの表示
plot(xC_vec,yC_vec,xlab="x",ylab="y",xlim=c(min(c(xA_vec,xB_vec,xC_vec,xD_vec)),max(c(xA_vec,xB_vec,xC_vec,xD_vec))),ylim=c(min(c(yA_vec,yB_vec,yC_vec,yD_vec)),max(c(yA_vec,yB_vec,yC_vec,yD_vec))),col=3,main="渦の表示")

par(new=T)

#Dの表示
plot(xD_vec,yD_vec,xlab="x",ylab="y",xlim=c(min(c(xA_vec,xB_vec,xC_vec,xD_vec)),max(c(xA_vec,xB_vec,xC_vec,xD_vec))),ylim=c(min(c(yA_vec,yB_vec,yC_vec,yD_vec)),max(c(yA_vec,yB_vec,yC_vec,yD_vec))),col=4,main="渦の表示")

```