■４．１　非線形方程式の種類

表4-1  グラフ描画のための関数（graph.pyとして保存）
from tkinter import *
import math
def initTk(Title):#Tkモジュール起動
    root=Tk();root.title(Title); return root
def initCanvas(root,W,H):#キャンパスの定義
    canvas=Canvas(root,width=W,height=H)
    canvas.pack(); return canvas
#■目盛り表示
def drawScale(canvas,X0,Y0,X1,Y1,Xsc,Ysc,Xmin,Ymin,dx,dy,xb,yb,
xl,yl,nx,ny,format="%.1f",yformat="%.1f"):
    #  canvas : canvasオブジェクト
    #  X0, Y0 : 描画グラフの左上隅座標 
    #　X1 　　: 左側からのY目盛ラベル表示相対位置
    #　Y1 　　: 最下端からのX目盛ラベル表示相対位置
    #　Xsc　　: X目盛の描画間隔
    #　Ysc　　: Y目盛の描画間隔
    #  Xmin   : Xの最小値
    #  Ymin   : Yの最小値
    #  dx, dy : X, Y の刻み幅
    #  xb, yb : 太線にする間隔(本数)
    #  xl, yl : ラベル表示する間隔(本数)
    #  nx, ny : 目盛りの本数
    #  xformat: Xラベルの表示書式
    #  yformat: Yラベルの表示書式
    CL3="#000000";
    V=Ymin+ny*dy; Xmax=X0+Xsc*nx; Ymax=Y0+Ysc*ny
    for i in range(ny+1):#縦方向目盛り表示
        Y=i*Ysc+Y0; CL="#777777"
        if i % yb: CL="#CCCCCC"
        canvas.create_line(X0,Y,Xmax,Y,width=1,fill=CL)  
        if i % yl==0:
            canvas.create_text(X0-X1,Y,text=(yformat % (V+1E-7)),
                               font=('Times',8),fill=CL3)
        V-=dy
    V=Xmin
    for i in range(nx+1):#横方向目盛り表示
        X=i*Xsc+X0; CL="#777777"
        if i % xb: CL="#CCCCCC"
        canvas.create_line(X,Y0,X,Ymax,width=1,fill=CL)
        if i % xl==0:
            canvas.create_text(X,Ymax+Y1,text=(xformat % V),
                               font=('Times',8),fill=CL3)
        V+=dx
    canvas.create_rectangle(X0,Y0,Xmax,Ymax,
                            width=1,fill="",outline="black")
#■値の描画
def drawV(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_line(XX1,YY1,XX2,YY2,width=width,fill=CL)

■４．２　実根の求め方

表4-2  f (x) = x2−4x + 1のグラフ化
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 20, 20, 0, -4,
          0.5, 1, 2, 2, 2,1, 8, 6,"%d","%.1f")
def fn(X):return X*X-4*X+1
X2=0; Y2=fn(X2); DX=0.1
for i in range(40):
    X1=X2; Y1=Y2; X2+=DX; Y2=fn(X2)
    drawV(canvas, X1,Y1,X2,Y2,"red",
          X0=30,Y0=60,Xscale=40,Yscale=20)

表4-3  繰返し法による実根
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 20, 10, 0, 0,
          1, 0.02, 5, 5, 1, 5, 10, 15,"%d","%.2f")
#以下繰り返し法
def fn(X):return (X*X + 1)/4   #この部分に漸化式を書く
def IterativeMethod(iterMax):
    X1 = 0; R=[]; R.append(X1)
    for i in range(iterMax):
        X2 = fn(X1);R.append(X2); X1 = X2
    return R
#以下グラフ化
R=IterativeMethod(10);TV=2-math.sqrt(3)#TV:厳密解
for i in range(1, len(R)):
    drawV(canvas, i-1,TV,i,TV,"blue",width=1,
          X0=30,Y0=170,Xscale=20,Yscale=500)
    drawV(canvas, i-1,R[i-1],i,R[i],"red",
          X0=30,Y0=170,Xscale=20,Yscale=500)

表4-4  繰返し法による実根
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 20, 20, 0, -1,
          1, 0.5, 5, 5, 1, 1, 10, 5,"%d","%.1f")
