■１２．２　山登り法／山下り法

表12-1  最も単純な山下り法
def fn(x):return 3*x*x-10*x+5 #目的関数が異なる場合はここを変更
def dawnHill(x0, dx, Nloop):
    x=x0; i=1; Fx=fn(x); FxP=fn(x+dx); FxM=fn(x-dx);R=[]
    while i<Nloop and (Fx>FxP or Fx>FxM):
        R.append([x, Fx, FxP, FxM]) #　チェック用
        if FxP <= FxM: x +=dx #目的関数の値が減少する方向に移動
        else: x -=dx
        i+=1;Fx=fn(x); FxP=fn(x+dx); FxM=fn(x-dx)
    R.append([x, Fx, FxP, FxM]) #　チェック用
    return R
R = dawnHill(0, 0.001, 5000)
print(R[len(R)-1][0])
print("i,x,Fx,FxP,FxM")
for i in range(len(R)):
    print(i+1,",", R[i][0],",", R[i][1],",", R[i][2],",", R[i][3])

表12-2  改良された山下り法
def fn(x):return 3*x*x-10*x+5 
def dawnHill2(x0, dxin, epsi, Nloop) :
    x = x0; dx = dxin; i=1;R=[]
    while i<Nloop and dx>epsi:
        Fx=fn(x); FxP=fn(x+dx); FxM=fn(x-dx)
        while i < Nloop and (Fx>FxP or Fx>FxM):
            R.append([x,dx,Fx, FxP, FxM]) #　チェック用
            if FxP <= FxM: x+=dx #目的関数の値が減少する方向に移動
            else: x -= dx
            i+=1; Fx=fn(x); FxP=fn(x+dx); FxM=fn(x-dx)
        dx /= 2
    R.append([x,dx, Fx, FxP, FxM]) #　チェック用
    return R
R = dawnHill2(0,1, 0.0005, 5000)
print("i,x,dx,Fx,FxP,FxM")
for i in range(len(R)):
    print(i+1,",", R[i][0],",", R[i][1],",", R[i][2],",",
          R[i][3],",", R[i][3])

表12-3  N次元の山下り法
from tkinter import *
import random, math
root=Tk()
canvas=Canvas(root, width=540,height=540, bg="white");canvas.pack()
def fn(x) :#目的関数が異なる場合はここを変更
    return 3*x[0]*x[0]-5*x[0]+5*x[1]*x[1]-20*x[1]-3*x[0]*x[1]+10
#■以下，カラーマップ上に収束の様子を表示するための関数
def fnMinMax():# 関数値の最小，最大を求める
    mn=9999999;mx=-99999999
    for i in range(401):
        x=[i/100,0]
        for j in range(401):
            x[1]=j/100; f=fn(x)
            if f<mn: mn=f
            if f>mx: mx=f
    return [mn,mx]
def codeCH(CD):# RGB値のチェック
    if CD<0: return 0
    if CD>255: return 255
    return int(CD)
def colorCode(H,minMax):# 関数値を色コードで表現
    V=(H-minMax[0])/(minMax[1]-minMax[0])
    if   V<0.25: R=0; G= int(127*V*4); B=int(127*(0.25-V)*4)
    elif V<0.5 : R=0; G=(127*(V-0.25)*4+127); B=0
    elif V<0.75: R=int((V-0.5)*4*255); G=255; B=0
    else       : R=255; G=int((1-V)*4*255)  ; B=0
    return "#%02X%02X%02X" %(codeCH(R),codeCH(G),codeCH(B))
def colorMap():# カラーマップ表示
    minMax=fnMinMax()
    for i in range(401):#カラーマップ
        x=[i/100,0]
        for j in range(401):
            x[1]=j/100; f=fn(x);Color=colorCode(f,minMax)
            canvas.create_rectangle(i+10,410-j,i+10,410-j,
                                    fill=Color,outline="")
    for i in range(0, 401,100):#横線
        canvas.create_line(10,i+10,410,i+10,fill="black")
    for i in range(0, 401,100):#縦線
        canvas.create_line(i+10,10,i+10,410,fill="black")
