第2章 振　動 ２．１　バネの運動 【リスト2-1】バネの運動 from tkinter import * import math,time Position=-100; SpringLength=200 Mass=10; Kfact=80; Myu=0.00 Velocity=0; DT=0.01; def execP(canvas): global Kfact,Position, SpringLength global Mass, Velocity, Myu # 計算 Force = -Kfact * Position Acc = Force / Mass * DT Velocity = (Velocity + Acc) * (1 - Myu) Position = Position + Velocity * DT Y=SpringLength + Position # 描画 canvas.delete("DT") canvas.create_line(100,50, 100,50+Y, fill="blue", tag="DT") canvas.create_oval(97, 47,103, 53, fill="yellow", tag="DT") canvas.create_oval(90, 40+Y,110, 60+Y, fill="red", tag="DT") #メイン処理 root=Tk() canvas=Canvas(root, width=500,height=500);canvas.pack() while(1): execP(canvas); root.update() time.sleep(0.001) 【リスト2-2】グラフ表示のための追加と変更 from tkinter import * import math,time IPos=Position;plotDataEnd=False; IT=0 Graph=[0 for i in range(5000)] def plotStart(root): global beforY,IT, SpringLength, IPos IT=0; beforY=50+SpringLength + IPos YY=50+ SpringLength Y1=YY+IPos*1.2;Y2=YY-IPos*1.2 ST="k = %5.2f μ = %7.4f m = %5.2f" % (Kfact, Myu, Mass) canvas.create_text(325,Y1-12, text=ST,font=('Times',12)) canvas.create_line(200,YY, 450,YY, fill="blue") for i in range(1,13): C="gray" if i % 5==0: C="black" Y=YY+i*10 canvas.create_line(200,Y, 450,Y, fill=C) Y=YY-i*10 canvas.create_line(200,Y, 450,Y, fill=C) for IT in range(5000): XX=IT*DT*5+200 if IT%1000==0: canvas.create_line(XX,Y1,XX, Y2, fill="gray") YY=Graph[IT]+50 canvas.create_line(XX,beforY, XX+1,YY, fill="red",width=2) beforY=YY;IT+=1 canvas.create_rectangle(200,Y1, 450,Y2,fill="", outline="blue", width=2) def plotData(Y): global Graph,IT, plotDataEnd Graph[IT]=Y IT+=1 if IT>=5000: plotDataEnd=True; return def execP(canvas): 　　・ 　　・(中略) 　　・ # 描画 　　plotData(Y) ・ 　 ・(中略) 　 ・ #メイン処理 root=Tk() canvas=Canvas(root, width=500,height=500);canvas.pack() while(1): if plotDataEnd: break execP(canvas); root.update() time.sleep(0.001) plotStart(root) ２．２　単振り子 【リスト2-5】振り子のシミュレーション from tkinter import * from tkinter import messagebox import time, math Myu=0.00; DT=0.01; Deg=-60 # Myu: 粘性抵抗, DT:時間刻み, Deg:初期角度 XX1=150; YY1=10 #■計算 def calc(): global V, TH, DT,Myu V=V-DT*(math.sin(TH)+V*Myu);TH=TH+DT*V if abs(V)+abs(TH)<0.01: return True return False #■振り子の描画 def drawLine(canvas): global XX1, YY1, XX2, YY2 canvas.delete("DT") XX2= XX1+100*math.sin(TH); YY2=YY1+100*math.cos(TH) canvas.create_line(XX1, YY1, XX2, YY2, fill="blue", tag="DT", width=1) canvas.create_oval(XX1-2, YY1-2, XX1+2, YY1+2, fill="yellow", tag="DT") canvas.create_oval(XX2-5, YY2-5, XX2+5, YY2+5, fill="red" , tag="DT") #■メイン処理 root=Tk(); root.title("振り子(pendulum)の運動") canvas=Canvas(root,width=500, height=500) canvas.pack() while True: TH=Deg*math.pi/180; V=0 # Radianによる角度と速度初期設定 while True: drawLine(canvas); root.update() if calc(): break time.sleep(0.002) if not messagebox.askokcancel("シミュレーション終了", "再開しいますか?"):break 【リスト2-6】グラフの表示(表示用のみ。他は【リスト2-5】と同じ) ・ ・(前略) ・ #グラフ表示用 IPos=math.pi*Deg/180;plotDataEnd=False; IT=0 Graph=[0 for i in range(5000)] def plotStart(root): global beforY,IT, IPos IT=0; beforY=200+IPos*180/math.pi; YY=200 Y1=YY-90;Y2=YY+90 canvas.create_text(325,Y1-12,text="μ = %7.4f" % Myu, font=('Times',12)) canvas.create_line(200,YY, 450,YY, fill="blue") for i in range(1,9): C="gray" if i % 5==0: C="black" Y=YY+i*10; canvas.create_line(200,Y, 450,Y, fill=C) Y=YY-i*10; canvas.create_line(200,Y, 450,Y, fill=C) for IT in range(5000): XX=IT*DT*5+200 if IT%1000==0: canvas.create_line(XX,Y1,XX, Y2, fill="gray") YY=Graph[IT]*180/math.pi+200 canvas.create_line(XX,beforY, XX+1,YY, fill="red",width=2) beforY=YY;IT+=1 canvas.create_rectangle(200,Y1, 450,Y2,fill="", outline="blue", width=2) def plotData(Y): global Graph,IT, plotDataEnd Graph[IT]=Y; IT+=1 if IT>=5000: plotDataEnd=True; return ・ ・(中略)関数calcは【リスト2-5】と同じ ・ def drawLine(canvas): global XX1, YY1, XX2, YY2 plotData(TH) 　　・ 　　・(中略)関数drawLineの後の部分は【リスト2-5】と同じ 　　・ #■メイン処理 ・ ・(中略)メイン処理の前部分は【リスト2-5】と同じ ・ while True: if plotDataEnd: break drawLine(canvas); root.update() if calc(): break time.sleep(0.002) plotStart(root);break 【リスト2-3】バネ形状の描画(描画部分のみ。他は【リスト2-1】と同じ) ・ ・（前略） ・ def drawSpring(canvas, L, X,Y, r, C): KM=6; X1 = X; Y1 = Y; M = 20; TH=math.pi; dth=TH/KM N = KM*2* M + KM+1; dy = L / (N + r); rr=r*0.5 X0=X1; Y0=Y1-rr*math.cos(TH); yH =Y0 for i in range(N): X2=r*math.sin(TH)+X0; Y2=rr*math.cos(TH)+yH canvas.create_line(X1,Y1, X2,Y2,fill=C,width=2) X1=X2; Y1=Y2; TH+=dth; yH+=dy canvas.create_line(X1,Y1, X1,L+50,fill=C,width=2)from tkinter import * Position=-100; SpringLength=200;Mass=10; Kfact=80; Myu=0.00 Velocity=0; DT=0.01; def execP(canvas): ・ ・（中略） ・ for i in range(20):#計算20回に1回描画する Force = -Kfact * Position Acc = Force / Mass * DT Velocity = (Velocity + Acc) * (1 - Myu) Position = Position + Velocity * DT canvas.delete(ALL) Y=SpringLength + Position drawSpring(canvas,Y-10, 100,50, 10, "black") canvas.create_oval(97, 47,103, 53, fill="yellow") canvas.create_oval(90, 40+Y,110, 60+Y, fill="red") ・ ・（後略） ・ 【リスト2-4】互い違いの線分によるバネ形状の表現 from tkinter import * import math,time SID=[]; BID=0 def drawSpring2Init(canvas,C="black", C2="red"): global SID,BID; M=20 SID=[0 for i in range(M*2+1)] for i in range(M*2+1): SID[i]=canvas.create_line(10,10, 20,20,fill=C,width=2) BID=canvas.create_oval(10, 10,30, 30, fill=C2) def updFig(canvas, L, X,Y, r, C, C2="yellow", C3="red"): global SID,BID X1 = X; Y1 = Y; M = 20; DY=L/(M*2);X11=X1-r; X12=X1+r Y2=Y1+DY/2;ID=0; canvas.coords(SID[ID], X1, Y1, X11,Y2) Y1=Y2 for i in range(M): Y2=Y1+DY; ID+=1;canvas.coords(SID[ID],X11,Y1, X12,Y2) Y1=Y2+DY;ID+=1; canvas.coords(SID[ID],X12,Y2, X11,Y1) canvas.coords(BID,X-10, Y+L-10,X+10, Y+L+10) Position=100; SpringLength=200 Mass=10; Kfact=40; Myu=0.0001;Velocity=0; DT=0.01; def execP2(canvas): global Kfact,Position, SpringLength, Mass, Velocity, Myu Force = -Kfact * Position; Acc = Force / Mass * DT Velocity = (Velocity + Acc) * (1 - Myu) Position = Position + Velocity * DT; Y=SpringLength + Position updFig(canvas, Y-10, 100,50, 10, "black") root=Tk() canvas=Canvas(root, width=500,height=500);canvas.pack() drawSpring2Init(canvas, "black") canvas.create_oval(100-3, 50-3,100+3, 50+3, fill="yellow") while(1):execP2(canvas); root.update(); time.sleep(0.001) ２．３　一次元波動 【リスト2-7】　１次元波動方程式 from tkinter import * import time import math def array(N1,N2): return [[0 for i in range(N2)] for j in range(N1)] Num=50 # 節点の数-1 Myu=0.001 # 粘性抵抗 M=2 # 質量 Kx=100 # X方向硬さ(ばね定数) Ky=0.005 # Y方向硬さ(バネ定数) DT=0.01 # 時間刻み DX=1 # X方向刻み幅 InitialType = 3 # Y座標値設定方法(この値を変えてみよう) FreeBoundary = False # 両端の境界条件(自由端のときTrue, 固定端のときFalse) DrawNode = True # 各節点を描画するときTrue，しないときFalse DrawSin = True # 近似用Sinカーブを描画するときTrue, しないときFalse X = array(2, Num+1); Y = array(2, Num+1) #　座標配列 VX = array(2, Num+1); VY = array(2, Num+1) #　速度配列 for i in range(Num+1):X[0][i]=i*DX # X座標値初期設定 X[1][0]=0; X[1][Num]=DX*Num # X座標の先頭と最後は変更しない if InitialType==1: DP = math.pi/Num; TH=0 # sin波 for i in range(Num+1):Y[0][i] = 8 * math.sin(TH); TH += DP elif InitialType==2: for i in range(Num+1): # 三角(弦を引っ張って離す/弾く) if i<=int(Num/2): Y[0][i]=16*i/Num else: Y[0][i]=16-16*i/Num elif InitialType==3: VY[0][int(Num/2)]=200 # 中央の初速度(弦を叩く) else : ii=int(Num/2);DP = math.pi/6; TH =-DP*2 #中央のみCos for i in range(ii-2,ii+3):Y[0][i]=-8*math.cos(TH); TH+=DP #■計算 def nextTime(i1): global FreeBoundary i2=1-i1; T=0 for i in range(1,Num): # X方向 D = X[i1][i-1] - 2 * X[i1][i] + X[i1][i + 1] A = D * M/Kx - VX[i1][i]*Myu VX[i2][i] = VX[i1][i] + A * DT X[i2][i] = X[i1][i] + VX[i2][i] * DT # X方向 D = Y[i1][i-1] - 2 * Y[i1][i] + Y[i1][i + 1] A = D * M/Ky-VY[i1][i]*Myu/DT VY[i2][i] = VY[i1][i] + A * DT Y[i2][i] = Y[i1][i] + VY[i2][i] *DT T+= abs(VY[i2][i])+abs(Y[i2][i]) if FreeBoundary: Y[i2][0]= Y[i2][1]; Y[i2][Num]= Y[i2][Num-1] if T<1: return 3 # 弦の停止を判定 return i2 #■弦の描画 def drawLine(canvas,i1,width=1): global DrawSin, DrawNode canvas.delete("DT") if DrawSin:# sin波近似曲線 DP = math.pi/Num; TH = 0; Y0 = Y[i1][int(Num/2)] for i in range(Num): # sin波近似曲線 YY1 = 100-10*Y0 * math.sin(TH); YY2 = 100-10 * Y0 * math.sin(TH+DP) XX1 = 10 + i * DX * 10; XX2=XX1+DX*10 canvas.create_line(XX1, YY1, XX2, YY2, fill="black", tag="DT", width=1) TH += DP for i in range(Num): #弦の描画 canvas.create_line(10+10 * X[i1][i] , 100-10*Y[i1][i] , 10+10 * X[i1][i+1], 100-10*Y[i1][i+1], fill="blue", tag="DT", width=2) if DrawNode:#弦の描画 for i in range(Num+1): XX = 10 + 10 * X[i1][i]; YY =100 - 10 * Y[i1][i] canvas.create_oval(XX-2, YY-2, XX+2, YY+2, fill="red", tag="DT") #■メイン処理 root=Tk() root.title("一次元の振動") WDH=10 + Num * DX * 10; canvas=Canvas(root,width=WDH+10, height=240);canvas.pack() i1=0 #for i in range(2): #決まった回数だけ実行するとき，次の行をコメント while True: drawLine(canvas,i1); root.update() i1=nextTime(i1) if i1==3: break time.sleep(0.001) ２．４　二次元波動 【リスト2-8】三次元グラフ(隠れ線処理) from tkinter import * import math,time #■初期設定 NumX = 50; NumY = 50 Z=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] ID=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] ZX=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] ZY=[[0 for j in range(NumY + 1)] for i in range(NumX + 1)] def initHidden(): XY=[[0,0] for i in range(4)] XY[0][0]=200;XY[0][1]=200; XY[1][0]=200;XY[1][1]=210 XY[2][0]=210;XY[2][1]=210; XY[3][0]=210;XY[3][1]=200 A=math.pi/9; B=math.pi/6; DX=5 cosA=DX*math.cos(A);cosB=DX*math.cos(B) sinA=DX*math.sin(A);sinB=DX*math.sin(B) for i in range(NumX+1): for j in range(NumY+1): ZX[i][j]=250+(i*cosA-j*cosB) ZY[i][j]=250-(i*sinA+j*sinB) for i in range(NumX): for j in range(NumY): ID[i][j]=canvas.create_polygon(XY,fill="#7777FF", outline="black") #■隠れ線処理 def hiddenLine(): global canvas, NumX, NumY XY=[0 for i in range(8)] for i in range(NumX): i2=NumX-i;i1=i2-1 for j in range(NumY): j2=NumX-j; j1=j2-1 XY[0]=ZX[i1][j1]; XY[1]=ZY[i1][j1]-Z[i1][j1] XY[2]=ZX[i1][j2]; XY[3]=ZY[i1][j2]-Z[i1][j2] XY[4]=ZX[i2][j2]; XY[5]=ZY[i2][j2]-Z[i2][j2] XY[6]=ZX[i2][j1]; XY[7]=ZY[i2][j1]-Z[i2][j1] canvas.coords(ID[i][j],XY) 　 【リスト2-9】二次元波動シミュレーション（【リスト2-8】の後に続けます） ・ ・(前略)【リスト2-8】のプログラム ・ ZS=[[[0 for j in range(NumY+1)] for i in range(NumX+1)] for k in range(2)] V =[[[0 for j in range(NumY+1)] for i in range(NumX+1)] for k in range(2)] simTime = 0; ID1 = 0; ID2 = 1; Myu=0.02;counter=0 FreeBoundary=False # 境界条件を自由端にするときTrueにする captMode= False #画面キャプチャするときTrueにする def setData(): #データを設定 for j in range(NumX+1): for k in range(NumY+1): Z[j][k]=80.0*ZS[ID1][j][k] def initWave(): global simTime, ID1, ID2, Myu,counter,XS,V simTime = 0; ID1 = 0; ID2 = 1; Myu=0.02;counter=0 for i in range(NumX+1): for j in range(NumY+1): DX = i - NumX/2; DY = j - NumY/2 R = math.sqrt(DX * DX + DY * DY) if R < 0.1: A=1 elif R > 3 : A=0 else : A = math.sin(R) / R; ZS[0][ i][ j] = -A; V[0][ i][ j] = 0; #print(i,j,ZS[0][ i][ j]) def waveEquation(SaveN, Alfa, DX, DT): global NumX, NumY, ZS, ID1,ID2,simTime,Myu,counter Beta = Alfa * Alfa / (DX * DX); for k in range(SaveN): for i in range(1, NumX): for j in range(1, NumY): ZS[ID2][i][j] = ZS[ID1][i][j] + DT * V[ID1][i][j] for i in range(1, NumX): for j in range(1, NumY): Acc = Beta * (ZS[ID2][i][j + 1] - 2 * ZS[ID2][i][j] + ZS[ID2][i][j - 1] + ZS[ID2][i + 1][j] - 2 * ZS[ID2][i][j] + ZS[ID2][i - 1][j]) V[ID2][i][j] = V[ID1][i][j] + DT *(Acc- Myu*V[ID1][i][j]) if FreeBoundary: #自由端境界条件 for i in range(NumX+1): V [ID2][i][0] = V [ID2][i][1] V [ID2][i][NumY] = V [ID2][i][NumY-1] ZS[ID2][i][0] = ZS[ID2][i][1] ZS[ID2][i][NumY] = ZS[ID2][i][NumY-1] for i in range(NumY+1): V [ID2][0][i] = V [ID2][1][i] V [ID2][NumX][i] = V [ID2][NumX-1][i] ZS[ID2][0][i] = ZS[ID2][1][i] ZS[ID2][NumX][i] = ZS[ID2][NumX-1][i] else: #固定端境界条件 for i in range(NumX+1): V[ID2][i][0] =0; V[ID2][i][NumY]= 0 ZS[ID2][i][0]=0; ZS[ID2][i][NumY]=0 for i in range(NumY+1): V[ID2][0][i] =0; V[ID2][NumX][i]=0 ZS[ID2][0][i]=0; ZS[ID2][NumX][i] = 0 simTime += DT; ID1 = ID2; ID2 = 1-ID2;counter+=1 return 0 def leftMouseDown(event):#ダウン waveEquation(50, 0.05, 0.05, 0.01) setData(); hiddenLine(); root.title("2次元波動 t = %f" % simTime);root.update() #■メイン処理 root=Tk();root.title("2次元波動");root.geometry("500x300") canvas=Canvas(root,width=500,height=300,bg="white"); canvas.pack() initHidden(); initWave() if captMode: waveEquation(50, 0.05, 0.05, 0.01) setData(); hiddenLine(); root.title("2次元波動 t=%f" % simTime);root.update() canvas.bind("",leftMouseDown) root.mainloop() else: while(1): waveEquation(50, 0.05, 0.05, 0.01) setData(); hiddenLine() root.title("2次元波動 t = %f" % simTime);root.update() time.sleep(0.01)