#■以下繰り返し法
def fn(X):return (-2*X*X-1)/5   #この部分に漸化式を書く
def IterativeMethod(iterMax):
    X1 = 1; R=[]; R.append(X1)
    for i in range(iterMax):
        X2 = fn(X1);R.append(X2); X1 = X2
    return R
def fn2(X): return 1/(2*(X+1))
#■以下グラフ化
R=IterativeMethod(10);TV=2-math.sqrt(3)#TV:厳密解
print(R)
TV=(1+math.sqrt(17))/8
for i in range(1, len(R)):
    drawV(canvas, i-1,TV,i,TV,"#770077",width=1,
          X0=30,Y0=80,Xscale=20,Yscale=40)
    drawV(canvas, i-1,R[i-1],i,R[i],"red",
          X0=30,Y0=80,Xscale=20,Yscale=40)
    drawV(canvas, i-1,fn2(R[i-1]),i,fn2(R[i]),"blue",
          X0=30,Y0=80,Xscale=20,Yscale=40)

表4-5 区間縮小法
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 20, 20, 0, -1,
          1, 1, 5, 5, 1, 1, 10, 6,"%d","%.f")
#■以下区分縮小法
def fn(X):return X*X - X - 2  # ここに解くべき関数を書く
def DiminishingInterval(iterMax,minX, maxM,EPS):
    X1 = minX; X2 = maxM; R=[]
    for i in range(iterMax):
        X = (X1+X2)/2; Y = fn(X)
        R.append([X1,X2,X,Y])#Xが求めるべき値
        print(X,X1,X2,Y)
        if abs(Y) < EPS: return R
        elif Y > 0 : X2 = X
        else       : X1 = X
        E = abs(X2 - X1)
        if E<EPS:return R
    return R
#■以下グラフ化
R=DiminishingInterval(10,1,4, 0.00000001);
print(R)
for i in range(1, len(R)):
    drawV(canvas, i-1,R[i-1][0],i,R[i][0],"red",
          X0=30,Y0=120,Xscale=20,Yscale=20)
    drawV(canvas, i-1,R[i-1][1],i,R[i][1],"blue",
          X0=30,Y0=120,Xscale=20,Yscale=20)
    drawV(canvas, i-1,R[i-1][2],i,R[i][2],"#770077",
          X0=30,Y0=120,Xscale=20,Yscale=20)
    drawV(canvas, i-1,R[i-1][3],i,R[i][3],"#007700",
          X0=30,Y0=120,Xscale=20,Yscale=20)

表4-6 ニュートン・ラプソン法
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 40, 10, 0, -2,
          1, 1, 5, 5, 1, 2, 5, 12,"%d","%.f")
#■以下ニュートン・ラプソン法
def fn(X):return X*X - X - 2  # ここに解くべき関数を書く
def dfn(X):return 2*X -1       # ここに微分式を書く
def NewtonRaphson(iterMax,startX,EPS):
    X = startX; R=[]
    for i in range(iterMax):
        Y = fn(X)
        R.append([X,Y])#Xが求めるべき値
        X=X-Y/dfn(X)
        if abs(Y) < EPS: return R
    return R
#■以下グラフ化
R=NewtonRaphson(10,4, 0.00000001)
print(R)
for i in range(1, len(R)):
    drawV(canvas, i-1,R[i-1][0],i,R[i][0],"red",
          X0=30,Y0=120,Xscale=40,Yscale=10)
    drawV(canvas, i-1,R[i-1][1],i,R[i][1],"blue",
          X0=30,Y0=120,Xscale=40,Yscale=10)

表4-7 セカント法
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 40, 10, 0, -2,
          1, 1, 5, 5, 1, 2, 8, 12,"%d","%.f")
#■以下セカント法
def fn(X):return X*X - X - 2  # ここに解くべき関数を書く
def secant(iterMax,minX,maxX,EPS):
    X1 = minX; X2=maxX; Y1=fn(X1);R=[]
    for i in range(iterMax):
        Y2 = fn(X2);Y=(Y1+Y2)/2
        X=X2-Y2*(X1-X2)/(Y1-Y2)
        R.append([X,Y])#Xが求めるべき値
        if abs(Y) < EPS: return R
        X1=X2; X2=X; Y1=Y2
    return R
