1 a波の動き

　波の動きをシミュレーションすることを考える。このレポートでは，ある波に対して他の波を生成，衝突させ，全体的な波を抑制させることを目的とする。
初めに，波をシミュレーションする方法をまとめる。一般的な波は以下の波動方程式を満たすことが知られている。

\[ 　\frac{\partial^2 u}{\partial t^2}(x,y,t) = \bigg(\frac{ \tau}{p} \bigg)^2\bigg(\frac{\partial^2 u}{\partial x^2}(x,y,t) + \frac{\partial^2 u}{\partial y^2}(x,y,t)\bigg) \\ 　\mbox{但し，$\tau$は膜の張力，$p$は面密度} \]

この微分方程式を数値解析的に解き，水の波のシミュレーションを作成する。ここでは解き方として差分法を採用する。以下を初期設定である。

N <- 30 #x,yの範囲 (m)
tau <- 72.75 #水の表面張力 (N/m)
p <- 1000 #面密度 (kg/m^2)
dt <- 0.1 #時間増分 (s)
dx <- 1 #x分割幅 (m)
dy <- 1 #y分割幅 (m)

N1 <- N%/%3
c <- tau/p
num <- N/dx 
num2 <- (num)^2 #総plot数

#x,y座標のデータ
Hx <- (1:num2)
Hy <- (1:num2)
k <- 1
for(i in 1:num){
  Hy[(num*(i-1)+1):(num*i)] <- i*dy
  for(j in 1:num){
   Hx[k] <- j*dx
   k <- k+1
  }
}

初めに波の初期状態を設定する。また，ここでは高さ30mの回転放物面を初期波とし観察を行う。以下，その回転放物面を設定するためのスクリプトである。

u1 <- rep(0, num2) #データを保存1
u2 <- rep(0, num2) #データを保存2
u3 <- rep(0, num2) #データを保存3

#初期波のデータを保存
O <- c(N/2, N/2) #初期波の進行方向

for(i in 1:num2){ 
  if(-3*(Hx[i]-O[1])^2-3*(Hy[i]-O[2])^2+30 > 0){
    u1[i] <- u1[i]-3*(Hx[i]-O[1])^2-3*(Hy[i]-O[2])^2+30
  }
}

以下が差分法を用いた波動方程式の数値解析である。また，この解析法では誤差が生じる。故に，ここでは「波が収まる」を以下のように定義する。:
試行によって得られた波の状態をu3とする。その試行においてu3の最も大きい値をmax，最も小さい値をminとし，これらの差heightが1(m)を下回ったとき「波が収まる」とする。
関数sabunは差分法を実行するための関数で，第1引数として初期波の情報を持つu1，第2引数として1か2を与え，第3引数として終了する高さ[m]または時間[s]を指定する値を与える。これらで差分法を止める条件を指定することができ，1ならば最も大きい高低差heightがTendで指定した高さに収まるまで試行を繰り返し，2ならば時間がTendで指定した時間に達したとき停止する。第4引数として0か1を与える。波の衰退を考慮するときは0，しないときは1を与える。第5引数には他に加えたい機能を関数として与える。即存の機能としては，差分法の他に中央での高さと端での高さの平均を測定するためのスクリプトが含まれており，戻り値として波が収まった時間とそれらの測定した高さをlistにまとめて返す。

# 引数は初期波u1
# a: 1 -> height, 2 -> time
# Tend: height > Tend || t < Tend
# dec : 0 -> 波の衰退あり, 1 -> 波の衰退なし
# Fun : Fun <- list(list(関数1, n1), list(関数2, n2)... n: 0 -> u1,u2に値を返す, 1 -> u1, u2に値を返さない

