■５．１　線形補間と多項式補間 表5-1 sin関数表 X(度） sin(X) 0 0.00000 10 0.17365 20 0.34202 30 0.50000 40 0.64279 表5-2 多項式補間プログラム def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] def copyArray(A):#配列複写 N1=len(A); N2=len(A[0]); C=array(N1,N2) for i in range(N1): for j in range(N2): C[i][j]=A[i][j] return C def GaussianELMPivot(MAT): A=copyArray(MAT);N=len(A);N1=N+1;ESP=0.0000001 for k in range(N):# 前進消去 k1 = k + 1;Max = abs(A[k][k]);IR = k #----ピボット入れ替え if k != N : for i in range(k1,N): ab=abs(A[i][k]) if ab>= Max : Max = ab; IR = i if Max < ESP : print("解を求めることができません") return [] if IR != k : D=A[k]; A[k]=A[IR]; A[IR]=D #---ピボット入れ替えここまで for i in range(k + 1,N): ALFA = A[i][k] / A[k][k] for j in range(k + 1,N1):A[i][j]-=ALFA*A[k][j] B=array(N,1) for i in range(N-1, -1, -1):# 後退代入 T = A[i][N] for j in range(i + 1,N):T -=A[i][j] * A[j][N] A[i][N] = T / A[i][i]; B[i][0]=A[i][N] return B #■以下，補間用プログラム def setParam(XT,YT): N=len(XT); N1 = N+1; A=array(N,N1) for i in range(N): X = XT[i]; A[i][0]=1 for j in range(1,N): A[i][j] = A[i][j-1]*X A[i][N]= YT[i] return GaussianELMPivot(A) def Interpolate(X,Param): N = len(Param); XX = 1;T = 0 for i in range(N): T += XX*Param[i][0]; XX*=X return T XT=[0,10,20,30,40] YT=[0, 0.17365, 0.34202, 0.5, 0.64279] Param=setParam(XT,YT);print("係数リスト", Param) print("補間結果",Interpolate(22,Param)) ■５．２　グレゴリー・ニュートンの補間式 表5-3 グレゴリー・ニュートンの補間式による補間 #■以下，補間用プログラム def Comb(N, R): #組合せ回数の計算 if R == 0 or R == N : return 1 if R == 1 : return N return Comb(N-1, R-1) + Comb(N-1, R) def setParam(F,UX):#係数の設定 a(i) N=len(F);A=array(N) for i in range(N): C=F[i];SG=-1 for j in range(i-1,0, -1): CB=Comb2(i, j); C += SG*CB*F[j];SG = -SG for i in range(1,i+1): C /= i # i! で除す処理 A[i]=C return A def GregoryTest(F, UX): A=setParam(F,UX);N=len(A) T=0 #各項目の計算とトータルの計算 for i in range(N): D=UX; FX=A[i] for j in range(i): FX *= D; D -= 1 T += FX; D -= 1 return T XT=[0,10,20,30,40] YT=[0, 0.17365, 0.34202, 0.5, 0.64279] X=GregoryTest(YT, 22/10) print("補間結果",X) 表5-4 組合せ個数の効率的な計算法 def Comb2(N,R):#組合せ回数計算の改善 A=array(N);NN=N; RR=R if (NN - RR) < RR : RR = N - R if RR == 0 : return 1 if RR == 1 : return N for i in range(2, RR+1): A[i-1 ]= i +1 for i in range(1, NN - R): A[0] = i +2 for j in range(2, RR+1): A[j-1] = A[j - 2] + A[j-1] return A[RR-1] ■５．３　ラグランジェの補間式 表5-5 ラグランジェの補間プログラム import math def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] #■以下，補間用プログラム def setTerm(i, N, F):#項設定 ii=1 for j in range(1,N-i+1): ii*=-1 T=ii * F[i-1] #T=ii * F[i-1] / Fact(i - 1) / Fact(N - i) for k in range(1, i):T/=k #階乗の計算でオーバフロー for k in range(1, N-i+1):T/=k #するため，この計算で代替 return T def interp(i, N, UX):#補間計算 U = UX; T = 1 for k in range(1,N+1): if i !=k : T *=U U -= 1 return T def LagrangeInter(F, UX): N=len(F); T=-F[0] for i in range(1,N+1): A = F[1-1] if i >= 2: A=setTerm(i, N, F) * interp(i, N, UX) T = T + A return T YT=[0, 0.17365, 0.34202, 0.5, 0.64279] X=LagrangeInter(YT, 22/10) print("補間結果",X) 表5-6 ロジスティックの成長曲線の近似 from graph import * #graph.pyについては表4-1参照 ・ ・ ・ ---(中略：その他の関数については表5－5参照)--- ・ ・ ・ def LogisticV(k,m,a,t):return k /(1 + m * math.exp(-a*t)) def setLogisticData(): Y=array(101);k=1; m=100; a=1.5; t=0 for i in range(101): Y[i]=LogisticV(k,m,a,t);t+=0.1 return Y def selectDat(A): N=int(len(A)/10+1); F=[] for i in range(N): F.append(A[i*10]) return F def loopLagrangeInter(F): N=len(F)*10+1; A=[];t=0 for i in range(N): X= LagrangeInter(F, t); A.append(X); t +=0.1 return Aroot=initTk("Numerical Analysis") canvas=initCanvas(root,440,240);canvas.pack() drawScale(canvas, 30,20, 15, 10, 30, 10, 0, -0.5, 1, 0.1, 5, 5, 1, 5, 10, 20,"%d","%.1f") A=setLogisticData();F=selectDat(A);X=loopLagrangeInter(F):t=0 for i in range(1, len(A)): drawV(canvas, t,X[i-1],t+0.1,X[i],"red", X0=30,Y0=170,Xscale=30,Yscale=100) drawV(canvas, t,A[i-1],t+0.1,A[i],"blue", width=2, X0=30,Y0=170,Xscale=30,Yscale=100) t+=0.1 表5-7 誤差の明確化 ---(前略：その他の部分については表5－6参照)--- drawScale(canvas, 30,20, 15, 10, 30, 10, 0, -0.05, 1, 0.01, 5, 5, 1, 5, 10, 20,"%d","%.2f") A=setLogisticData();F=selectDat(A);X=loopLagrangeInter(F); t=0 for i in range(1, len(A)): drawV(canvas, t,A[i-1]- X[i-1],t+0.1,A[i]- X[i],"red", X0=30,Y0=160,Xscale=30,Yscale=1000) t+=0.1 【データ生成用】 def setData1(): Y=array(101) for i in range(101): if i<=30: Y[i]=10 elif i>70: Y[i]=5 else :Y[i]=-5*(i-30)/40+10 return Y def setData2(): Y=array(101) for i in range(101): if i<=50: Y[i]=i*0.02+5 else: Y[i]=(i-50)*0.02+5 print(Y[i]) return Y ■５．４　スプライン補間 表5-8 Bスプライン基底関数を求めるプログラム例 from graph import * import math def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] #■以下，Bスプライン基底関数 def setQ(): Q=array(6) for i in range(6): Q[i] = (i+1) * 0.2 return Q def Bspline(Q, i, K, X): # Kは次数+1で与えます if K == 1: if X < Q[i-1]: return 0 if X < Q[i]: return 1 return 0 B1 = Bspline(Q, i, K - 1, X) B2 = Bspline(Q, i + 1, K - 1, X) return ((X - Q[i-1]) * B1 / (Q[i+K-2] - Q[i-1]) + (Q[i+K-1] - X) * B2 / (Q[i+K-1] - Q[i])) def setBsplineBase(N):#基底関数を求める（引数は次数） Q=setQ(); L=[] for i in range(121): X=i*0.01;Y=Bspline(Q, 1, N+1, X); L.append([X,Y]) return L root=initTk("Numerical Analysis") canvas=initCanvas(root,440,240);canvas.pack() drawScale(canvas, 30,20, 15, 10, 30, 20, 0, 0, 0.1, 0.1, 5, 5, 1, 5, 12, 10,"%.1f","%.1f") L=setBsplineBase(4) #引数に次数を指定する for i in range(1, len(L)): drawV(canvas, L[i-1][0],L[i-1][1],L[i][0],L[i][1],"red", X0=30,Y0=220,Xscale=300,Yscale=200) 表5-9 スプライン補間 from graph import * import math def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] def copyArray(A): N1=len(A); N2=len(A[0]); C=array(N1,N2) for i in range(N1): for j in range(N2): C[i][j]=A[i][j] return C def GaussianELMPivot(MAT):#連立方程式の解法 A=copyArray(MAT);N=len(A);N1=N+1;ESP=0.0000001 for k in range(N):# 前進消去 k1 = k + 1;Max = abs(A[k][k]);IR = k #----ピボット入れ替え if k != N : for i in range(k1,N): ab=abs(A[i][k]) if ab>= Max : Max = ab; IR = i if Max < ESP : print("解を求めることができません") return [] if IR != k : D=A[k]; A[k]=A[IR]; A[IR]=D #---ピボット入れ替えここまで for i in range(k + 1,N): ALFA = A[i][k] / A[k][k] for j in range(k + 1,N1):A[i][j]-=ALFA*A[k][j] B=array(N,1) for i in range(N-1, -1, -1):# 後退代入 T = A[i][N] for j in range(i + 1,N):T -=A[i][j] * A[j][N] A[i][N] = T / A[i][i]; B[i][0]=A[i][N] return B #■以下，Bスプライン関連 ・ ・ ・ --(中略：Bsplineについては表5-8参照)---- ・ ・ ・ def setQ2(X, K):#補間用 Q 設定 N=len(X); Q=array(N*2+1); DX=X[1]-X[0]; Q[K-2]=X[0]-DX for i in range(K - 3, -1,-1): Q[i] = Q[i+1]-DX for i in range(N) : Q[i+K-1] = X[i] DX = X[N-1]-X[N-2] for i in range(K-1): Q[i+N+3] = Q[i+N+2]+DX return Q def setArray(Q, K, Y, N):#配列設定 A=array(N,N+1) for i in range(N): for j in range(N): A[i][j] = Bspline(Q, j+1, K, Q[i + K - 2]) A[i][N]=Y[i] return A def setX(Xstart, NN, DX):# Xの設定(先頭を0とする） XX=array(NN); XX[0] = Xstart- 1 for j in range(1,NN): XX[j] = XX[j-1] + DX return XX def compute(B, XX, Q, K,N):#補間計算 NN=len(XX); YY=array(NN) for j in range(NN): T = 0 for i in range(N): T += B[i][0] * Bspline(Q, i+1, K, XX[j]) YY[j]=T return YY def BSplineInterpo(K, X,Y):#Bスプライン補間 N =len(X); Q=setQ2(X,K); A=setArray(Q, K, Y, N) B=GaussianELMPivot(A) #連立方程式の解法 XX=setX(X[0], 41, 0.1); YY=compute(B,XX, Q, K,N) return [XX,YY] X=[1,2,3,4,5]; Y=[3,5,4,6,4] RL=BSplineInterpo(4, X,Y) root=initTk("Numerical Analysis") canvas=initCanvas(root,440,240);canvas.pack() drawScale(canvas, 30,20, 15, 10, 50, 20, 1, 0, 1,1, 1, 1, 1, 1, 4, 7,"%.f","%.f") for i in range(1, len(RL[0])): drawV(canvas, RL[0][i-1],RL[1][i-1],RL[0][i],RL[1][i],"red", X0=30,Y0=160,Xscale=50,Yscale=20) 表5-10 毎回基底関数を計算しない方法 from graph import * import math def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] def createTable(X, Y, DY1, DYN):#■ 表の生成 N=len(X); H=array(N); D=array(N); Table=array(N) for i in range(N-1) : H[i]= X[i+1] - X[i] for i in range(1,N - 1): D[i]= 2*(X[i+1]-X[i-1]) U = (Y[1]-Y[0])/H[0] for i in range(1, N-1): T= (Y[i+1]-Y[i])/ H[i]; Table[i]=T-U; U=T Table[0] = DY1/6; Table[N-1] = DYN/6; N2=N-2 Table[1] = Table[1] - H[0]*Table[0] Table[N2]= Table[N2]-H[N2]*Table[N-1] for i in range(1,N2): ii=i+1; T = H[i]/D[i]; Table[ii]-=Table[i]*T; D[ii] -= H[i]*T for i in range(N2, 0, -1): Table[i] = (Table[i]-H[i]* Table[i+1]) / D[i] return Table def Fval(X): return (X*X-1)*X #■ F値 def unitInt(XN, X, Y, Table): #■ 単一補間 N=len(X); i=0; j=N-1 while i < j:#補間する区間を見つける K = int((i + j) / 2) # 中央値が if X[K] < XN : i = K+1# 小さければ中央+1をiとする else : j = K # 大きければ中央をjとする if i > 0 : i -= 1 # iを差し引く H = X[i+1]-X[i]; T=(XN-X[i])/H #補間計算 return (T*Y[i+1] + (1-T)*Y[i] + H*H*(Fval(T)*Table[i+1]+Fval(1-T)*Table[i])) def interp(X, Y, Table): #■ 補間 N=len(X); NN=(N-1)*10+1;XX=array(NN); YY=array(NN) XN = X[0] for i in range(NN): XX[i] = XN; YY[i] = unitInt(XN, X, Y, Table) XN += 0.1 return [XX,YY] def BSplineInterpo2(X,Y): Table=createTable(X, Y, (Y[1]-Y[0])/(X[1]-X[0]), (Y[4]-Y[3])/(X[4]-X[3])) return interp(X, Y, Table) X=[1,2,3,4,5]; Y=[3,5,4,6,4] RL=BSplineInterpo2(X,Y) root=initTk("Numerical Analysis") canvas=initCanvas(root,440,240);canvas.pack() drawScale(canvas, 30,20, 15, 10, 50, 20, 1, 0, 1,1, 1, 1, 1, 1, 4, 7,"%.f","%.f") for i in range(1, len(RL[0])): drawV(canvas, RL[0][i-1],RL[1][i-1],RL[0][i],RL[1][i],"red", X0=-20,Y0=160,Xscale=50,Yscale=20)   表5-11 簡便法によるスプライン曲線 from graph import * import math def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] def BSplineSimpleMethod(X,Y): N=len(X); NN=(N-1)*10+1; N1=N-1; N2=N-2 XX=array(NN); YY=array(NN); DY=array(N) for i in range(N): ip=i+1; im=i-1 if i ==0 : DY[i]= (Y[1 ]-Y[0 ])/(X[1 ]-X[0 ]) elif i==N1: DY[i]= (Y[N1]-Y[N2])/(X[N1]-X[N2]) else : DY[i]= (Y[ip]-Y[im])/(X[ip]-X[im]) for k in range(N-1): DX = X[1] - X[0];kp=k+1 DX2 = DX * DX; DX3 = DX * DX2 A1 = 2*(Y[k ]-Y[kp])/DX3+( DY[k]+DY[kp])/DX2 A2 = 3*(Y[kp]-Y[k ])/DX2-(2*DY[k]+DY[kp])/DX A3 = DY[k]; A4 = Y[k]; XN = 0; DX /=10 for j in range(10): ID = k*10+j XX[ID] = X[k]+XN; XN2 = XN*XN; XN3 = XN2*XN YY[ID] = A1*XN3 + A2*XN2 + A3*XN + A4 XN += DX # 変動量を変化させる XX[NN-1]=X[N-1];YY[NN-1]=Y[N-1] return [XX,YY] X=[1,2,3,4,5]; Y=[3,5,4,6,4] RL=BSplineSimpleMethod(X,Y) root=initTk("Numerical Analysis") canvas=initCanvas(root,440,240);canvas.pack() drawScale(canvas, 30,20, 15, 10, 50, 20, 1, 0, 1,1, 1, 1, 1, 1, 4, 7,"%.f","%.f") for i in range(1, len(RL[0])): drawV(canvas, RL[0][i-1],RL[1][i-1],RL[0][i],RL[1][i],"red", X0=-20,Y0=160,Xscale=50,Yscale=20) 表5-12 アルキメデスの螺旋の描画 from tkinter import * import math root=Tk();root.title("Numerical Analysis") canvas=Canvas(root,width=600,height=600);canvas.pack() X1=0; Y1=0; ds=math.pi/32; s=0 for i in range(200): s+=ds; X2=s*math.cos(s); Y2=s*math.sin(s)#パラメトリックな定義 canvas.create_line(300+X1*10, 300-Y1*10, 300+X2*10, 300-Y2*10,width=3,fill='red') X1=X2; Y1=Y2 表5-13 3次スプライン曲線 from tkinter import * import math #■以下，共通的な処理 root=Tk();root.title("Numerical Analysis") canvas=Canvas(root,width=300,height=300);canvas.pack() def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] def TriDiag(A, D):#三重対角行列による連立方程式の解法 N=len(A); X=array(N); CC=array(N) btemp = A[0][1]; X[0] = D[0]/btemp for j in range(1,N): CC[j] = A[j-1][2]/btemp; btemp = A[j][1]-A[j][0]* CC[j] X[j] = (D[j]-A[j][0]*X[j-1])/btemp for j in range(N-2,-1,-1): X[j]-=CC[j+1]*X[j+1] return X #■以下，図形表示用 current=[0,0]#カレント位置 def Move(X,Y): current[0]=X; current[1]=Y def canvasX(X):return 100+X*50 #キャンバス上のX def canvasY(Y):return 200-Y*50 #キャンバス上のY def Draw(canvas, X,Y):#スプライン曲線の線素を描く X1=canvasX(current[0]); Y1=canvasY(current[1]) X2=canvasX(X) ; Y2=canvasY(Y) canvas.create_line(X1,Y1, X2,Y2, width=3,fill='red') current[0]=X; current[1]=Y def drawPoint(canvas,Point):#補間前の線を描く N=len(Point) X1=canvasX(Point[0][0]); Y1=canvasY(Point[0][1]) for i in range(1,N): X2=canvasX(Point[i][0]); Y2=canvasY(Point[i][1]) canvas.create_line(X1,Y1, X2,Y2, width=1,fill='blue') X1=X2; Y1=Y2 #■以下，スプライン曲線計算用 def func(P1, P2, P3, Sigma1, Sigma2): return 3*((P2-P1)*Sigma2/Sigma1+(P3-P2)*Sigma1/Sigma2) def clampedConditio(Sigma,N, Point,P1d,Pnd): #固定条件 A=array(N,3); Pd=array(2,N); Dd=array(2,N); Dtemp=array(N) A[0][1] =2*(Sigma[0]+Sigma[1]); #A[0][2]= Sigma[0]#端の境界条件 A[N-1][0] = Sigma[N-2]; A[N-1][1] = 2*(Sigma[N-3] + Sigma[N-2]) for k in range(1,N-1):#端以外の境界条件 A[k][0] = Sigma[k]; A[k][1] = 2*(Sigma[k - 1]+Sigma[k]) A[k][2] = Sigma[k - 1] for j in range(2):#パラメータ設定 for k in range(1,N-1): Dd[j][k-1]=func(Point[k-1][j], Point[k][j], Point[k+1][j], Sigma[k-1], Sigma[k]) Dd[j][0] = Dd[j][0] - Sigma[1] * P1d[j] Dd[j][N-3]= Dd[j][N-3]-Sigma[N-3] * Pnd[j] for k in range(N): Dtemp[k] = Dd[j][k] Ptemp=TriDiag(A, Dtemp)#連立方程式を解く for k in range(N-1): Pd[j][k + 1] = Ptemp[k] Pd[j][0] = P1d[j]; Pd[j][N-1] = Pnd[j] return [Pd,Dd,A] def relaxedCondition(Sigma, N,Point):#自然条件 A=array(N,3); Pd=array(2,N); Dd=array(2,N); Dtemp=array(N) A[0][1]= 2*Sigma[0];A[0][2] = Sigma[0]#端の境界条件 A[N-1][0]= Sigma[N-2]; A[N-1][1]= 2*Sigma[N - 2]; A[N-1][2]=0 for k in range(1,N-1):#端以外の境界条件 A[k][0] = Sigma[k]; A[k][1] = 2*(Sigma[k]+Sigma[k-1]) A[k][2] = Sigma[k-1] for j in range(2):#パラメータ設定 Dd[j][0 ] = 3*(Point[1 ][j]-Point[0 ][j]) Dd[j][N-1] = 3*(Point[N-1][j]-Point[N-2][j]) for k in range(1,N - 1): Dd[j][k] = func(Point[k - 1][j], Point[k][j], Point[k + 1][j], Sigma[k - 1], Sigma[k]) for k in range(N): Dtemp[k] = Dd[j][k] Ptemp=TriDiag(A, Dtemp)#連立方程式を解く for k in range(N): Pd[j][k] = Ptemp[k] return [Pd,Dd,A] def Spline(canvas, Point, Sigma, P1d, Pnd, ds, Cond):#スプライン曲線 N=len(Point);Ak0=array(2);Ak1=array(2);Ak2=array(2);Ak3=array(2) if Cond : R=clampedConditio(Sigma,N, Point,P1d,Pnd)#固定条件 else: R= relaxedCondition(Sigma, N,Point) #自然条件 Pd=R[0]; Dd=R[1]; A=R[2] Move(Point[0][0], Point[0][1])#最初の点位置に移動 for k in range(N-1):#補間と描画 for j in range(2):#補間用パラメータ設定 Sigma2= Sigma[k]*Sigma[k]; Sigma3=Sigma[k]*Sigma2 Ak0[j]= (2*(Point[k][j]-Point[k+1][j])/Sigma3 +(Pd[j][k]+ Pd[j][k+1])/Sigma2) Ak1[j]= (3*(Point[k+1][j]-Point[k][j])/Sigma2 -(2*Pd[j][k]+Pd[j][k+1])/Sigma[k]) Ak2[j]= Pd[j][k]; Ak3[j]= Point[k][j] sk = 0 while (sk+ds) <= Sigma[k]:#描画 sk += ds PX = ((Ak0[0]*sk+Ak1[0])*sk+Ak2[0])*sk+Ak3[0] PY = ((Ak0[1]*sk+Ak1[1])*sk+Ak2[1])*sk+Ak3[1] Draw(canvas,PX,PY)#線を描く Draw(canvas,Point[N-1][0], Point[N-1][1])#最後の点位置 Point=[[0,0],[0,2],[2,0],[2,2]]; drawPoint(canvas,Point) #Point=[[0,0],[0,2],[2,2],[2,0]]; drawPoint(canvas,Point) Sigma=[1,1,1] #Sigma=[1,2,1] #P1d=[1,1]; Pnd=[1,1] P1d=[1,1]; Pnd=[1,1] Spline(canvas, Point, Sigma, P1d, Pnd, 0.01, True) ■５．５　波形データの補間 表5-14 sinc関数の値 from graph import * #graph.pyについては表4-1参照 import math def sinc(X): if abs(X)<1E-30: return 1 return math.sin(X)/X root=initTk("Numerical Analysis");canvas=initCanvas(root,440,240) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,20, 15, 10, math.pi*10, 20, st, -0.4, DX*20, 0.2, 5, 6, 1, 1, 10, 8,"%.1f","%.1f") X1=st; Y1=sinc(X1); DX=math.pi/10 for i in range(1, 200): X2=X1+DX; Y2=sinc(X2) drawV(canvas, X1,Y1,X2,Y2,"red", X0=math.pi*50+30,Y0=140,Xscale=5,Yscale=100) X1=X2;Y1=Y2 表5-15 sinc関数合成による補間 from graph import * #graph.pyについては表4-1参照 import math def sinc(X): if abs(X)<1E-30: return 1 return math.sin(X)/X def drawCircle(canvas,X1,Y1,X2,Y2,CL,width=3,#円の表示 X0=30,Y0=120,Xscale=100, Yscale=100): XX1=X1*Xscale+X0;YY1=Y0-Y1*Yscale XX2=X2*Xscale+X0;YY2=Y0-Y2*Yscale canvas.create_oval(XX1,YY1,XX2,YY2,fill=CL) def dspWave(canvas, V, TH, CL):#サンプル値別sinc関数の値表示 X1=0; Y1=V*sinc(X1-TH); DX=math.pi/10 for i in range(0, 60): X2=X1+DX;Y2=V*sinc(X2-TH) drawV(canvas, X1,Y1,X2,Y2,CL, X0=30,Y0=120,width=2,Xscale=20,Yscale=100) X1=X2;Y1=Y2 def dspSample(canvas, V, TH):#サンプル値表示 drawV(canvas, TH,0,TH,V,"black", X0=30,Y0=120,width=1,Xscale=20,Yscale=100) drawCircle(canvas, TH-0.2,V-0.05, TH+0.2,V+0.05,"red", X0=30,Y0=120,width=2,Xscale=20,Yscale=100) def synthesis(dt,X):#サンプル値の合成 T=0 for i in range(len(dt)): T+= dt[i]*sinc(X-math.pi*i) return T def dspSynthe(canvas, dt, CL):#合成波の表示 X1=0; Y1=synthesis(dt,X1); DX=math.pi/10 for i in range(0, 60): X2=X1+DX;Y2=synthesis(dt,X1) drawV(canvas, X1,Y1,X2,Y2,CL, X0=30,Y0=120,width=3,Xscale=20,Yscale=100) X1=X2;Y1=Y2 root=initTk("Numerical Analysis");canvas=initCanvas(root,440,240) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,20, 15, 10, math.pi*10, 20, 0, -1, DX*5, 0.2, 5, 5, 1, 1, 12, 10,"%.1f","%.1f") CL=["#777700","#007777","#770077","#007700", "#3333FF","#FF3333","#33FF33","#33FF33"] dt=[0, 0.625,0.625, 0, -0.625,-0.625, 0,0.625] dspSynthe(canvas, dt, "red") for i in range(len(CL)): dspWave(canvas, dt[i], math.pi*i,CL[i]) for i in range(len(CL)): dspSample(canvas, dt[i], math.pi*i) 表5-16 標本化のシミュレーション from graph import * #graph.pyについては表4-1参照 import math def fun(X): return math.sin(X/2)-math.sin(X) def sampling(canvas,W,CL,flag):#データ表示 X1=0; Y1=fun(X1); DX=math.pi/20;R=[] for i in range(0, 120): if (i % 5) ==0: R.append(Y1) X2=X1+DX;Y2=fun(X2) if flag: drawV(canvas, X1,Y1,X2,Y2,CL, X0=30,Y0=120,width=W,Xscale=20,Yscale=50) X1=X2;Y1=Y2 return R def dspDigital(canvas,R): X1=0; DX=math.pi/20*5 for i in range(0, 120,5): X2=X1+DX drawV(canvas, X1, 0, X1, R[int(i/5)],'blue', X0=30,Y0=120,width=2,Xscale=20,Yscale=50) X1=X2 root=initTk("Numerical Analysis");canvas=initCanvas(root,440,240) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,20, 15, 10, math.pi*10, 10, 0, -2, DX*5, 0.2, 5, 5, 1, 5, 12, 20,"%.1f","%.1f") R=sampling(canvas,2,"red",True) dspDigital(canvas,R) 表5-17 階段関数化 ---(前略：未定義部分については表5-16を参照)--- def sampleHold(canvas,R,flag): X1=0; DX=math.pi/20; SH=[] for i in range(0, 120): SH.append(R[int(i/5)]) if flag: CL= "#33FF33" if (i % 5)==0: CL="blue" drawV(canvas, X1, 0, X1, R[int(i/5)],CL, X0=30,Y0=120,width=4,Xscale=20,Yscale=50) X1+=DX return SH root=initTk("Numerical Analysis");canvas=initCanvas(root,440,240) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,20, 15, 10, math.pi*10, 10, 0, -2, DX*5, 0.2, 5, 5, 1, 5, 12, 20,"%.1f","%.1f") R=sampling(canvas,2,"red",True) SH=sampleHold(canvas,R,True) 表5-18 平均化による零次ホールドの補間 ---(前略：未定義部分については表5-16，表5-17を参照)--- def average(SH,i,N): ist=i-N+1; T=0 if ist<0: ist=0; for k in range(ist,i+1):T+=SH[k] return T/ N def dspSH(canvas, S):#多段階平均化による低域フィルタ X1=0; DX=math.pi/20; A=[]; A.append(S[0]); B=[]; B.append(S[0]) C=[]; C.append(S[0]); D=[]; D.append(S[0]) for i in range(1, 120): A.append(average(S,i,4)); B.append(average(A,i,4)) C.append(average(B,i,4)); D.append(average(C,i,4)) X2=X1+DX; drawDT(canvas, X1,S[i-1],X2,S[i],2,'#770077') drawDT(canvas, X1,A[i-1],X2,A[i],2,'green') drawDT(canvas, X1,B[i-1],X2,B[i],2,'#555500') drawDT(canvas, X1,C[i-1],X2,C[i],3,'blue') X1=X2 root=initTk("Numerical Analysis");canvas=initCanvas(root,440,240) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,20, 15, 10, math.pi*10, 10, 0, -2, DX*5, 0.2, 5, 5, 1, 5, 12, 20,"%.1f","%.1f") R=sampling(canvas,2,"red",True) SH=sampleHold(canvas,R,False) dspSH(canvas, SH) 表5-19 二値平均によるローパスフィルタ from graph import * #graph.pyについては表4-1参照 import math def array(N1, N2=0):#配列 if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] 　　・ 　　・ 　　・ ---(中略：その他の関数については表5-18参照)　　 　　・ 　　・ 　　・ def dspSH2(canvas, S):#多段階平均化によるLPF(その2) X1=0; DX=math.pi/20;M=20; A=array(M,2) for i in range(1, 120): A[0][1]=(S[i-1]+S[i])/2 for j in range(1,M): A[j][1]=(A[j-1][0]+A[j-1][1])/2 X2=X1+DX; drawDT(canvas, X1,S[i-1],X2,S[i],2,'#770077') for j in range(0,M-1,3): drawDT(canvas, X1,A[j][0],X2,A[j][1],1,'green') drawDT(canvas, X1,A[M-1][0],X2,A[M-1][1],3,'Blue') X1=X2 for j in range(M):A[j][0]=A[j][1] root=initTk("Numerical Analysis");canvas=initCanvas(root,440,240) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,20, 15, 10, math.pi*10, 10, 0, -2, DX*5, 0.2, 5, 5, 1, 5, 12, 20,"%.1f","%.1f") R=sampling(canvas,2,"red",True) SH=sampleHold(canvas,R,False) dspSH2(canvas, SH) 表5-20 ３次畳み込み関数による補間 from graph import * #graph.pyについては表4-1参照 import math def setDT(N):#データ設定 DT=[] for i in range(N+1): DT.append(math.sin(math.pi*i/3)) return DT def cubicF(t,alfa):#３次畳み込み関数による近似 t1=abs(t); t2=t1*t1; t3=t2*t1 if t1 >= 2: return 0 if t1 >= 1: return alfa*t3 - 5*alfa*t2 + 8*alfa*t1 - 4*alfa; return (alfa + 2)*t3 - (alfa + 3)*t2 + 1 def cubicInterpo(DT,alfa):#三次畳み込み関数による補間 N=len(DT); t=0;V=0;R=[] for i in range(N): for j in range(10): t=0.1*j; V = DT[i ]*cubicF(t ,alfa) if i>=1: V+= DT[i-1]*cubicF(t+1,alfa) if i<=N-3: V+= DT[i+2]*cubicF(t-2,alfa) if i<=N-2: V+= DT[i+1]*cubicF(t-1,alfa) #print(t+i,",",V) R.append([t+i,V]) return R def dspV(canvas, S):#グラフ化 N=len(S) for i in range(1, N): drawV(canvas, S[i-1][0],S[i-1][1],S[i][0],S[i][1], "red", X0=30,Y0=100,width=3,Xscale=100/3,Yscale=50) def dspSin(canvas):#sin関数表示 X1=0; Y1=0 for i in range(1,90): X2=i*0.1;Y2=math.sin(X2*math.pi/3) drawV(canvas, X1,Y1,X2,Y2, "blue", X0=30,Y0=100,width=3,Xscale=100/3,Yscale=50) X1=X2;Y1=Y2 def standardDev(S):#標準偏差計算 N=len(S); T=0 for i in range(N): DY=math.sin(i*0.1*math.pi/3)-S[i][1]; T += DY*DY return math.sqrt(T/N) def dspSample(canvas,DT):#標本値表示 N=len(DT) for i in range(N): drawV(canvas, i,0,i,DT[i], "green", X0=30,Y0=100,width=2,Xscale=100/3,Yscale=50) root=initTk("Numerical Analysis");canvas=initCanvas(root,440,200) st=-10*math.pi;DX=math.pi/10 drawScale(canvas, 30,40, 15, 10, 10, 5, 0, -1.2, 0.3, 0.1, 5, 12, 5, 2, 30, 24,"%.1f","%.1f") DT=setDT(8); alfa=-0.5; S=cubicInterpo(DT,alfa) dspSample(canvas,DT); dspSin(canvas); dspV(canvas, S) D=standardDev(S) canvas.create_text(180,30, text="α = %4.1f 標準偏差=%8.5f" %(alfa, D), font=('Courier',12),fill='black')