#■以下グラフ化
R=secant(10,1, 4, 0.00000001)
print(R)
for i in range(1, len(R)):
    drawV(canvas, i-1,R[i-1][0],i,R[i][0],"red",
          X0=30,Y0=120,Xscale=40,Yscale=10)
    drawV(canvas, i-1,R[i-1][1],i,R[i][1],"blue",
          X0=30,Y0=120,Xscale=40,Yscale=10)

■４．３　連立非線形方程式

表4-8  ニュートン・ラプソン法による連立非線形方程式の解
from graph import *
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240);canvas.pack()
drawScale(canvas, 30,20, 15, 10, 40, 20, 0, -0.4,
          1, 0.2, 2, 2, 1, 1, 4, 8,"%d","%.1f")
#■以下連立方程式用
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):#配列複写---copy.copyを使う方がよいかも
    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 setJACOBI(A):
    JA=array(2,3) #ヤコビアン設定
    JA[0][0] = 2*A[0]+A[1]; JA[0][1] = A[0]+2* A[1]
    JA[1][0] = 1          ; JA[1][1] = -1
    JA[0][2] = -(A[0]*A[0] + A[0]*A[1] + A[1]*A[1] - 1)# - f(X)
    JA[1][2] = -(A[0] - A[1])
    return JA
def nonLinearEq(MAT, iterMax,EPS):
    N=len(MAT);A=array(N)
    for i in range(N): A[i]=MAT[i]
    R=[]; R.append([A[0],-0.4]) #h1初期値=-0.4は適当な値
    for i in range(iterMax):
        JACOBI=setJACOBI(A); B=GaussianELMPivot(JACOBI)
        if B==[]: return R
        E = 0    #補正値の加算，収束判定
        for i in range(N):
            E += B[i][0]*B[i][0]; A[i] = A[i]+B[i][0]
        R.append([A[0],B[i][0]])
        if E<EPS:return R
    return R
#■以下グラフ化
R=nonLinearEq([1,1],10, 0.00000001)
print(R)
for i in range(1, len(R)):
    drawV(canvas, i-1,R[i-1][0],i,R[i][0],"red",
          X0=30,Y0=140,Xscale=40,Yscale=100)
    drawV(canvas, i-1,R[i-1][1],i,R[i][1],"blue",
          X0=30,Y0=140,Xscale=40,Yscale=100)

■４．４　複素数解

表4-9  複素数によるニュートン・ラプソン法
from graph import *
import math
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240)
canvas.pack()
drawScale(canvas, 30,20, 15, 10, 40, 10, 0, -1,
          1, 0.1, 5, 5, 1, 5, 5, 20,"%d","%.1f")
#■以下ニュートン・ラプソン法
def fn(X):return X*X + X + 1  # ここに解くべき関数を書く
def dfn(X):return 2*X + 1     # ここに微分式を書く
def cabs(X): return math.sqrt(X.real*X.real+X.imag*X.imag)
def NewtonRaphson(iterMax,startX,EPS):
    X = startX; R=[]
    for i in range(iterMax):
        Y = fn(X)
        R.append([X,Y])#Xが求めるべき値
        X=X-Y/dfn(X)
        if cabs(Y) < EPS: return R
    return R
#■以下グラフ化
R=NewtonRaphson(10,complex(0,1), 0.00000001)
print(R)
for i in range(1, len(R)):
    drawV(canvas, i-1,R[i-1][0].real,i,R[i][0].real,"red",
          X0=30,Y0=120,Xscale=40,Yscale=100)
    drawV(canvas, i-1,R[i-1][0].imag,i,R[i][0].imag,"#007700",
          X0=30,Y0=120,Xscale=40,Yscale=100)
    drawV(canvas, i-1,R[i-1][1].real,i,R[i][1].real,"blue",
          X0=30,Y0=120,Xscale=40,Yscale=100)
    drawV(canvas, i-1,R[i-1][1].imag,i,R[i][1].imag,"#770077",
          X0=30,Y0=120,Xscale=40,Yscale=100)