sabun <- function(u1, a, Tend, dec, Fun){
  t <- 0 #経過時間 (s)
  w <- 0 #試行回数
  u2 <- u1 #u1をu2 に複製する
  height <- 3
  max <- 0
  min <- 1
  O1h <- c() #中心の高さを収納
  Ave <- c() #端での高さの平均を収納
  X <- c()
  Y <- c()
  
  if(a == 1){
    b <- height > Tend
  }else{
    b <- t < Tend
  }
  
  while(b){
    height <- 0
    max <- 0
    min <- 1
    ave <- 0
    
    #加えたい機能を処理する
    if(length(Fun) != 0){
      for(i in 1:length(Fun)){
        if(Fun[[i]][[2]] == 0){
          F <- Fun[[i]][[1]](u1, u2, w, t)
          u1 <- F[[1]]
          u2 <- F[[2]]
        }else{
         Fun[[i]][[1]](u1, u2, w, t)
        }
      }
    }
    
    for(i in 1:num2){ 
      #端の4点を除く点の動きを計算
      if((Hx[i] == dx && Hy[i] == dy) || (Hx[i] == N && Hy[i] == dy) || (Hx[i] == dx && Hy[i] == N) || (Hx[i] == N && Hy[i] == N)){
        
      }else if(Hx[i] == dx){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((u2[i+1]-2*u2[i])/dx^2+(u2[i+num]-2*u2[i]+u2[i-num])/dy^2)
        
      }else if(Hx[i] == N){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((-2*u2[i]+u2[i-1])/dx^2+(u2[i+num]-2*u2[i]+u2[i-num])/dy^2) 
        
      }else if(Hy[i] == dy){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((u2[i+1]-2*u2[i]+u2[i-1])/dx^2+(u2[i+num]-2*u2[i])/dy^2)
        
      }else if(Hy[i] == N){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((u2[i+1]-2*u2[i]+u2[i-1])/dx^2+(-2*u2[i]+u2[i-num])/dy^2) 
        
      }else{ 
        u3[i] <- 2*u2[i] -u1[i] + dt*c^2*((u2[i+1]-2*u2[i]+u2[i-1])/dx^2 + (u2[i+num]-2*u2[i]+u2[i-num])/dy^2)  
      }
      
      #端の4点の動きを計算
      if(Hx[i] == dx && Hy[i] == dy){　
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((u2[i+1]-2*u2[i])/dx^2+(u2[i+num]-2*u2[i])/dy^2)
        
      }else if(Hx[i] == N && Hy[i] == dy){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((-2*u2[i]+u2[i-1])/dx^2+(u2[i+num]-2*u2[i])/dy^2)
        
      }else if(Hx[i] == dx && Hy[i] == N){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((u2[i+1]-2*u2[i])/dx^2+(-2*u2[i]+u2[i-num])/dy^2)
        
      }else if(Hx[i] == N && Hy[i] == N){
        u3[i] <- 2*u2[i]-u1[i]+dt*c^2*((-2*u2[i]+u2[i-1])/dx^2+(-2*u2[i]+u2[i-num])/dy^2)
      }
      
      #波の衰退を計算
      if(dec != 1){
        u3[i] <- u3[i]*0.999
      }
      
      #高低差の計算
      if(u3[i] > max){
        max <- u3[i]
      }
      
      if(u3[i] < min){
        min <- u3[i]
      }
      
      #端での波の高さの平均を計算
      if(Hx[i] == dx || Hx[i] == N || Hy[i] == dy || Hy[i] == N){
        ave <- ave + u3[i]
      }
    }
    
    height <- abs(max - min)
    u1 <- u2
    u2 <- u3
    t <- t+dt
    w <- w+1
    O1h <- append(O1h, u3[num2/2-N/2])
    Ave <- append(Ave, ave)
    
    if(a == 1){
      b <- height > Tend
    }else{
      b <- t < Tend
    }
  }
  return(list(t, O1h, Ave))
}

以上のスクリプトを使い解析を進める。
また，初期波は先ほど定義した回転放物面とする。
関数sabunに3次元グラフを作成するscatterplot3を実行する機能を加える。そのための関数は以下のスクリプトである。