def drawXY(XY):# N次元山下り法の軌跡を描く
    X1=XY[0][1][0]*100+10;Y1=410-XY[1][0][1]# dxの変化を描画
    for i in range(1,len(XY)):
        X2=XY[i][1][0]*100+10;Y2=410-XY[i][1][1]*100
        canvas.create_line(X1,Y1,X2,Y2,fill="indigo",width=2)
        X1=X2;Y1=Y2
    X1=XY[0][0][0]*100+10;Y1=410-XY[0][0][1]# xの変化を描画
    for i in range(1,len(XY)):
        X2=XY[i][0][0]*100+10;Y2=410-XY[i][0][1]*100
        canvas.create_line(X1,Y1,X2,Y2,fill="yellow",width=2)
        canvas.create_text(X1+8,Y1+8,text="%3d" %i,font=("Times",8))     
        X1=X2;Y1=Y2
colorMap()#カラーマップ表示
#■以下 N次元山下り法
def dFx(x, i, ddx):# ΔxまたはΔy変化したときの値(i番目の変数)
    XTemp = x[i];x[i]+=ddx; dF=fn(x);x[i]= XTemp
    return dF
def copyX(x):#Pythonのリスト内変数の値変化を避ける処理
    X=[]
    for i in range(len(x)):X.append(x[i])
    return X
def checkEpsi(N, dx, epsi):# dxが微小範囲を超えるかどうか
    for i in range(N):
        if dx[i] > epsi : return True
    return False
def NDdawnHill(x0, dx0, epsi, Nloop):#N次元山下り法
    N=len(x0); x=[]; dx=[];R=[]
    for i in range(N):#初期値の複写
        x.append(x0[i]);dx.append(dx0[i])
    i = 1;Fx=fn(x)
    R.append([copyX(x), copyX(dx),Fx]); i+=1 
    while i< Nloop and checkEpsi(N, dx, epsi):
        Fx = fn(x)
        for j in range(N):# J番目変数を変化させる
            FxP=dFx(x,j,dx[j]); FxM=dFx(x,j,-dx[j])
            while i < Nloop and (Fx>FxP or Fx>FxM):
                if FxP<=FxM: x[j]+=dx[j]; Fx=FxP
                else       : x[j]-=dx[j]; Fx=FxM
                R.append([copyX(x),copyX(dx),Fx]); i += 1 
                FxP=dFx(x,j,dx[j]); FxM=dFx(x,j,-dx[j])
            dx[j]/= 2
    R.append([copyX(x), copyX(dx),Fx])
    return R
R = NDdawnHill([0,0],[1,1], 0.0005, 5000); drawXY(R)
print("  i,      x  ,      y  ,     dx  ,     dy  ,    Fx")
fmt= "%3d"+(", %8.4f"*5) 
for i in range(len(R)):
    print(fmt %((i+1),R[i][0][0],R[i][0][1],R[i][1][0],
                R[i][1][1],R[i][2]))

表12-4  目的関数の変更
def fnP(x, ip):
    xx=[]
    for i in range(len(x)):xx.append(x[i][ip])
    return fn(xx)
def setFuncValue(x):
    R=[]
    for ip in range(len(x[0])):R.append(fnP(x,ip))
    return R

表12-5  初期単体頂点の設定
def array(N1, N2=0):#配列
    if N2 == 0 :return [0 for i in range(N1)]
    else: return [[0 for i in range(N2)] for j in range(N1)]
def setInitVertex(N,NP, x, h):#初期単体頂点の設定
    x=array(N,NP)
    x[0][0] = 1; x[1][0] = 1 #初期値は0でない値を設定すること
    if N == 1 : x[0][1] = x[0][0]*h
    else:
        for ip in range(1,N+1):
            for i in range(N):
                if   i<(ip-1): x[i][ip] = x[i][0]*h
                elif i>(ip-1): x[i][ip] = x[i][0]/h
                else         : x[i][ip] = x[i][0]
    return x

表12-6  収束判定値の計算
def sqrSigma(NP, fvalue):#平均との差の平方二乗和
   s = 0             
   for ip in range(NP):s += fvalue[ip]
   av = s / NP; s = 0 #平方二乗和
   for ip in range(NP):A=fvalue[ip]-av; s+= A*A
   return s

