NonLineearSystemEquations

2018年12月02日

Numerical solutions of nonlinear systems of equations.

A system of nonlinear equations has the form

f1(z1,x2,x3,…..xn)=0

f2(z1,x2,x3,…..xn)=0

fn(z1,x2,x3,…..xn)=0

The system nonlinear equations can be solved by Newton’s method or Steepest

Decent method( please refer “Numerical analysis” by Richard L. Burdoen & J.Douglas Faires). Newton’s method or Steepest Decent methods all need an initial approximate solution. Because the Newton’s method can converge very rapidly when the approximation is near the true solution, while Steepest Decent method can get an approximation solution for any poor initial approximation. In this article we use Steepest Decent to get the first approximation and then use the solution as the initial solution of Nwton’s method to solve the non-linear equation. The snippet code of Newton’s and Steepest Decent method are modification from the code listed in “Numerical Analysis”.

Public Sub SteepEstMethod(canvas As PictureBox, AAA As String, nEq As Integer, stEqs() As String, VvarsSt() As String, VvarsI() As Double, NTry As Integer, Tol As Double, XAns() As Double)

'C***********************************************************************

'C *

'C STEEPEST DESCENT ALGORITHM 10.3 *

'C *

'C***********************************************************************

'C TO APPROXIMATE A SOLUTION P TO THE MINIMIZATION PROBLEM

'C G(P) = MIN ( G(X) : X IN R(n) )

'C GIVEN AN INITIAL APPROXIMATION X:

'C INPUT NUMBER N OF VARIABLES; INITIAL APPROXIMATION X;

'C TOLERANCE TOL; MAXIMUM NUMBER OF ITERATIONS N.

'C OUTPUT APPROXIMATE SOLUTION X OR A MESSAGE OF FAILURE.

Dim FLAG As Boolean

Dim Xx() As Double, Z() As Double, C() As Double, G() As Double, A() As Double

ReDim Xx(1 To nEq), Z(1 To nEq), C(1 To nEq), G(1 To nEq), A(1 To nEq)

'Dim FLAG1 As Boolean

Dim OK As Boolean, k As Integer, i As Integer, j As Integer

Dim Z0 As Double, X0 As Double, G0 As Double, G1 As Double, A0 As Double

Dim Zero As Double

Dim H1 As Double, H2 As Double, H3 As Double

'C THE FUNCTION F USED HERE CORRESPONDS TO G USED IN THE ALGORITHM

'pi = 4 * Atan(1#)

Zero = 1E-20

'OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')

'OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')

'WRITE(6,*) 'This is the Steepest Descent Method'

'WRITE(6,*) 'NOTE THAT THE FUNCTIONS ARE VERY COMPLICATED.

'WRITE(6,*) 'Have the functions F1,...,Fn and F (G) been'

'WRITE(6,*) 'defined at the beginning of the program and.'

'WRITE(6,*) 'has the gradient been implemented at line 118.'

'WRITE(6,*) 'Enter Y or N '

'WRITE(6,*) ' '

'READ(5,*) AAA

If LCase(AAA) = "y" Then

OK = False

101 If OK Then GoTo 11

'WRITE(6,*) 'Input the number n of equations.'

'WRITE(6,*) ' '

'READ(5,*) N

If nEq >= 2 Then

OK = True

Else

MsgBox ("Neq must be greater than 1.")

End If

GoTo 101

11 OK = False

14 If (OK) Then GoTo 15

'WRITE(6,*) 'Input tolerance '

'WRITE(6,*) ' '

'READ(5,*) TOL

If Tol <= 0# Then

MsgBox ("Tolerance must be positive ")

Else

OK = True

End If

GoTo 14

15 OK = False

12 If OK Then GoTo 13

'WRITE(6,*) 'Input the maximum number of iterations.'

'WRITE(6,*) ' '

'READ(5,*) M

If (NTry <= 0) Then

MsgBox ("Niteration Must be positive integer ")

Else

OK = True

End If

GoTo 12

'13 For I = 1 To NTry

'WRITE(6,*) 'Input initial approximation X(',I,')'

'WRITE(6,*) ' '

'READ(5,*) X(I)

'16 CONTINUE

Else

MsgBox ("The program will end so that the functions can be created ")

OK = False

End If

If Not (OK) Then GoTo 400

'WRITE(6,*) 'Select output destination: '

'WRITE(6,*) '1. Screen '

'WRITE(6,*) '2. Text file '

' WRITE(6,*) 'Enter 1 or 2 '

'WRITE(6,*) ' '

'READ(5,*) FLAG1

