www.pudn.com > 牛顿法解方程之混沌情况1.32源代码.zip > Module2.bas
Attribute VB_Name = "Module2"
Option Explicit
'======================================================================
'实现不同方程的牛顿迭代法 和 一些基本复变函数
'======================================================================
Public Const PI As Double = 3.14159265358979
Public Const e As Double = 2.718281828
'+++++++++++实现不同方程的牛顿迭代法，并返回方程的各种基本性质+++++++++++++
Function MMi(x0 As Double, y0 As Double, nN As Long, M As Long, _
nM As Long, Hx As Double, Hy As Double, dL1 As Double, _
dL2 As Double, dL3 As Double, dL4 As Double) As Long
'调用时MMi(x0, y0, Int(SeData(0, 13)), M, nM, dX, dY, dL1, dL2, dL3, dL4)
' 牛顿迭代法解方程原理：
' 不失一般性设方程为： f(Z) = 0 （关于复数Z的函数）
' 对 f(Z) 求导函数得： f'(Z)
' 对任意复数Z0可以有Z1 , Z1 = Z0 - f(Z0)/f'(Z0)
' 同样，对复数Z1可以有Z2 , Z2 = Z1 - f(Z1)/f'(Z1)
' …… …… ……
' 则，有迭代式：Z(n+1)=Z(n)-f(Z(n))/f'(Z(n))
' 对于选定的起始点，迭代大多都会收敛于方程f(z) = 0 的某个根，
' 这就是牛顿迭代法解方程的基本方式；
' 但也可能存在许多点，使迭代根本就不会收敛，
' 甚至可能出现混沌的状态。
'dX 第一次迭代x轴的变化率
'dY 第一次迭代y轴的变化率
'dL1 第一次迭代移动距离
'dL2 第nM次迭代移动距离
'dL3 第nM次迭代距离（0，0）点的距离
'dL4 迭代得到解以后距离解的大概距离
Dim x1 As Double, x2 As Double, y1 As Double, y2 As Double
Dim m0 As Double, i As Long, k As Long
Dim SeTa1 As Double, SeTa2 As Double, P1 As Double
Dim P2 As Double, A As Double, B As Double
Dim temp1 As Double, temp2 As Double, temp3 As Double, temp4 As Double, temp5 As Double
Dim temp6 As Double, temp7 As Double, temp8 As Double, temp9 As Double, temp0 As Double
Dim R1 As Double, R2 As Double
Dim Sum As Long 'Sum值用来判断迭代是否已经到达根（方程的解），这是一种比较好的方式
Dim bx As Double, by As Double, xb As Double, yb As Double
Dim N As Long, N7 As Double
On Error GoTo aaa:
x1 = x0: y1 = y0
N = Int(SeData(0, 14))
N7 = 0.01
Sum = 0
Select Case nN
Case 1 '方程 Z^N-1=0
For i = 1 To M
P1 = Sqr(x1 * x1 + y1 * y1)
SeTa1 = ZArg(x1, y1, 0)
x2 = ((N - 1) * P1 * Cos(SeTa1) + P1 ^ (1 - N) * Cos((1 - N) * SeTa1)) / N
y2 = ((N - 1) * P1 * Sin(SeTa1) + P1 ^ (1 - N) * Sin((1 - N) * SeTa1)) / N
temp0 = Abs(Abs(x1) - Abs(x2)) ^ 2 + Abs(Abs(y1) - Abs(y2)) ^ 2
If temp0 < N7 Then
Sum = Sum + 1
If Sum > 2 Then
dL4 = (temp0) / N7
If i > nM Then
Exit For
End If
End If
End If
If i = 1 Then
Hx = Abs(x1 - x2): Hy = Abs(y1 - y2)
dL1 = (Hx ^ 2 + Hy ^ 2) ^ 0.5
End If
If i = nM Then
dL2 = ((x1 - x2) ^ 2 + (y1 - y2) ^ 2) ^ 0.5
dL3 = (x2 ^ 2 + y2 ^ 2) ^ 0.5
End If
x1 = x2: y1 = y2
Next i
MMi = (k - 1) * M / N + i
Case 2, 0 '方程 Z^3-1=0 的特解
x1 = -x1 '其实没有必要，这里做了一下水平翻转
For i = 1 To M
x2 = x1: y2 = y1: m0 = (x1 * x1 + y1 * y1) ^ 2
If (x1 * x1 * x1 - x1 * y1 * y1 * 3 - 1) ^ 2 _
+ (x1 * x1 * y1 * 3 - y1 * y1 * y1) ^ 2 < N7 Then
dL4 = Abs(Sqr(x1 * x1 + y1 * y1) - 1) / N7
If i > nM Then
Exit For
End If
End If
x1 = x2 - (x2 + y2) / 6 + (x2 * x2 - y2 * y2 - x2 * y2 * 2) / 6 / m0
y1 = y2 + (x2 - y2) / 6 + (y2 * y2 - x2 * x2 - x2 * y2 * 2) / 6 / m0
If i = 1 Then
Hx = Abs(x1 - x2): Hy = Abs(y1 - y2)
dL1 = (Hx ^ 2 + Hy ^ 2) ^ 0.5
End If
If i = nM Then
dL2 = ((x1 - x2) ^ 2 + (y1 - y2) ^ 2) ^ 0.5
dL3 = (x1 ^ 2 + y1 ^ 2) ^ 0.5
End If
Next i
If x1 > 0.9 Then
MMi = i + 1 * 0.33 * M
ElseIf y1 > 0.8 Then
MMi = i + (0.33 + 1 * 0.33) * M
ElseIf y1 < -0.8 Then
MMi = i + (0.66 + 0.33 * 1) * M
End If
Case 3 '方程 Z*(1+Z^A)/(1-Z^A)=R
bx = x0: by = y0
A = SeData(0, 15) ': B = SeData(0, 16)
R1 = SeData(0, 17): R2 = SeData(0, 18)
For i = 1 To M
temp1 = bx: temp2 = by
Call Zfang(temp1, temp2, A, temp3, temp4)
Call Zji(temp3, temp4, temp1, temp2, temp5, temp6)
Call Zji(temp3, temp4, R1, R2, temp7, temp8)
temp5 = temp5 + temp7 - R1 + temp1
temp6 = temp6 + temp8 - R2 + temp2
Call Zshang(temp3, temp4, temp1, temp2, temp7, temp8)
temp7 = temp7 + 1 + (A + 1) * temp3
temp8 = temp8 + (A + 1) * temp4
Call Zshang(temp5, temp6, temp7, temp8, temp3, temp4)
bx = temp1 - temp3
by = temp2 - temp4
temp0 = Abs(Abs(bx) - Abs(temp1)) ^ 2 + Abs(Abs(by) - Abs(temp2)) ^ 2
If temp0 < N7 Then
Sum = Sum + 1
If Sum > 2 Then
dL4 = (temp0) / N7
If i > nM Then
Exit For
End If
End If
End If
If i = 1 Then
Hx = Abs(bx - temp1): Hy = Abs(by - temp2)
dL1 = (Hx ^ 2 + Hy ^ 2) ^ 0.5
End If
If i = nM Then
dL2 = ((bx - temp1) ^ 2 + (by - temp2) ^ 2) ^ 0.5
dL3 = (bx ^ 2 + by ^ 2) ^ 0.5
End If
Next i
MMi = i
End Select
Exit Function
aaa:
MMi = i
End Function
'一些复变函数
'============================
'Z1^Z2 （复数的复数次方）
Public Sub ZZFang(x1 As Double, y1 As Double, x2 As Double, y2 As Double, k As Double, x As Double, y As Double)
Dim T As Double, TT As Double
Dim P As Double, Fai As Double
On Error Resume Next
P = Sqr(x1 * x1 + y1 * y1)
If P = 0 Then
x = 0: y = 0
Exit Sub
End If
Fai = ZArg(x1, y1, k)
T = P ^ x2 * e ^ (-y2 * Fai)
TT = Log(P) * y2 + x2 * Fai
x = T * Cos(TT)
y = T * Sin(TT)
End Sub
'Z1*Z2 （复数乘积）
Public Sub Zji(x1 As Double, y1 As Double, x2 As Double, y2 As Double, x As Double, y As Double)
x = x1 * x2 - y1 * y2
y = x1 * y2 + y1 * x2
End Sub
'Z1/Z2 （复数商）
Public Sub Zshang(x1 As Double, y1 As Double, x2 As Double, y2 As Double, x As Double, y As Double)
Dim T As Double
T = x2 * x2 + y2 * y2
If T = 0 Then
If x1 = 0 Then
x = 1
y = 1
Else
x = Sgn(x1) * 1E+50
y = Sgn(y1) * 1E+50
End If
Else
x = (x1 * x2 + y1 * y2) / T
y = (-x1 * y2 + y1 * x2) / T
End If
End Sub
'Z^N （复数的实数次方）
Public Sub Zfang(x1 As Double, y1 As Double, N As Double, x As Double, y As Double)
Dim T As Double, TT As Double, AtnYX As Double
T = Sqr(x1 * x1 + y1 * y1)
AtnYX = Atn(y1 / x1)
If x1 < 0 Then
TT = PI + AtnYX
ElseIf y1 > 0 Then
TT = AtnYX
Else 'if y1<=0 then
TT = PI * 2# + AtnYX
End If
T = T ^ N
TT = TT * N
x = T * Cos(TT)
y = T * Sin(TT)
End Sub
'Arg(Z) （复数的辐角）
Public Function ZArg(x As Double, y As Double, k As Double) As Double
If x = 0 Then
If y = 0 Then
ZArg = 0
ElseIf y > 0 Then
ZArg = PI / 2#
Else 'if y<0 then
ZArg = PI * 1.5
End If
ElseIf x > 0 Then
ZArg = Atn(y / x)
If ZArg < 0 Then ZArg = PI * 2 + ZArg
Else 'if x<0 then
ZArg = Atn(y / x) + PI
End If
ZArg = ZArg + PI * 2 * k
End Function
'为了实现 Mandelbrot 特效定义的函数 （其实就是Mandelbrot函数迭代）
Public Sub fz2(x0 As Double, y0 As Double, xx As Double, yy As Double, x As Double, y As Double, N As Long)
Dim t1 As Double, t2 As Double, t3 As Double, t4 As Double
Dim i As Long
t1 = x0: t2 = y0
For i = 1 To N
t3 = t1 * t1 - t2 * t2 + xx
t4 = 2 * t1 * t2 + yy
t1 = t3: t2 = t4
Next i
x = t1: y = t2
End Sub
'=======================================================