表12-7  目的関数の値による並替え
def sortVertex(NP, N, fvalue, x): #頂点の並替え
    for ip in range(NP):
        minV = fvalue[ip]; minIP = ip
        for j in range(ip+1,NP):#着目点より後尾の最小値
            if fvalue[j]<minV: minV=fvalue[j]; minIP=j
        if minIP != ip : #着目点と最小値の入替え
            T=fvalue[minIP];
fvalue[minIP]=fvalue[ip];fvalue[ip]=T
            for i in range(N):
              T=x[i][minIP];
x[i][minIP]=x[i][ip];x[i][ip]=T
    return [fvalue, x]

表12-8  変化位置の計算と移動量
def setCenterOfGravity(N, x) :#最不適点以外の重心座標を設定
    centerOfGravity=[]
    for i in range(N):# 最不適点座標はx(1,N+1)～x(N,N+1))なので
        s = 0         # x(1,1)～x(1,N)から重心座標を求める。
        for ip in range(N):s += x[i][ip]
        centerOfGravity.append(s/N)
    return centerOfGravity
def setMirrorCand1(NP,N,xMirror,alpha,centerOfGravity,x):#鏡像点設定
    for i in range(N):
        #print("i=",i,"NP=",NP,"N=", N)
        xMirror[i][0]=(1+alpha)*centerOfGravity[i]-alpha*x[i][NP-1]
    return xMirror
def setMirrorCand2(N, xMirror, beta, centerOfGravity):#拡大点設定
    for i in range(N):
        xMirror[i][1]=beta*xMirror[i][0]+(1-beta)*centerOfGravity[i]
    return xMirror
def setMirrorCand3(NP,N,xMirror,gamma,centerOfGravity,x):#収縮点設定
    for i in range(N):
        xMirror[i][1]=gamma*x[i][NP-1]+(1-gamma)*centerOfGravity[i]
    return xMirror
def moveArray(N, x, ip, xm, ID): #点の移動
    for i in range(N):x[i][ip-1] = xm[i][ID]
    return x
def reduction(NP, N, x):#縮小処理
    for ip in range(1,NP):
        for i in range(N):
            x[i][ip] = (x[i][0]+x[i][ip])/2
    return x

表12-9  その他の処理
def copyXP(x):
    DT=[]
    for i in range(len(x)):
        for j in range(len(x[0])):DT.append(x[i][j])
    return DT
def forCheck(k,NP, N, x, fvalue):
    print("k= %5d" % k)
    for ip  in range(NP):
        S="ip=%2d" % ip
        for i in range(N): S+="%8.5f"%x[i][ip]
        S+=" %8.5f" % fvalue[ip]
        print(S)
def drawXP(XY):# 滑降シンプレックス法の軌跡を描く
    #X1=XY[0][0]*100+10;Y1=410-XY[0][3]# xの変化を描画
    for i in range(0,len(XY)):
        X1=XY[i][2]*100+10;Y1=410-XY[i][5]*100
        for k in range(3):
            X2=XY[i][k]*100+10;Y2=410-XY[i][k+3]*100
            canvas.create_line(X1,Y1,X2,Y2,fill="yellow",width=2)
            X1=X2;Y1=Y2
    X1=XY[0][0]*100+10;Y1=410-XY[0][3]*100# xの変化を描画
    for i in range(1,len(XY)):
        X2=XY[i][0]*100+10;Y2=410-XY[i][3]*100# xの変化を描画
        canvas.create_line(X1,Y1,X2,Y2,fill="red",width=2)
        X1=X2;Y1=Y2
    root.update()