表4-10  複素数によるセカント法
from graph import *
import math
root=initTk("Numerical Analysis")
canvas=initCanvas(root,440,240);canvas.pack()
drawScale(canvas, 30,20, 15, 10, 40, 10, 0, -10,
          1, 1, 2, 5, 1, 5, 7, 20,"%d","%.1f")
#■以下, 複素数のセカント法
def fn(X):return X*X + X + 1  # ここに解くべき関数を書く
def cabs(X): return math.sqrt(X.real*X.real+X.imag*X.imag)
def secant(iterMax,minX,maxX,EPS):
    X1 = minX; X2=maxX; Y1=fn(X1);R=[]
    for i in range(iterMax):
        Y2 = fn(X2); Y=(Y1+Y2)/2
        X=X2-Y2*(X1-X2)/(Y1-Y2)
        R.append([X,Y])#Xが求めるべき値
        if cabs(Y) < EPS: return R
        X1=X2; X2=X; Y1=Y2
    return R
#■以下グラフ化
R=secant(10,complex(0,1),complex(0,2), 0.00000001)
print(R)
for i in range(1, len(R)):
    drawV(canvas, i-1,R[i-1][0].real,i,R[i][0].real,"red",
          X0=30,Y0=120,Xscale=40,Yscale=40)
    drawV(canvas, i-1,R[i-1][0].imag,i,R[i][0].imag,"#007700",
          X0=30,Y0=120,Xscale=40,Yscale=40)
    drawV(canvas, i-1,R[i-1][1].real,i,R[i][1].real,"blue",
          X0=30,Y0=120,Xscale=40,Yscale=40)
    drawV(canvas, i-1,R[i-1][1].imag,i,R[i][1].imag,"#770077",
          X0=30,Y0=120,Xscale=40,Yscale=40)

表4-11  ミュラー法
import math, cmath
#■以下, Muller法
def fn(X):return X*X + X + 1  # ここに解くべき関数を書く
def cabs(X): return math.sqrt(X.real*X.real+X.imag*X.imag)
def root2EQ(A, B, C): #根データ
    A2 = A* 2.0; B2 = -B; BAC = cmath.sqrt(B * B - 4 * A * C)
    return [(B2 + BAC)/A2, (B2 - BAC)/A2]
def Muller(iterMax, EPS):
    X1=complex(0,0); X2=complex(1,0); X3=complex(2,0)#初期値
    F1 = fn(X1); F2 = fn(X2); F3 = fn(X3)
    DX1 = X1 - X2; DX2 = X2 - X3; DX3 = X3 - X1
    DF1 = F1 - F2; DF2 = F2-F3
    for Iter in range(iterMax):
        E = cabs(F3)
        if E < EPS : break
        # 連立方程式を解く
        XX = DF2*DX1-DF1*DX2; YY = DX1*DX2*DX3;
        A = XX / YY; B = DF1/DX1- A*(X1+X2); C = F1-B*X1-A*X1*X1
        R = root2EQ(A, B, C)#二次方程式の根
        X1 = X2; X2 = X3; F1 = F2; F2 = F3
        DX1 = DX2; DX2 = DX3; DF1 = DF2
        E1 = cabs(R[0]-X3); E2 = cabs(R[1]-X3)
        if E1 <= E2 : X3 = R[0]
        else:         X3 = R[1]
        F3 = fn(X3)
        DX3 = X3-X1; DF2 = F2-F3
    return R
print(Muller(10,0.00000001))

表4-12  ニュートン・ラプソン法による複数の解
import math, cmath
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 printCarray(L, form="%9.5f"):
    for i in range(len(L)):
        RE=L[i].real; IM=L[i].imag
        if   IM==0: S = form % RE
        elif RE==0: S= (form % IM)+" i"
        elif IM>0 : S= form % RE + " +" + (form % IM)+" i"
        else: S= form % RE + " -" + (form % abs(IM))+" i"
        print(S)
