例題8-1 例題8-2 依牛頓方法解 x^2 + y^2 - x =0 x^2 - y^2 - y =0

例題8-1  例題8-2
依牛頓方法解
x^2 + y^2 - x =0
x^2 - y^2  - y =0

TOL=1X10^-7  , x0=0.76 , y0=0.41

'''
/* ex8-2.py based on Newton's method
 * for solving the n x n system of nonlinear
 * algebra equations
 * f1(x1,x2,...,xn)=0
 * f2(x1,x2,...,xn)=0
 * ...
 * fn(x1,x2,...,xn)=0
 * Output solution x1,x2,x3,...,xn.
 */
 '''
import math

MAX =  20
def f(x,y):
     return (math.pow(x,2)+math.pow(y,2)-x)
 
def g(x,y):
    return (math.pow(x,2)-math.pow(y,2)-y)

def fx(x,y):
    return (2*x-1)

def fy(x,y):
    return (2*y)

def gx(x,y):
    return (2*x)

def gy(x,y):
     return (-2*y-1)
 

def gaussh(n, a ,  x ):
    for k in range (1,n):
        #/* check if a[k][k]=0 is true then interchange */
        #/* E(k) and E(k+1).............................*/
        if(a[k][k]==0):
            for m in range (1, n+2):
                temp=a[k][m]
                a[k][m]=a[k+1][m]
                a[k+1][m]=temp

        #/* To reduce the matrix to triangular form */
        for i in range (k , n):
            bb=a[i+1][k]/a[k][k];
            for j in range (k , n+2):
                a[i+1][j]=a[i+1][j]-bb*a[k][j];

    if(abs(a[n][n])==0.0):
        print("NO UNIQUE SOLUTION!!!\n");
        exit(1);
    #/* To start backward substitution */
    x[n]=a[n][n+1]/a[n][n];
    for i in range (n-1 , 0, -1):
        cc=0.0
        for j in range (i+1, n+1):
            cc=cc+a[i][j]*x[j]
        x[i]=(a[i][n+1]-cc)/a[i][i]

#++++++++++++++++++++++++++++++
TOL = 0.00000001
a=  [ [0.0 for i in range(11)]  for j in range(11) ]       #[n][n],
dxy= [0.0 for i in range(11)]
n=2
i=1
x=0.76
y=0.41
print(" i\t\tx\t\t\t y\t\t\t  f(x,y)\t\t\t  g(x,y)\n")
while(i<=MAX):
    a[1][1]=fx(x,y);a[1][2]=fy(x,y);a[1][3]=-f(x,y);
    a[2][1]=gx(x,y);a[2][2]=gy(x,y);a[2][3]=-g(x,y);
    print("{%2d} \t\t{%9.6f} \t\t{%9.6f} \t\t {%10.7f} \t\t{%10.7f}" %(i-1,x,y,f(x,y),g(x,y)))
    gaussh(n,a,dxy)   #/* call the function gaussh() */

    if(abs(f(x,y))<TOL  and  abs(g(x,y))<TOL):
            break
    x=x+dxy[1]
    y=y+dxy[2]
    i+=1

print("\nroot_x={%9.6f}  root_y={%9.6f}\n" %(x,y))

輸出畫面
========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX8-2.py ============
 i x  y  f(x,y)             g(x,y)

{ 0} { 0.760000} { 0.410000} {-0.0143000} {-0.0005000}
{ 1} { 0.772056} { 0.419794} { 0.0002413} { 0.0000494}
{ 2} { 0.771845} { 0.419643} { 0.0000001} { 0.0000000}
{ 3} { 0.771845} { 0.419643} { 0.0000000} { 0.0000000}

root_x={ 0.771845}  root_y={ 0.419643}

>>>