canvas.Print "STEEPEST DESCENT METHOD FOR NONLINEAR SYSTEMS"

canvas.Print "Iter. No. Vector G(vector)"

'C STEP 1

k = 0

For i = 1 To nEq 'DO 16 I=1,NEq

'WRITE(6,*) 'Input initial approximation X(',I,')'

'WRITE(6,*) ' '

'READ(5,*) X(I)

Xx(i) = VvarsI(i)

' canvas.Print "i=" & i & "x(" & i & ")=" & Xx(i)

16 Next i ' CONTINUE

'canvas.Print (" K= " & k & ";X(1)= " & X(1) & ";X(2)= " & X(2) & ";X(3)= " & X(3))

k = 1

'C STEP 2

100 If (k > NTry) Then GoTo 200

'MsgBox ("kstep = " & K)

'C STEP 3

'////////////////////////////////

'Dim Fxyz1 As Double, Fxyz2 As Double, Fxyz3 As Double, fxyz As Double

'Fxyz1 = 3 * Xx(1) - Cos(Xx(2) * Xx(3)) - 0.5 '在點(x1,x2,x3)處方程式1值

'Fxyz2 = Xx(1) * Xx(1) - 81 * (Xx(2) + 0.1) ^ 2 + Sin(Xx(3)) + 1.06

'Fxyz3 = Exp(-Xx(1) * Xx(2)) + 20 * Xx(3) + (10 * PI - 3) / 3

'F(X1,X2,X3)=F1(X1,X2,X3)**2+F2(X1,X2,X3)**2+F3(X1,X2,X3)**2 '在點(x1,x2,x3)處方程式1,方程式2,方程式3平方和

'fxyz = Fxyz1 ^ 2 + Fxyz2 ^ 2 + Fxyz3 ^ 2

'canvas.Print "fxyz1~3 A=" & Fxyz1 & ";" & Fxyz2 & ";" & Fxyz3 & ";f^2 A=" & fxyz

'//////////////////////////////

'C(1) = F1(X(1), X(2), X(3)) '在點(x1,x2,x3)處方程式1值

'C(2) = F2(X(1), X(2), X(3)) '在點(x1,x2,x3)處方程式2值

'C(3) = F3(X(1), X(2), X(3)) '在點(x1,x2,x3)處方程式3值

'F(X1,X2,X3)=F1(X1,X2,X3)**2+F2(X1,X2,X3)**2+F3(X1,X2,X3)**2 '在點(x1,x2,x3)處方程式1,方程式2,方程式3平方和

Dim sumF As Double

sumF = 0

For i = 1 To nEq

C(i) = FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, Xx)

sumF = sumF + C(i) * C(i)

Next i

'canvas.Print "fxyz1~3 B =" & C(1) & ";" & C(2) & ";" & C(3) & ";f^2 B=" & sumF

'C THE VECTOR C HOLDS F(X) FROM THE EXAMPLE AND G1 HOLDS G(X)

' G1 = F(X(1), X(2), X(3))

G(1) = sumF 'G(1) = G1

'Dim Q As Double, P As Double, Z1 As Double, Z2 As Double, Z3 As Double

'Q = Exp(-Xx(1) * Xx(2))

'P = Sin(Xx(2) * Xx(3))

'Z1 = 6 * c(1) + 4 * Xx(1) * c(2) - 2 * Xx(2) * Q * c(3)

'Z2 = 2 * Xx(3) * P * c(1) - 324 * (Xx(2) + 0.1) * c(2) - 2 * Xx(1) * Q * c(3)

'Z3 = 2 * Xx(2) * P * c(1) + 2 * Cos(Xx(3)) * c(2) + 40 * c(3)

' canvas.Print "Z1,Z2,Z3 = " & Z1 & ";" & Z2 & ";" & Z3

Dim vTpt As Double, zoSum2 As Double

zoSum2 = 0

For i = 1 To nEq

vTpt = 0

For j = 1 To nEq

vTpt = vTpt + 2 * C(j) * PartialDif(glVbScriptA(j - 1), stEqs(j), VvarsSt, Xx, i)

Next j

Z(i) = vTpt

zoSum2 = zoSum2 + Z(i) * Z(i)

Next i

' canvas.Print "Z(1),Z(2),Z(3) =" & Z(1) & ";" & Z(2) & ";" & Z(3)

vTpt = 0

'Dim Z0A As Double

'Z0A = Sqr(Z1 ^ 2 + Z2 ^ 2 + Z3 ^ 2)

Z0 = Sqr(zoSum2)