scat <- function(u1, u2, w, t){
   if(w %% 20 == 0){
     scatterplot3d(Hx, Hy, u1, xlim=c(1, N), ylim=c(1, N), zlim=c(-30, 30), xlab="width [m]", ylab="depth [m]",  zlab="height [m]", angle=100, sub=round(t,1))
   }
}

Fun <- list(list(scat, 1))

初期波は以下の様な曲面になる。

plot3d(Hx, Hy, u1, xlim=c(1, N), ylim=c(1, N), zlim=c(-30, 30), xlab="width [m]", ylab="depth [m]",  zlab="height [m]")

You must enable Javascript to view this page properly.

また，今回は波が収まるまで試行を行うようにするため，第2引数として1，第3引数として1[m]を与える。第5引数としては上記で作成したFunを与える。波が収まるまで以下のような動きをする。

kihon <- sabun(u1, 1, 1, 0, Fun)

収まるまで次の時間が掛かった。

print(kihon[[1]])

## [1] 432.5

以下は，中央の高さと端における高さの平均をとったグラフである。また，時間は300秒までとする。

plot(1:length(kihon[[2]])*0.1,　kihon[[2]], xlim=c(1, 300), ylim=c(-45, 45), lab=c(30, 14, 0), type="l", xlab="time [s]", ylab="height [m]",sub="中央の高さ (回転放物面)")

plot(1:length(kihon[[3]])*0.1, kihon[[3]], xlim=c(1, 300), ylim=c(-80, 80), lab=c(30, 14, 0), type="l", xlab="time [s]", ylab="height [m]", sub="端の高さの平均 (回転放物面)")

2 b波の性質

2.1 b1周期性

　前セクションにて作成した回転放物面による中央の高さと端における高さの平均をとったグラフより，上下運動においては完全ではないが周期性が存在することが分かる。ここでは４角錐及び４面体をそれぞれ初期波とし，高さに関して周期性が存在するのかを検討する。

まず，４角錐を扱う。以下が４角錐を設定するためのスクリプトである。関数heiは引数として３次元座標Z1，Z2，Z3，Z4を受けとる。Z1，Z2，Z3を通るような平面を用意し，Z4のx,y座標における平面上のz座標を返す関数である。

hei <- function(Z1, Z2, Z3, Z4){
  z12 <- Z2-Z1
  z13 <- Z3-Z1
  
  hou <- c(z12[2]*z13[3]-z12[3]*z13[2], z12[3]*z13[1] - z12[1]*z13[3], z12[1]*z13[2]-z12[2]*z13[1])
  
  return ((-hou[1]*(Z4[1]-Z1[1])-hou[2]*(Z4[2]-Z1[2]))/hou[3]+Z1[3])#z座標を出力
}

Z1 <- c(8, 23, 0)
Z2 <- c(23, 23, 0)
Z3 <- c(8, 8, 0)
Z4 <- c(23, 8, 0)
Z5 <- c(15, 15, 20)

for(i in 1:num2){ #初期の波のデータを保存
  if(Hx[i] > 7 && Hx[i] < 24 && Hy[i] > 7 && Hy[i] < 24){
    if(Hx[i]==15 && Hy[i]==15){
      u1[i] <- 20
    }else if(Hy[i] >= Hx[i] && Hy[i] >= -Hx[i]+31){
      Z <- c(Hx[i], Hy[i]) 
      u1[i] <- hei(Z1, Z2, Z5, Z)
      
    }else if(Hy[i] <= Hx[i] && Hy[i] >= -Hx[i]+31){
      Z <- c(Hx[i], Hy[i]) 
      u1[i] <- hei(Z2, Z4, Z5, Z)
      
    }else if(Hy[i] <= Hx[i] && Hy[i] <= -Hx[i]+31){
      Z <- c(Hx[i], Hy[i]) 
      u1[i] <- hei(Z4, Z3, Z5, Z)   
      
    }else if(Hy[i] >= Hx[i] && Hy[i] <= -Hx[i]+31){
      Z <- c(Hx[i], Hy[i]) 
      u1[i] <- hei(Z1, Z3, Z5, Z) 
      
    }else{
      u1[i] <- 0
    }
  }
}

