例題8-3 利用牛頓法求解非線性聯立方程式 x1 x2 x3 - x1^2 + x2^2 = 1.34 x1 x2 - x3^2 = 0.09 exp(x1) - exp(x2) + x3 = 0.41

例題8-3 利用牛頓法求解非線性聯立方程式

x1 x2 x3 -  x1^2 + x2^2 = 1.34
x1 x2 - x3^2 = 0.09
exp(x1) - exp(x2) + x3 = 0.41

初始值為(1,1,1)  ,  TOL=10^(-8)


'''
/* ex8-3.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
 * Input number of unknowns and equations n
 * with coefficent a11,a12,...,ann and b1,b2,
 * ...bn. Output solution x1,x2,x3,...,xn.
 */
'''
import math
def f1(x1,x2,x3):
    return (x1*x2*x3-math.pow(x1,2)+math.pow(x2,2)-1.34)

def f2(x1,x2,x3):
    return (x1*x2-math.pow(x3,2)-0.09)

def f3(x1,x2,x3):
    return (math.exp(x1)-math.exp(x2)+x3-0.41)

def f1x1(x1,x2,x3):
    return (x2*x3-2*x1)

def  f1x2(x1,x2,x3):
    return (x1*x3+2*x2)

def f1x3(x1,x2,x3):
    return (x1*x2)

def f2x1(x1,x2,x3):
    return (x2)

def f2x2(x1,x2,x3):
    return (x1)

def f2x3(x1,x2,x3):
    return (-2*x3)

def f3x1(x1,x2,x3):
    return (math.exp(x1))

def f3x2(x1,x2,x3):
    return (-math.exp(x2))

def f3x3(x1,x2,x3):
    return (1)

TOL = 0.00000001
MAX = 20

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]

#++++++++++++++++++++++++++++++
a=  [ [0.0 for i in range(11)]  for j in range(11) ]       #[n][n],
dx= [0.0 for i in range(11)]
n=3;
i=1;
x1=1.0;
x2=1.0;
x3=1.0;
print(" i\tx1\t\tx2\t\tx3\t\tf1()\t\tf2()\t\t f3()\n")
while(i<=MAX):
    a[1][1]=f1x1(x1,x2,x3)
    a[1][2]=f1x2(x1,x2,x3)
    a[1][3]=f1x3(x1,x2,x3)
    a[1][4]=-f1(x1,x2,x3)
    a[2][1]=f2x1(x1,x2,x3)
    a[2][2]=f2x2(x1,x2,x3)
    a[2][3]=f2x3(x1,x2,x3)
    a[2][4]=-f2(x1,x2,x3)
    a[3][1]=f3x1(x1,x2,x3)
    a[3][2]=f3x2(x1,x2,x3)
    a[3][3]=f3x3(x1,x2,x3)
    a[3][4]=-f3(x1,x2,x3)
    print("{%2d}\t{%11.8f}\t{%11.8f}\t{%11.8lf}" %( i-1,x1,x2,x3),end='')
    print("\t{%11.8f}\t{%11.8f}\t{%11.8f}" %(f1(x1,x2,x3),f2(x1,x2,x3),f3(x1,x2,x3)))
    gaussh(n,a,dx)    # /* call the function gaussh() */
    if(abs(f1(x1,x2,x3))<TOL  and  abs(f2(x1,x2,x3))<TOL and abs(f3(x1,x2,x3))<TOL):
        break
    x1=x1+dx[1];
    x2=x2+dx[2];
    x3=x3+dx[3];
    i+=1
 
#bound:
print("\nroot_x1={%10.8f}\n" %(x1))
print("root_x2={%10.8f}\n" %(x2))
print("root_x3={%10.8f}\n" %(x3))

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX8-3.py ============
 i x1 x2 x3 f1() f2()              f3()

{ 0} { 1.00000000} { 1.00000000} { 1.00000000} {-0.34000000} {-0.09000000} { 0.59000000}
{ 1} { 0.89626021} { 1.09518003} { 0.95072012} {-0.01066825} {-0.01230246} { 0.00142115}
{ 2} { 0.90223257} { 1.10035236} { 0.95012799} { 0.00001382} { 0.00003054} { 0.00000373}
{ 3} { 0.90221844} { 1.10034319} { 0.95013153} {-0.00000000} { 0.00000000} { 0.00000000}

root_x1={0.90221844}

root_x2={1.10034319}

root_x3={0.95013153}

>>>