'?canvas.Print "Z0A,Z0=" & Z0A & ";" & Z0

'C STEP 4

If (Abs(Z0) < 1E-20) Then

MsgBox ("ZERO GRADIENT: MAY HAVE MIN")

GoTo 400

End If

'C STEP 5

For j = 1 To nEq 'DO 10 J=1,N

10 Z(j) = Z(j) / Z0

Next j

'C A IS USED FOR ALPHA AND G0 FOR G(3)

A(1) = 0

X0 = 1

'????????????????????????????????????////

'Dim G0a As Double

'G0a = Fxyz_T(Xx(1) - X0 * Z(1), Xx(2) - X0 * Z(2), Xx(3) - X0 * Z(3)) '在點 x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3)處方程式1,方程式2,方程式3平方和

Dim xXtpt() As Double

ReDim xXtpt(1 To nEq)

For i = 1 To nEq

xXtpt(i) = Xx(i) - X0 * Z(i)

Next i

sumF = 0

For i = 1 To nEq

sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2

Next i

G0 = sumF

'?canvas.Print "g0a,g0=" & G0a & ";" & G0

'????????????????????????????????????????

'C STEP 6

FLAG = True

If (Abs(G0) < Abs(G1)) Then FLAG = False

500 If Not (FLAG) Then GoTo 600

'C STEPS 7 AND 8

X0 = X0 / 2

' MsgBox ("x0=" & X0 & "zero= " & Zero)

If (X0 < Zero) Then MsgBox ("NO LIKELY IMPROVEMENT: MAY HAVE MIN")

If (X0 < Zero) Then GoTo 400

'G0 = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))

'????????????????????????????????????////

'Dim G0b As Double

'G0b = Fxyz_T(Xx(1) - X0 * Z(1), Xx(2) - X0 * Z(2), Xx(3) - X0 * Z(3))

For i = 1 To nEq

xXtpt(i) = Xx(i) - X0 * Z(i)

Next i

'G0 = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3)) '在點 x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3)處方程式1,方程式2,方程式3平方和

sumF = 0

For i = 1 To nEq

sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2

Next i

G0 = sumF

' canvas.Print "g0b,g0=" & G0b & ";" & G0

'????????????????????????????????????????

If (Abs(G0) < Abs(G1)) Then FLAG = False

If (Abs(G0) < Abs(G1)) Then GoTo 500

600 A(3) = X0

G(3) = G0

'C STEP 9

X0 = X0 / 2

'G(2) = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))

For i = 1 To nEq

xXtpt(i) = Xx(i) - X0 * Z(i)

Next i

sumF = 0

For i = 1 To nEq

sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2

Next i

G(2) = sumF

A(2) = X0

'C STEP 10

H1 = (G(2) - G(1)) / (A(2) - A(1))

H2 = (G(3) - G(2)) / (A(3) - A(2))

H3 = (H2 - H1) / (A(3) - A(1))

'C STEP 11

X0 = 0.5 * (A(1) + A(2) - H1 / H3)

'G0 = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))

'????????????????????????????????????////

For i = 1 To nEq

xXtpt(i) = Xx(i) - X0 * Z(i)

Next i

sumF = 0

For i = 1 To nEq

sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2

Next i

G0 = sumF

'????????????????????????????????????????

'C STEP 12

A0 = X0

For j = 1 To nEq 'DO 20 J=1,N

If (Abs(G(j)) < Abs(G0)) Then

A0 = A(j)

G0 = G(j)

End If

20 'CONTINUE

Next j

'MsgBox ("pass 20;k= " & K)

If Abs(A0) < Zero Then

MsgBox ("NO CHANGE LIKELY: PROBABLY ROUNDING DIFFICULTIES")

GoTo 400

End If

'C STEP 13

For j = 1 To nEq 'DO 30 J=1,N

30 Xx(j) = Xx(j) - A0 * Z(j) '

' canvas.Print " K= " & K & ";X(j)= " & Xx(j)

Next j

'canvas.Print " K= " & K & ";X(1)= " & Xx(1) & ";X(2)= " & Xx(2) & ";X(3)= " & Xx(3)

'C STEP 14

If (Abs(G0) < Tol) Or (Abs(G0 - G1) < Tol) Then

MsgBox ("PROCEDURE COMPLETED SUCCESSFULLY")

GoTo 400

End If

'C STEP 15

k = k + 1

'MsgBox ("k=" & K)

GoTo 100

'C STEP 16

Exit Sub

200 MsgBox ("MAXIMUM ITERATIONS EXCEEDED")

For i = 1 To nEq

ReDim Preserve XAns(1 To i)