表12-10  全体制御
def downHillSimplex():
    # x:座標データ, fvalue[]: 目的関数の値
    # centerOfGravity[] :最不適点以外の重心座標
    # xMirror[] :鏡像点の候補
    N = 2; x=array(N,N+1); fvalue=array(N + 1);xMirror=array(N,2)
    eps=0.1; errDT=0.00001;alpha=1; beta=2; gamma=0.5; h=0.5
    Nloop = 100;RET=[]; NP=N+1; sqrErr=(errDT*errDT)*NP
    x = setInitVertex(N, NP,x,h);fvalue=setFuncValue(x)
    RET.append(copyXP(x))
    for k in range(Nloop):
        fx = sortVertex(NP, N, fvalue, x)#単体頂点値の順序並び替え
        fvalue=fx[0]; x=fx[1]
        RET.append(copyXP(x)); forCheck(k,NP, N, x, fvalue)
        if sqrSigma(NP, fvalue) < sqrErr : break #収束判定
        centerOfGravity=setCenterOfGravity(N,x)#最不適点以外の重心
        xMirror=setMirrorCand1(NP,N,xMirror,alpha,
                               centerOfGravity, x)#鏡像点(添字0)
        fwm1 = fnP(xMirror, 0)
        if fwm1 < fvalue[N-1]:#鏡像点が最不適点の次の点より小さい
            if fwm1 < fvalue[0]: #鏡像点が最適点より小さい
                xMirror=setMirrorCand2(N,xMirror,beta,
                                       centerOfGravity)#拡大点(1)
                fwm2 = fnP(xMirror, 1)
                if fwm2 < fwm1 :
                    x= moveArray(N, x, NP, xMirror, 1)
                    fvalue[NP-1] = fwm2
                else:
                    x=moveArray(N, x, NP, xMirror, 0)
                    fvalue[NP-1] = fwm1
            else:#鏡像点が最不適点の次の点より小さく最適点より大
               x= moveArray(N, x, NP, xMirror, 0)
               fvalue[NP-1] = fwm1
        else:
            if fwm1 < fvalue[NP-1]:#鏡像点が最不適点より小さい
                x= moveArray(N, x, NP, xMirror, 0)
                fvalue[NP-1] = fwm1
            xMirror=setMirrorCand3(NP, N, xMirror, gamma,
                                   centerOfGravity, x)#収縮点計算
            fwm3 = fnP(xMirror, 1)
            if fwm3 < fvalue[NP-1]:#収縮点が最不適点より小さい
                x=moveArray(N, x, NP, xMirror, 1)
                fvalue[NP-1] = fwm3
            else:#収縮点が最不適点より大きい時
                x=reduction(NP, N, x)#縮小処理で次の鏡像変化に備える
    return RET
R=downHillSimplex()
drawXP(R)

■１２．３　勾配法

表12-11  1次元の勾配法
def fn(x):return 3*x*x-10*x+5 #'目的関数が異なる場合はここを変更
def dfn(x): return 6*x-10 #勾配法ではこの関数を追加
def forCheck1DG(i, x, Fx, FxP, dx):#繰り返しの様子を観察するための関数
    fm=("%5d"+", %8.5f"*4)
    print(("%5d"+", %8.5f"*4)%(i,x,dx, Fx,FxP))
def gradient1D(x0, dxin, epsi, Nloop):
    x = x0; dx = dxin; i = 1
    while i < Nloop and dx > epsi:
        Fx = fn(x); dFx = dfn(x);
        ddx =-dFx*dx/abs(dFx); FxP = fn(x+ddx)
        while i < Nloop and Fx > FxP:
            forCheck1DG(i, x, Fx, FxP, dx)
            i+=1; x +=ddx
            Fx=fn(x); dFx=dfn(x);
            ddx=-dFx*dx/abs(dFx); FxP=fn(x+ddx)
        dx/=2
    forCheck1DG(i, x, Fx, FxP, dx)
    return x
gradient1D(0,1,0.0005,5000)

表12-12  微分値を直接求めることができない場合

def dfnX1(x): EDX=0.001; return fn(x+EDX)-fn(x-EDX)

def dfnX2(x):
    n=int(x); eps=x-n
    df1=(F[n+1]-F[n-1])/2; df2=(F[n+2]-F[n])/2
    return eps*(df2-df1)+df1
 
表12-13  1次元勾配法の拡張による方法
def fnP(x, ip):
    xx=[]
    for i in range(len(x)):xx.append(x[i][ip])
    return fn(xx)
def deltaFn(x, i):
    if i ==0: return 6*x[0]-5-3*x[1] #∂f/∂x
    return 10*x[1]-20-3*x[0]         #∂f/∂y