初期波は以下の様な曲面になる。

plot3d(Hx, Hy, u1, xlim=c(1, N), ylim=c(1, N), zlim=c(-30, 30), xlab="width [m]", ylab="depth [m]",  zlab="height [m]")

You must enable Javascript to view this page properly.

これらを用いて高さに関するグラフを作成する。但し，ここでは関数sabunをそのまま使用する。

Fun <- list()
sikakusui <- sabun(u1, 1, 1, 0, Fun)

以下は，中央の高さと端における高さの平均をとったグラフである。また，時間は300秒までとする。

plot(1:length(sikakusui[[2]])*0.1,　sikakusui[[2]], xlim=c(1, 300), ylim=c(-45, 45), lab=c(30, 14, 0), type="l", xlab="time [s]", ylab="height [m]", sub="中央の高さ (４角錐)")

plot(1:length(sikakusui[[3]])*0.1, sikakusui[[3]], xlim=c(1, 300), ylim=c(-80, 80), lab=c(30, 14, 0), type="l", xlab="time [s]", ylab="height [m]", sub="端の高さの平均 (４角錐)")

次に４面体を以下のように設定して，同様にしてグラフを作成する。

Z1 <- c(O[1], O[2]-14/3*sqrt(3), 0)
Z2 <- c(O[1]-7, O[2]+7/3*sqrt(3), 0)
Z3 <- c(O[1]+7, O[2]+7/3*sqrt(3), 0)
Z4 <- c(O[1], O[2], 14/3*sqrt(6))

#傾きを計算
a1 <- (Z1[2]-Z2[2])/(Z1[1]-Z2[1])
a2 <- (Z1[2]-Z3[2])/(Z1[1]-Z3[1])
a3 <- (Z4[2]-Z2[2])/(Z4[1]-Z2[1])
a4 <- (Z4[2]-Z3[2])/(Z4[1]-Z3[1])


for(i in 1:num2){ #初期の波のデータを保存
  if(Hy[i] > a1*(Hx[i]-Z1[1])+Z1[2] && Hy[i] > a2*(Hx[i]-Z1[1])+Z1[2] && Hy[i] < Z2[2]){
    if(Hx[i]==Z4[1] && Hy[i]==Z4[2]){
      u1[i] <- 14/3*sqrt(6)
      
    }else if(Hx[i] <= O[1]){
      if(Hy[i] < a3*(Hx[i]-O[1])+O[2]){
        Z <- c(Hx[i], Hy[i]) 
        u1[i] <- hei(Z1, Z2, Z4, Z)
      }else{
        Z <- c(Hx[i], Hy[i]) 
        u1[i] <- hei(Z2, Z3, Z4, Z)
      }
      
    }else{
      if(Hy[i] < a4*(Hx[i]-O[1])+O[2]){
        Z <- c(Hx[i], Hy[i]) 
        u1[i] <- hei(Z1, Z3, Z4, Z)
      }else{
        Z <- c(Hx[i], Hy[i]) 
        u1[i] <- hei(Z2, Z3, Z4, Z)
      }
    }
    
  }else{
      u1[i] <- 0
  }
}

plot3d(Hx, Hy, u1, xlim=c(1, N), ylim=c(1, N), zlim=c(-30, 30), xlab="width [m]", ylab="depth [m]",  zlab="height [m]")

You must enable Javascript to view this page properly.

simen <- sabun(u1, 1, 1, 0, Fun)