XAns(i) = 0 'Xx(i)

Next i

400

canvas.Print "sta 400 K= " & k & ";X(1)= " & Xx(1) & ";X(2)= " & Xx(2) & ";X(3)= " & Xx(3)

For i = 1 To nEq

ReDim Preserve XAns(1 To i)

XAns(i) = Xx(i)

Next i

End Sub

Public Sub NewtonMethod(canvas As PictureBox, AAA As String, nEq As Integer, stEqs() As String, VvarsSt() As String, VvarsI() As Double, NTry As Integer, Tol As Double, XAns() As Double)

'C***********************************************************************

'C *

'C NEWTON 'S METHOD FOR SYSTEMS ALGORITHM 10.1 *

'C *

'C***********************************************************************

'C TO APPROXIMATE THE SOLUTION OF THE NONLINEAR SYSTEM F(X)=0 GIVEN

'C AN INITIAL APPROXIMATION X:

'C INPUT: NUMBER n OF EQUATIONS AND UNKNOWNS; INITIAL APPROXIMATION

'C X=(X(1),...,X(n)); TOLERANCE TOL; MAXIMUM NUMBER OF

'C ITERATIONS N.

'C OUTPUT: APPROXIMATE SOLUTION X=(X(1),...,X(n)) OR A MESSAGE

'C THAT THE NUMBER OF ITERATIONS WAS EXCEEDED.

'C INITIALIZATION

Dim AA() As Double, yY() As Double, Xx() As Double

ReDim AA(1 To nEq, 1 To nEq + 1), yY(1 To nEq), Xx(1 To nEq)

Dim i As Integer, k As Integer, j As Integer

Dim R

'C USE AA FOR J; NN FOR MAXIMUM NUMBER OF ITERATIONS N

Dim OK As Boolean

'OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')

'OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')

'WRITE(6,*) 'This is the Newton Method for Nonlinear Systems.'

'WRITE(6,*) 'Has the function F been defined and has the'

'WRITE(6,*) 'Jacobian of partial derivatives been defined'

'WRITE(6,*) 'at lines 105 and 115 respectively.'

''WRITE(6,*) 'Enter Y or N '

'WRITE(6,*) ' '

'READ(5,*) AAA

If (LCase(AAA) = "y") Then

OK = False

101 If (OK) Then GoTo 11

'WRITE(6,*) 'Input the number n of equations.'

'WRITE(6,*) ' '

'EAD(5,*) NEq

If (nEq >= 2) Then

OK = True

Else

MsgBox ("NEq must be greater than 1.")

End If

GoTo 101

11 OK = False

14 If OK Then GoTo 15

'WRITE(6,*) 'Input tolerance '

'WRITE(6,*) ' '

'READ(5,*) TOL

If (Tol <= 0#) Then

MsgBox ("Tolerance must be positive ")

Else

OK = True

End If

GoTo 14

15 OK = False

12 If OK Then GoTo 13

' WRITE(6,*) 'Input the maximum number of iterations.'

'WRITE(6,*) ' '

'READ(5,*) Ntry

If (NTry <= 0) Then

MsgBox ("NTry Must be positive integer ")

Else

OK = True

End If

GoTo 12

13 For i = 1 To nEq 'DO 16 I=1,NEq

'WRITE(6,*) 'Input initial approximation X(',I,')'

'WRITE(6,*) ' '

'READ(5,*) X(I)

Xx(i) = VvarsI(i)

canvas.Print "i=" & i & "x(" & i & ")=" & Xx(i)

16 Next i ' CONTINUE

Else

MsgBox ("program will end so that the functions F(X) and J(X) can be created ")

OK = False

End If

If Not OK Then GoTo 400

canvas.Print "NEWTON METHOD FOR NONLINEAR SYSTEMS"

'Pi = 4 * ATAN(1#)

'C STEP 1

k = 1

'C STEP 2

100 If (k > NTry) Then GoTo 200

'C STEP3

'C COMPUTE J(X)

'AA(1, 1) = 3#

'AA(1, 2) = x(3) * Sin(x(2) * x(3))

'AA(1, 3) = x(2) * Sin(x(2) * x(3))

'AA(2, 1) = 2 * x(1)

'AA(2, 2) = -162 * (x(2) + 0.1)

'AA(2, 3) = Cos(x(3))

'AA(3, 1) = -x(2) * Exp(-x(1) * x(2))

'AA(3, 2) = -x(1) * Exp(-x(1) * x(2))

'AA(3, 3) = 20

For i = 1 To nEq

' MsgBox ("steq= " & stEqs(i) & " ;xvst= " & VvarsSt(i) & " ; xv= " & xX(i))

For j = 1 To nEq

AA(i, j) = PartialDif(glVbScriptA(i - 1), stEqs(i), VvarsSt, Xx, j)

'canvas.Print "AA(" & I & "," & J & ")= " & AA(I, J)

Next j

Next i

'C COMPUTE - F(X)

'AA(1, 4) = -(3 * x(1) - Cos(x(2) * x(3)) - 0.5)

'AA(2, 4) = -(x(1) ^ 2 - 81 * (x(2) + 0.1) ^ 2 + Sin(x(3)) + 1.06)

'AA(3, 4) = -(Exp(-x(1) * x(2)) + 20 * x(3) + (10 * PI - 3) / 3)

For i = 1 To nEq

AA(i, nEq + 1) = -FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, Xx)

