module2.bas

一个用C语言编写的求混沌序列的源程序,挺好的!

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


'=======================================================