plot(1:length(simen[[2]])*0.1,　simen[[2]], xlim=c(1, 300), ylim=c(-45, 45), lab=c(30, 14, 0), type="l", xlab="time [s]", ylab="height [m]", sub="中央の高さ (４面体)")

plot(1:length(simen[[3]])*0.1, simen[[3]], xlim=c(1, 300), ylim=c(-80, 80), lab=c(30, 14, 0), type="l", xlab="time [s]", ylab="height [m]", sub="端の高さの平均 (４面体)")

以上の結果より，この２つの初期波においても上下運動において周期性が存在することが分かった。さらに，高さについては周期性は確認できないが緩やかな増減である。よって，緩やかな増減であることに注意すれば，周期を半分ずらした同じ初期波を持つ波を衝突させることにより，波の高さが相殺され抑制することが可能である。

2.2 b2独立性

　次に，波が持つ「独立性」について検証を行う。波の独立性とは二つの波を衝突させたとき，直後には二つの波は合わさり大きな波を形成するが，時間が経つとそれぞれの波形に戻る性質のことである。
ここでは前セクションで扱った回転放物面を初期波にもつ波を2つそれぞれ２点の地点(5,15)，(25,15)で同時に発生させ，15秒後中心周辺にある波を残し他の波を消す作業行う。また，平面y=15の断面の経過を動画に起こす。
関数sabunに平面y=15における断面をplotする機能入れるための関数docを引数として与え，波の衰退が起きないように設定し，さらに，70秒で終了するように設定する。動画は波発生後17秒経過後のグラフである。

u1 <- rep(0, num2)

O1 <- c(5, N/2)
O2 <- c(N-5, N/2)

for(i in 1:num2){ 
  if(-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30 > 0){
    u1[i] <- u1[i]-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30
  }
  
  if(-3*(Hx[i]-O2[1])^2-3*(Hy[i]-O2[2])^2+30 > 0){
    u1[i] <- u1[i]-3*(Hx[i]-O2[1])^2-3*(Hy[i]-O2[2])^2+30
  }
}

#機能と加えるためのスクリプト
doc <- function(u1, u2, w, t){
  
  X <- c()
  Y <- c()
  
  for(i in 1:num2){
   if(Hy[i]==15){
     X <- append(X, Hx[i])
     Y <- append(Y, u1[i])
   }
  }
  
  #17秒後記録を開始
  if(round(t*10)/10 > 17){
    if(w %% 6 == 0){
      plot(X, Y, xlim=c(0, 30), ylim=c(-10, 10), xlab="width [m]", ylab="height [m]", sub=round(t,1), type="l")
    }
  }
  
  if(round(t*10)/10 == 15){
    for(i in 1:num2){
      if(u1[i] < 0 ){
        u1[i] <- 0
        u2[i] <- 0
      }
  
      if(Hx[i] < 5 || Hy[i] < 5 || Hx[i] > 25 || Hy[i] > 25){
        u1[i] <- 0
        u2[i] <- 0
      }
    }
  }
  return(list(u1, u2))
}

Fun <- list(list(doc, 0))

ind <- sabun(u1, 2, 70, 1, Fun)

以上より，波の独立性を確認することができた。
これより，ある波に対して異なる波を衝突させ瞬間的に相殺することができたとしても，独立性により再び波が発生することが分かる。よって，完全に相殺させるためには同じような周期を持つ波を衝突させる必要がある。

3 c抑制の応用

ここでは前セクションで得られた結論を参考にし，ある波に対して初期波の半分を持つ波を衝突させたとき同じように抑制することが可能かどうかを観察する。また，抑制したい波の初期波はセクション1で設定した回転放物面とする。この半分の波は以下のような局面である。

u1 <- rep(0, num2)
u2 <- rep(0, num2)

O1 <- c(O[1]+1, O[2])     
for(i in 1:num2){ 
  if(-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30 > 0 && Hx[i] <= O1[1] ){
    u1[i] <- u1[i]-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30
    u2[i] <- u2[i]-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30
  }
}