#■以下, ニュートン法による複数の解
def fn(X):return X*X*X + 2*X*X + 2*X+ 1  # ここに解くべき関数を書く
def dfn(X):return 3*X*X +4*X+2           # ここに微分式を書く
def cabs(X): return math.sqrt(X.real*X.real+X.imag*X.imag)
def updDiff(DF,F,X,k,Z):
    if k >= 2 :
        DFD = complex(0, 0) # F'-F*Σ(1/(X-Zi)) (i=1〜N-1)の計算
        for i in range(k-1): DFD +=1/(X-Z[i])
        DF-=F*DFD
    return DF
def multiComplexNewton(N):
    Z=array(N)
    EPS = 0.00001; iterMax = 500
    for k in range(N):
        X = complex(1, 1)
        for Iter in range(iterMax):
            F = fn(X); E = cabs(F)
            if E<EPS: break
            DF = updDiff(dfn(X), F, X, k,Z); X-= F/DF
        Z[k] = X
    return Z
printCarray(multiComplexNewton(3))

表4-13  DKA法による複数の解
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 printCarray(L, form="%9.5f"):
    for i in range(len(L)):
        RE=L[i].real; IM=L[i].imag
        if   IM==0: S = form % RE
        elif RE==0: S= (form % IM)+" i"
        elif IM>0 : S= form % RE + " +" + (form % IM)+" i"
        else: S= form % RE + " -" + (form % abs(IM))+" i"
        print(S)
#----------------------以下、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 maxFx<delta : return X
    print("収束しませんでした")
    return X
#A=[1, -9, -7, -35, 50] #Aには，次元数+1個の値を設定
A=[1, 0, 1, 0, 1]
printCarray(DKA(A, 0.001, 500))

表4-14  ベアストウ法による解
import math	
import math, copy
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 root(P,Q):#■二次方程式の解
    D=P*P-4*Q              #判別式
    if D>=0:
        V=math.sqrt(D)
        if P<0:r1=(-P+V)/2 #減算による桁落ち対策
        else:  r1=(-P-V)/2
        r2 = Q / r1
        return [complex(r1,0),complex(r2,0)]
    else:
        V=math.sqrt(-D)    #虚数項
        return [complex(-P/2,V/2),complex(-P/2,-V/2)]
def bairstow(a,n):#■ベアストウ法によるPQ計算
    B=array(n+1); C=array(n+1);EPS=0.000001
    P=1; Q=1; DP=1; DQ=1   #最初の差は大きな値とする
    while abs(DP)>EPS or abs(DQ)>EPS:#収れんの判定
        B[0] = A[0]; B[1] = A[1] - P*B[0]
        for k in range(2,n+1): B[k] = A[k] - P*B[k-1] - Q*B[k-2]
        C[0] = B[0]; C[1] = B[1] - P*C[0]
        for k in range(2,n+1): C[k] = B[k] - P*C[k-1] - Q*C[k-2]
        e =C[n-2 ]*C[n-2] - C[n-3]*(C[n-1]-B[n-1])
        DP=(B[n-1]*C[n-2] - B[n  ]*C[n-3])/e
        DQ=(B[n  ]*C[n-2] - B[n-1]*(C[n-1]-B[n-1]))/e
        P+=DP; Q+=DQ
    return [P,Q,B]
def solEq(A,n):#■Ｎ次方程式の解
    A=copy.copy(A); n=len(A)-1; r=[]
    while n>0:
        if n==1: r.append(complex(-A[1]/A[0],0))#一次式
        if n==2:                                #二次式
            rr=root(A[1],A[2])
            r.append(rr[0]); r.append(rr[1])
        else:                                   #Ｎ次式
            rr=bairstow(A,n); A=rr[2]
            rr=root(rr[0],rr[1])
            r.append(rr[0]); r.append(rr[1])
        n-=2
    return r
def printCL(A, form="%8.5f"):
    for i in range(len(A)):
        if   A[i].imag==0: S= form % A[i].real
        elif A[i].imag>0 : 
S= form % A[i].real+" +"+form % A[i].imag+"i"
        else: S=form % A[i].real+" -"+form % abs(A[i].imag)+"i"
        print(S)
A=[1,-8,-6,-32,40]
printCL(solEq(A,4))