def dFx(x, i, ddx):
    XTemp = x[i]; x[i]+=ddx;dF = fn(x)
    x[i] = XTemp               #x[i]を元に戻す
    return dF
def forNDG(i, x, Fx, pDX):#繰り返しを表示
    fm="%3d"+", %8.5f"*4
    print(fm % (i,x[0],x[1],Fx,pDX))
def copyX(x):
    R=[]
    for i in range(len(x)):R.append(x[i])
    return R
def NDgradient(x0, dxin, eps, Nloop):
    N=len(x0);x=copyX(x0);i=1; pDX = 1;RET=[]
    forNDG(i,x,fn(x), pDX); i+=1;RET.append(copyX(x))
    for k in range(Nloop):
        for j in range(N):
            Fx=fn(x);flag = True; dx = dxin
            while dx>eps: # flag:移動しなかったことを示す
                dFdX=deltaFn(x,j);pDX =-dx*dFdX/abs(dFdX)
                FxP =dFx(x,j,pDX)
                while Fx>FxP:
                    x[j]+=pDX; Fx=FxP;flag = False
                dx /= 2
            RET.append(copyX(x))
            forNDG(i,x, Fx, pDX);i+=1
        if flag : break #'移動しなかったら終わらせる。
    RET.append(copyX(x))
    forNDG(i, x, Fx, pDX)
    return RET
def drawNDG(XY):
    X1=XY[0][0]*100+10;Y1=410-XY[0][1]*100
    for i in range(1,len(XY)):
        X2=XY[i][0]*100+10;Y2=410-XY[i][1]*100
        canvas.create_line(X1,Y1,X2,Y2,fill="red",width=2)
        canvas.create_text(X1+8,Y1+8,text="%3d" %i,font=("Times",8))     
        X1=X2;Y1=Y2
R=NDgradient([0,0],1,0.00000005,5000)
drawNDG(R)

表12-14  最急降下法
（省略している関数については，表12-3．表12-13を参照）
def array(N1, N2=0):#配列
    if N2 == 0 :return [0 for i in range(N1)]
    else: return [[0 for i in range(N2)] for j in range(N1)]
def NDsteepest(x0, dxin, eps, Nloop):#最急降下法
    N=len(x0);x=copyX(x0);i=1; ddx = dxin;RET=[]
    dx=array(N);Xn=array(N); Fx = fn(x); flag = True
    forNDG(i,x,Fx, ddx); i+=1;RET.append(copyX(x))
    while ddx > eps:
        for j in range(N):#　移動点の計算
            dx[j]=-deltaFn(x,j)*ddx;Xn[j]=x[j]+dx[j]
        FxP = fn(Xn)
        while Fx>FxP:# 現的関数値が移動点より大きい場合
            for j in range(N):x[j]= Xn[j] #新しい点に移動
            Fx = FxP
            for j in range(N):# 移動点の計算
                dx[j]=-deltaFn(x,j)*ddx; Xn[j]=x[j]+dx[j]
            RET.append(copyX(x))
            forNDG(i,x, Fx, ddx);i+=1
            FxP = fn(Xn)
        ddx /= 2
    RET.append(copyX(x))
    forNDG(i,x, Fx, ddx);i+=1
    return RET
R=NDsteepest([0,0],1,0.00000005,5000)
drawNDG(R)

■偏微分値を直接求めることができない場合

表12-12と同様，降下方向を微小幅の差で代替することができますが，最急降下法では，XY方向の傾きが問題となりますので，以下のように中心差分値を求める必要があります。

def deltaFnX(x, j):
    Edx=0.001; Tm = x[j];x[j]=Tm+Edx;
    dfP=fn(x);x[j]= Tm-Edx;dfM = fn(x)
    RET = (dfP-dfM)/(2*Edx);x[j]=Tm
    return RET

■１２．４　ニュートン・ラプソン法

表12-15  1次元のニュートン・ラプソン法
def fn(x):#収束判定には使わない
    return x*x*x - 2*x*x - 7*x + 30
def dfn(x):return 3*x*x - 4*x - 7 #微分
def ddfn(x): return 6*x - 4       #2階微分
def NewtonRaphson(Xint, eps):
    x2 = Xint; E = 100
    while E > eps:
        print("%8.4f, %8.4f" %(x2, fn(x2)))
        x1=x2;x2=x1-dfn(x1)/ddfn(x1);E=abs(x2-x1)
