■９．１　数値微分 表9-1 刻み幅による誤差の違い def fn(x) :x2=x*x; x3=x2*x; x4=x3*x; return x4+x3+x2+x+1 def dfn(x):x2=x*x; x3=x2*x; return 4*x3+3*x2+2*x+1 def errDiff(): h=0.1 for i in range(21): x=h*i;df0=dfn(x) #微分値 df1=(fn(x )-fn(x-h))/ h #後退差分 df2=(fn(x+h)-fn(x ))/ h #前進差分 df3=(fn(x+h)-fn(x-h))/(2*h)#中心差分 er1=df1-df0; er2=df2-df0; er3=df3-df0 print(("%3.1f "+("%7.3e " * 7)) % (x, df0, df1, df2,df3, er1,er2,er3)) errDiff() 表9-2 刻み幅による誤差の変化 def errStepDiff(): h = 1 for i in range(30): df0 = dfn(1)# 真値 df1 = (fn(1+h)-fn(1))/h #前進差分 e = abs(df1-df0) print("%5.0e, %f" %(h, e)) h /=10 errStepDiff() ■９．２　数値積分「 表9-3 単純数値積分 import math def integV(x):#真値計算用 return -(math.sin(x)+math.cos(x))*math.exp(-x)/2 def fn(x) :return math.sin(x)/math.exp(x)#積分対象関数 def simpleInteg(n,b,a):#単純数値積分 dx=(b-a)/n; x=a;T=0 for i in range(n):x=a+i*dx; T+=fn(x)*dx return T tV= integV(0.35)-integV(0) aV= simpleInteg(10,0.35,0) print("単純積分", aV,"真値", tV,"誤差", aV-tV) 表9-4 台形の公式 ---(前略：　fn(x)とintegVについては表9-3参照)--- def trapezoidal(n,b,a):#台形の公式 dx=(b-a)/n; x2=a;y2=fn(x2);T=0 for i in range(1,n+1): x1=x2;y1=y2;x2=a+i*dx;y2=fn(x2) T+=(y1+y2)*dx/2 return T tV= integV(0.35)-integV(0) aV= trapezoidal(10,0.35,0) print("台形の公式", aV,"真値", tV,"誤差", aV-tV)import math 表9-5 シンプソンの公式 ---(前略：　fn(x)とintegVについては表9-3参照)--- def Simpson(n,b,a):#シンプソンの公式 dx=(b-a)/n; x2=a;y2=fn(x2);T1=0;T2=0 for i in range(1,n,2): T1+=fn(a+dx*i) for i in range(2,n-1,2): T2+=fn(a+dx*i) return (fn(a)+fn(b)+T1*4+T2*2)*dx/3 tV= integV(0.35)-integV(0) aV= Simpson(10,0.35,0) print("シンプソンの公式", aV,"真値", tV,"誤差", aV-tV) 表9-6 ガウス・ルジャンドルの公式 n i ti Ci 2 1 -sqrt(1/3) 1 2 sqrt(1/3) 1 3 1 -sqrt(3/5) 5/9 2 0 8/9 3 sqrt(3/5) 5/9 4 1 -0.8611363116 0.3478548451 2 -0.3399810436 0.6521451549 3 0.3399810436 0.6521451549 4 0.8611363116 0.3478548451 5 1 -0.9061798459 0.2369268851 2 -0.5384693101 0.4786286705 3 0 0.5688888889 4 0.5384693101 0.4786286705 5 -0.9061798459 0.2369268851 表9-7 ガウス積分 import math ti=[-0.9061798459,-0.5384693101, 0, 0.5384693101, 0.9061798459] Ci=[0.2369268851,0.4786286705,0.5688888889, 0.4786286705,0.2369268851] def GausInteg(fn, B_min, B_max): M=len(ti); M1=(B_max-B_min)/2; M2=(B_min+B_max)/2 S = 0 for i in range(M):X=M1*ti[i]+M2; S+=Ci[i]*fn(X) return S*M1 def trueV(B_min, B_max):# 真値（∫(1/(1+x)dx= log(1+x)） return math.log(1+B_max)-math.log(1+B_min) def Fn(X): return 1 / (1+X) #式が変わったらここを変更する。 G=GausInteg(Fn, 0,1); T=trueV(0,1) print("ガウス積分", G,"真値",T,"誤差",G-T) 表9-8 ルジャンドルの多項式を求める def Legen(n):#ルジャンドルの多項式を求める if n==0: return [1] if n==1: return [1,0] P1=Legen(n-1); R1=[]; a1=(2*n-1)/n P2=Legen(n-2); R2=[]; a2=(n-1)/n for i in range(len(P1)): R1.append(P1[i]*a1) R1.append(0) for i in range(len(P2)): R2.append(P2[i]*a2) N1=len(R1); N2=len(R2);NN=N1-N2 for i in range(N1): if i>=NN: R1[i]-=R2[i-NN] return R1 def strExp(P):#多項式の文字列変換 N=len(P);M=N-1; S="" for i in range(N): if P[i]!=0: if i== M : B="" elif i==(M-1): B="x" else: B="x^%d" %(M-i) if i==0 : S+=("%6.3f" %P[i])+B elif P[i]>0: S+=( " +%6.3f"%P[i])+B else:S+=( " %6.3f"%P[i])+B return S print(strExp(Legen(5))) 表9-9 ガウスの積分点を求める 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)] #■以下, DKA法よる複数の解 def cabs(X): return math.sqrt(X.real*X.real+X.imag*X.imag) def initFactor(A, N):# 係数の正規化 for j in range(1,N+1):A[j] /= A[0] return A def initRad(A, N): #初期値用半径の計算 R0 = 0 for j in range(2,N+1): Rn = N * math.pow(abs(A[j]),1.0/j) if Rn > R0: R0 = Rn return R0 def compInit(R0, N): #初期値計算 X=array(N); Pi=math.pi for j in range(N): Th = 2.0*Pi*j/N + 0.5*Pi/N X[j] = R0*complex(math.cos(Th),math.sin(Th)) return X def DKA(A, delta,iterMax):#DKA法 N=len(A)-1; A=initFactor(A,N);X=compInit(initRad(A, N), N) for i in range(iterMax): maxFx = 0 for j in range(N): Fx = 1; Xkj = 1; beforX= X[j] for k in range(N): Fx = Fx* beforX+ A[k + 1] if k != j: Xkj *= (beforX- X[k]) Adj = Fx/Xkj; X[j] = beforX-Adj; E = cabs(Fx) if E > maxFx: maxFx = E if maxFxR[j]: T=R[j-1];R[j-1]=R[j];R[j]=T return R P=Legen(5);print(strExp(P)) Z=DKA(P,0.0001,500);Ti=getTi(Z) 表9-10 ガウスの積分表の作成 def LegenV(n,x):#ルジャンドル多項式の値を求める if n==0: return 1 if n==1: return x P1=LegenV(n-1,x) *x*(2*n-1)/n P2=LegenV(n-2,x)*(n-1)/n return P1-P2 def getCi(n,x):#Ciを求める V=n*LegenV(n-1,x); return 2*(1-x*x)/(V*V) def genTab(N): P=Legen(N);print("\nN=",N, strExp(P)) Z=DKA(P,0.0001,500);Ti=getTi(Z) for i in range(len(Ti)): print("%3d %13.10f %13.10f"%(i,Ti[i],getCi(N,Ti[i]))) for i in range(2,7): genTab(i) 表9-11 多項式の値を求める def exeExp(P,x):#多項式の計算 N=len(P);T=0 for i in range(N):T=T*x+P[i] return T ■９．３　常微分方程式の数値解法 表9-12 前進差分による近似 import math def deffEq1(fn):#前進差分による方法(真値を計算する関数名を引数に指定) dt=0.1; y=1; R=[[0,1,1,0]]#(R[i][2]は真値, R[i][3]は真値との差) for i in range(1,11): t=dt*i;g=(t+1)*y; y+=g*dt;tV=fn(t);R.append([t,y,tV, tV-y]) return R def trueV(t):return math.exp(t*t/2+t)#真値の計算 def printR(L):#結果表示 print("時間,差分近似, 真値 , 誤差") print("-"*35) for i in range(len(L)): print("%4.2f,%8.5f,%8.5f, %8.5e" % (L[i][0],L[i][1],L[i][2],L[i][3])) printR(deffEq1(trueV)) 表9-13 修正オイラー法（修正１）による近似 def deffEq1(fn):#修正オイラー法（修正１） dt=0.1; y=1; R=[[0,1,1,0]]# R:[[時刻,計算値，真値, 誤差]] for i in range(1,11):# k1, k2を求めるときは前時刻で求める t=dt*i;k1=(t-dt+1)*y*dt; k2=(t+1)*(y+k1)*dt y+=(k1+k2)/2;tV=fn(t);R.append([t,y,tV, tV-y]) return R 表9-14 修正オイラー法（修正2）による近似 def deffEq1(fn):#修正オイラー法（修正2） dt=0.1; y=1; R=[[0,1,1,0]]# R:[[時刻,計算値，真値, 誤差]] for i in range(1,11):# k1, k2を求めるときは前時刻で求める t=dt*i;k1=(t-dt+1)*y*dt; k2=(t-dt/2+1)*(y+k1/2)*dt y+=k2;tV=fn(t);R.append([t,y,tV, tV-y]) return R 表9-15 オイラー法と修正オイラー法の誤差比較 時間 オイラー法 修正１ 修正２ ---------------------------------- 0.0 0.00E+00 0.00E+00 0.00E+00 0.1 7.11E-04 2.11E-04 4.61E-04 0.2 2.88E-03 5.40E-04 1.13E-03 0.3 7.17E-03 1.05E-03 2.08E-03 0.4 1.46E-02 1.81E-03 3.46E-03 0.5 2.65E-02 2.97E-03 5.43E-03 0.6 4.51E-02 4.69E-03 8.27E-03 0.7 7.32E-02 7.26E-03 1.23E-02 0.8 1.15E-01 1.11E-02 1.82E-02 0.9 1.78E-01 1.67E-02 2.67E-02 1.0 2.70E-01 2.52E-02 3.90E-02 表9-16 ルンゲ・クッタ法 import math def deffEq1(fn,g,ystart, tstart,tend,N):# 4次のRunge-Kutta法 dt=(tend-tstart)/N; y=ystart R=[[0,1,1,0]]# R:[[時刻,計算値，真値, 誤差]] for i in range(1,N+1):# k1,k2,…を求めるときは前時刻で求める tp=tstart+dt*i;t=tp-dt; tV=fn(tp) #真値 k1 = g(t, y); k2 = g(t+dt/2,y+dt*k1/2) k3 = g(t+dt/2, y+dt*k2/2); k4 = g(t+dt,y+dt*k3) y+=(k1 + 2*k2 + 2*k3 + k4)*dt/6 R.append([tp,y,tV,tV-y]) return R def g(t,y): return (t+1)*y def trueV(t):return math.exp(t*t/2+t)#真値の計算 def printR(L):#結果表示 print("時間,差分近似, 真値 , 誤差") print("-"*37) for i in range(len(L)): print("%4.1f,%8.6f,%8.6f, %8.6e" % (L[i][0],L[i][1],L[i][2],L[i][3])) printR(deffEq1(trueV,g,1, 0,1,10)) 表9-17 ルンゲ・クッタ・ギル法 def deffEq1(fn,g,ystart, tstart,tend,N):# 4次のRunge-Kutta法 dt=(tend-tstart)/N ; y=ystart; sqr2=math.sqrt(2) A1 = -1/2 + 1/sqr2 ; A2 = 1 - 1/sqr2 B1 = -1/sqr2 ; B2 = 1 + 1/sqr2 C1 = 2*(1 - 1/sqr2); C2 = 2*(1 + 1/sqr2) R=[[0,1,1,0]]# R:[[時刻,計算値，真値, 誤差]] for i in range(1,N+1):# k1,k2,…を求めるときは前時刻で求める tp=tstart+dt*i;t=tp-dt; tV=fn(tp) #真値 k1 = g(t, y); k2 = g(t+dt/2,y+dt*k1/2) k3 = g(t+dt/2, y+A1*dt*k1+A2*dt*k2) k4 = g(t+dt,y+B1*dt*k2+B2*dt*k3) y+=(k1 + C1*k2 + C2*k3 + k4)*dt/6 R.append([tp,y,tV,tV-y]) return R 表9-18 連立ルンゲ・クッタ法 import math def array(N1, N2=0):#■配列宣言（表2-1 再掲） if N2!=0: return [[0 for i in range(N2)] for i in range(N1)] return [0 for i in range(N1)] #■連立Runge-Kutta法 def CoalitionRungeKutta(g,Y0, xstart,xend,N):#連立Runge-Kutta法 M=len(Y0); K1=array(M); K2=array(M); K3=array(M) K4=array(M); YY=array(M); Y=array(N+1,M+1) DX=(xend-xstart)/N; X = xstart Y[0][M]=X # Y[i][M]にはXの値を入れる for k in range(M): Y[0][k]=Y0[k]#先頭データには初期値を複写 for i in range(1,N+1): for k in range(M): YY[k] = Y[i-1][k] for k in range(M): K1[k] = g(k, X, YY) for k in range(M): YY[k] = Y[i-1][k] + DX*K1[k]/2 for k in range(M): K2[k] = g(k, X+DX/2, YY) for k in range(M): YY[k] = Y[i-1][k] + DX*K2[k]/2 for k in range(M): K3[k] = g(k, X+DX/2, YY) for k in range(M): YY[k] = Y[i-1][k] + DX*K3[k] for k in range(M): K4[k] = g(k,X+DX, YY) for k in range(M):#Yの計算 Y[i][k] = (Y[i-1][k]+ (K1[k]+ 2*K2[k]+2*K3[k]+K4[k])*DX/6) X+=DX;Y[i][M]=X return Y def g(ID, X, Y):#IDでY1/Y2を判別 if ID : return X #IDが0でないとき return X+Y[1] #IDが0のとき def printR(L):#結果表示 print("時間, Y1 , Y2 "); print("-"*24) for i in range(len(L)): print("%4.1f,%9.6f,%9.6f" % (L[i][2],L[i][0],L[i][1])) Y0=[0,0]#初期値をY1=0, Y2=0とする printR(CoalitionRungeKutta(g,Y0,0,1,11)) 表9-19 連立ルンゲ・クッタ法の変更 ---(前略：他の部分については表9-14参照)----- def g(ID, X, Y): if ID : return -Y[0] #IDが0でないとき return Y[1] #ID=0のとき def trueV(t):return math.exp(t*t/2+t)#真値の計算 def printR(L):#結果表示 print("時間, Y1 , Y2 , 真値 , 誤差") print("-"*47) for i in range(len(L)): XT=4*math.sin(L[i][2]) print("%4.1f,%9.6f,%9.6f,%9.6f, %+9.4e" % (L[i][2],L[i][0],L[i][1],XT,L[i][0]-XT )) #print(L[i][2],",",L[i][0],",",L[i][1],",",XT,",",L[i][0]-XT) Y0=[0,4] printR(CoalitionRungeKutta(g,Y0,0,2*math.pi,60)) 表9-20 ルンゲ・クッタ法を使わない方法（オイラー法） import math def euler(v0,x0): t=0; dt=0.1; v=v0; x=x0 for i in range(51): t=dt*i; tx=4*math.sin(t) print("%3.1f %8.4f %8.4f %8.4f %8.4e" % (t,v,x,tx,tx-x)) #表示用 #print(t,",",v,",",x,",",tx,",",tx-x)#他ツール用 vb=v;xb=x; v=v-xb*dt; x=x+vb*dt 表9-22 テーラ展開による方法 import math def euler(v0,x0): t=0; dt=0.1; v=v0; x=x0 for i in range(51): t=dt*i; tx=4*math.sin(t) print("%3.1f%8.4f%8.4f%8.4f% +12.4e" % (t,v,x,tx,tx-x)) #表示用 #print(t,",",v,",",x,",",tx,",",tx-x)#他ツール用 vb=v; xb=x x=xb*math.cos(dt)+vb*math.sin(dt) v=vb*math.cos(dt)-xb*math.sin(dt) euler(4,0)