plot3d(Hx, Hy, u1, xlim=c(1, N), ylim=c(1, N), zlim=c(-30, 30), xlab="width [m]", ylab="depth [m]",  zlab="height [m]")

You must enable Javascript to view this page properly.

セクション１の最後に書いた中央での高さのグラフより，発生から１０秒経過後，波は最も低い値を持つことが分かる。よって，その時間に合わせて初期波が半分の波を中央から1mずらした場所に発生させる。
関数sabunに衝突させる波を10秒経過後に発生させる機能を加える。

u1 <- rep(0, num2) 

#初期波を設定する
for(i in 1:num2){ 
  if(-3*(Hx[i]-O[1])^2-3*(Hy[i]-O[2])^2+30 > 0){
    u1[i] <- u1[i]-3*(Hx[i]-O[1])^2-3*(Hy[i]-O[2])^2+30
  }
}

#機能と加えるためのスクリプト
syou <- function(u1, u2, w, t){
  #10秒後に衝突波を発生させる
  if(round(t*10)/10 == 10){
    O1 <- c(O[1]+1, O[2])   
    for(i in 1:num2){ 
      if(-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30 > 0 && Hx[i] <= O1[1] ){
        u1[i] <- u1[i]-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30
        u2[i] <- u2[i]-3*(Hx[i]-O1[1])^2-3*(Hy[i]-O1[2])^2+30
      }
    }
  }
  return(list(u1, u2))
}

Fun <- list(list(syou, 0))

sab <- sabun(u1, 1, 1, 0, Fun)

波が収まるまで次のような時間が掛かった。

print(sab[[1]])

## [1] 393.1

セクション1において計測した時間は

print(kihon[[1]])

## [1] 432.5

であったから，波が収まるまでの時間が短くなっていることが分かる。次に高さに関するグラフを作成する。時間は同様に300秒までとし，セクション1のグラフにそれぞれ重ねて表示することにする。

plot(1:length(sab[[2]])*0.1, sab[[2]], xlim=c(1, 300), ylim=c(-45, 45), lab=c(30, 14, 0), type="l", col="blue", xlab="time [s]", ylab="height [m]", sub="中央の高さ (回転放物面と衝突波)")
par(new=T)
plot(1:length(kihon[[2]])*0.1, kihon[[2]], xlim=c(1, 300), ylim=c(-45, 45), lab=c(30, 14, 0), type="l", col="red", xlab="time [s]", ylab="height [m]", sub="中央の高さ (回転放物面と衝突波)")
legend(250, 40, legend=c("衝突波なし", "衝突波あり"),col=c("red","blue"),lty=c(1,1))#, pch=c(1,2))

plot(1:length(sab[[3]])*0.1, sab[[3]], xlim=c(1, 300), ylim=c(-80, 80), lab=c(30, 14, 0), type="l", col="blue", xlab="time [s]", ylab="height [m]", sub="端の高さの平均 (回転放物面と衝突波)")
par(new=T)
plot(1:length(kihon[[3]])*0.1, kihon[[3]], xlim=c(1, 300), ylim=c(-80, 80), lab=c(30, 14, 0), type="l", col="red", xlab="time [s]", ylab="height [m]", sub="端の高さの平均 (回転放物面と衝突波)")
legend(250, 75, legend=c("衝突波なし", "衝突波あり"),col=c("red","blue"),lty=c(1,1))#, pch=c(1,2))

以上のグラフより，初期波が半分である衝突波によって，中央の高さは全体的に半分程度の高さに抑制されており，端の高さにおいても後半は同じような結果が得られていることが分かる。よって，波を抑制するには，同じ周期，波形を持つ波を衝突させるのが効果的であることが検証できた。

波の衝突による抑制

宮野朋浩

2016年1月19日

1 a波の動き

2 b波の性質

2.1 b1周期性

2.2 b2独立性

3 c抑制の応用