'canvas.Print "AA(" & I & "," & (nEq + 1) & ")= " & AA(I, nEq + 1)

Next i

'C STEP 4

'C SOLVES THE N X N LINEAR SYSTEM J(X) Y = -F(X)

Call LIN(nEq, nEq + 1, AA, yY)

'C STEPS 5 AND 6

'C R = INFINITY NORM OF Y

R = 0

For i = 1 To nEq 'DO 10 I=1,N

If (Abs(yY(i)) > R) Then R = Abs(yY(i))

' canvas.Print "前k=" & K & "; x(" & I & ")=" & xX(I)

10 Xx(i) = Xx(i) + yY(i)

'canvas.Print "後k=" & K & "; x(" & I & ")=" & xX(I)

Next i

canvas.Print "IterNo=" & k & "Error=" & R

'Form1.Caption = "IterNo=" & K & " ;Error=" & R

For i = 1 To nEq

canvas.Print "x(" & i & ")= " & Xx(i)

ReDim Preserve XAns(1 To i)

XAns(i) = Xx(i)

Next i

'C STEP 6

If (R < Tol) Then

'C PROCESS Is COMPLETE

'?MsgBox ("SUCCESS WITHIN TOLERANCE 1.0E-05=" & R)

GoTo 400

End If

'C STEP 7

k = k + 1

GoTo 100

Exit Sub

'C STEP 8

'C DIVERGENCE

200 MsgBox ("DIVERGENCE - STOPPED AFTER ITERNN" & NTry)

400

End Sub

Sub LIN(n As Integer, M As Integer, A() As Double, x() As Double)

'C LIN SOLVES THE LINEAR SYSTEM: J(X) Y = -F(X) USING PARTIAL

'C PIVOTING

ReDim Preserve A(1 To n, 1 To M), x(1 To n)

Dim i As Integer

Dim Zero As Double, L As Integer, k As Integer, IR As Integer, IA As Integer, j As Integer, JA As Integer

Dim y As Double, C As Double

Zero = 1E-20

k = n - 1

For i = 1 To k 'DO 10 I=1,K

y = Abs(A(i, i))

IR = i

IA = i + 1

For j = IA To n ' DO 20 J=IA,N

If (Abs(A(j, i)) > y) Then

IR = j

y = Abs(A(j, i))

End If

20 Next j ' CONTINUE

If (y <= Zero) Then

MsgBox ("LINEAR SYSTEM HAS NO SOLUTION") 'WRITE(6,1)

Exit Sub

End If

If (IR <> i) Then

For j = i To M 'DO 30 J=I,M

C = A(i, j)

A(i, j) = A(IR, j)

30 A(IR, j) = C

Next j

End If

For j = IA To n ' DO 40 J=IA,N

C = A(j, i) / A(i, i)

If (Abs(C) <= Zero) Then C = 0

For L = i To M 'DO 40 L=I,M

40 A(j, L) = A(j, L) - C * A(i, L)

Next L

Next j

10 Next i 'CONTINUE

If (Abs(A(n, n)) < 1E-20) Then

MsgBox ("LINEAR SYSTEM HAS NO SOLUTION") 'WRITE(6,1)WRITE(6,1)

Exit Sub

End If

x(n) = A(n, M) / A(n, n)

For i = 1 To k 'DO 50 I=1,K

j = n - i

JA = j + 1

C = A(j, M)

For L = JA To n 'DO 60 L=JA,N

60 C = C - A(j, L) * x(L)

Next L

50 x(j) = C / A(j, j)

Next i

'1 FORMAT(1X,'LINEAR SYSTEM HAS NO SOLUTION')

End Sub

VB_VB NET: chday169

上次修改此網站的日期： 2018年12月02日