NewtonRaphson(4,0.0000000001)
 

表12-16  1次元の拡張によるニュートン・ラプソン法

from tkinter import *
import random, math
root=Tk()
canvas=Canvas(root, width=540,height=540, bg="white");canvas.pack()
def fn(xy) :#収束判定には使わない（描画用）
    x=xy[0];y=xy[1];x2=x*x;x3=x2*x; y2=y*y; y3=y2*y
    return 5*x3+5*y3-7*x2-10*y2-5*x-20*y-3*x*y
---(中略):表示用関数については表12-3，表12-13を参照------
def dfnx(x, j): #1階微分値の計算
    if j==0: return 15*x[0]*x[0]-14*x[0]-5-3*x[1]#∂f∂x
    if j==1: return 15*x[1]*x[1]-20*x[1]-20-3*x[0]#∂f/∂y
    return 0
def ddfn(x, i, j):#2階微分値の計算
    h = i * 100 + j #以下の判定を単純化するためこの計算を行う。
    if h==  0: return 30*x[0]- 14#∂^2 f/∂x^2
    if h==  1: return -3         #∂^2 f/∂x∂y
    if h==100: return -3         #∂^2 f/∂y∂x
    if h==101: return 30*x[1]-20 #∂^2 f/∂y^2
    return 0
def forNewR(ipt, x):#表示
    print("%2d, %8.5f, %8.5f"%(ipt,x[0],x[1]))
def copyX(x):
    R=[]
    for i in range(len(x)):R.append(x[i])
    return R
def array(N1, N2=0):#配列
    if N2 == 0 :return [0 for i in range(N1)]
    else: return [[0 for i in range(N2)] for j in range(N1)]
def ND1NewtonRaphson(Xint, eps, Nloop):
    N=len(Xint);X2=copyX(Xint);x=array(N);DE=eps*eps;RET=[]
    ipt=1; forNewR(ipt,X2);RET.append(copyX(X2))
    for k in range(Nloop):
        x0=copyX(X2)#現在の座標位置を保存
        for i in range(N):#各次元別に収れん計算
            E = 100
            while E > eps:
                x[i] =X2[i]
                X2[i]=x[i]-dfnx(x,i)/ddfn(x,i,i)
                E = abs(X2[i]-x[i])
            ipt+=1;forNewR(ipt,X2);RET.append(copyX(X2))
        s = 0#最初の座標位置との差を計算
        for i in range(N): ds=X2[i]-x0[i]; s+=ds*ds
        if s < DE :break #動いていなければ繰返し終わり
    RET.append(copyX(X2)); return RET

表12-17  N次元ニュートン・ラプソン法
なお，関数dfnx，ddfn，手続きforCheckは表12-16と同じ，ガウスの消去法の関数GaussianELMPivotは2.4節の表2-10に示した関数です。
def NDNewtonRaphson(Xint,eps, Nloop):
    N=len(Xint);X2=copyX(Xint);x=array(N);DE=eps*eps;RET=[]
    ipt=1;forNewR(ipt,X2);RET.append(copyX(X2))
    NP=N+1; Hessen=array(N, NP)
    for K in range(Nloop):
        x0=copyX(X2)#現在の座標位置を保存
        for i  in range(N ):#Hessen行列設定
            for j in range(N):
                Hessen[i][j]= ddfn(X2,i,j)
        for i in range(N):#右辺の微分値設定
            Hessen[i][N]=-dfnx(X2, i)
        B=GaussianELMPivot(Hessen)#連立方程式を解く
        for i in range(N):
            X2[i] = x0[i]+B[i][0]
        ipt+=1;forNewR(ipt,X2);RET.append(copyX(X2))
        s = 0#動いた距離の2乗を計算
        for i in range(N): ds=X2[i]-x0[i]; s+=ds*ds
        if s<DE : break #動いた距離が短ければ繰返し終了
    RET.append(copyX(X2)); return RET
R=NDNewtonRaphson([4,4],0.00000005,5000)
#R=ND1NewtonRaphson([4,4],0.00000005,5000)
drawNDG(R)

表12-18  微分の替りに差分を用いる方法
以下の関数と手続き以外は，表12-16，表12-17と同じです。
def dfnx(x, j):
    esp=0.000001; h=esp/2; xj = x[j]
    x[j] = xj+h; fnPh=fn(x);x[j]=xj-h;fnMh=fn(x)
    x[j] = xj
    return (fnPh-fnMh)/esp      #esp の値は 2*h
def ddfn(x, i, j):
    esp=0.000001; h=esp/2; h2 = h*h
    if i == j :
        xi=x[i]; fn0=fn(x); x[i]=xi+h;
        fnPh=fn(x); x[i]=xi-h;fnMh=fn(x);x[i]=xi
        return (fnPh - 2*fn0 + fnMh) / h2
    else:
        xi = x[i]; xj=x[j] # x, y を保存
        x[i]=xi+h; x[j]=xj+h; fnPP=fn(x)# f(x + h, y + h)
        x[i]=xi+h; x[j]=xj-h; fnPM=fn(x)# f(x + h, y - h)
        x[i]=xi-h; x[j]=xj+h; fnMP=fn(x)# f(x - h, y + h)
        x[i]=xi-h; x[j]=xj-h; fnMM=fn(x)# f(x - h, y - h)
        x[i]=xi  ; x[j]=xj #x,y を元に戻す
        return (fnPP - fnPM - fnMP + fnMM) / (4*h2)

■１２．５　各種の注意点


表12-19  意味のない範囲を持つとき
def fn(x):
    # x[0]:同軸ケーブルの内部導体内径,x[1]:同軸ケーブルの外部導体外径
    # 形状的には x[1] > x[0] でなければならない。
    if x[1] <= x[0] or x[0]<0.00001:return 1E32 #意味のない範囲のとき
    return (x[1]/x[0]+1) / math.log(x[1]/x[0])

表12-20  線形補間
def linearInt(table, xc, k):
    numDT=len(table);eps=0.000000001;ii=0
    for i in range(numDT):
        if abs(xc-table[i][0])>eps and xc<table[i][0]:break
        ii=i
    if ii>=(numDT-1):ii=numDT-2
    x0=table[ii+1][0]-table[ii][0]; y0=table[ii+1][k]-table[ii][k]
    d=xc-table[ii][0];return (y0/x0)*d+table[ii][k]
table=[[1,10,90],[2,50,75],[3,70,65],[4,80,60],[5,85,58]]
for i in range(0,51):
    x=i/10+1;y1=linearInt(table,x,1);y2=linearInt(table,x,2)
    print(x,y1,y2)

表12-21  スプライン補間
def cmpSpline(startB, startV, endB, endV, X, table, ID):
    N=len(table); XID = int(X); IID = -1
    if XID>table[N-1][0]: IID=N-1
    else:
        for i in range(N):
            if table[i][0]==XID: IID = i; break
        if IID < 0 : return 0
    F2 = table[IID][ID]
    if IID == 0:
        F1 = startV; F3 = table[1][ID];F4=table[2][ID]
        if not(startB): F1=2*F2-F3
    else:
        F1 = table[IID-1][ID];F3=endV;F4=endV
        if IID ==(N-1):
            if not(endB):DF=F1-F2;F3=F2-DF;F4=F3-DF
        elif IID ==(N-2):
            F3=table[IID+1][ID]
            if not(endB):DF=F2-F3;F4=F3-DF
        else:
            F3=table[IID+1][ID];F4=table[IID+2][ID]
    DF1=(F3-F1)/2; DF2=(F4-F2)/2
    A11=2*(F2-F3)+(DF1+DF2); A12=3*(F3-F2)-(2*DF1+DF2)
    A13=DF1; A14=F2; H=X-XID
    return (A11*H*H*H + A12*H*H + A13*H + A14)
table=[[1,10,90],[2,50,75],[3,70,65],[4,80,60],[5,85,58]]
for i in range(50):
    x=i/10+1;y1=cmpSpline(False,0,False,0,x,table,1)
    y2=cmpSpline(False,0,False,0,x,table,2)
    print(x,",",y1,",",y2)