📄 nlmodule.bas
字号:
' 模块名：NLModule.bas
' 函数名：NLNdhRoot
' 功能： 使用牛顿－下山法求解实系数代数方程的全部根
' 参数 n - 多项式方程的次数
' a - Double型一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数
' xr - Double型一维数组，长度为n，返回n个根的实部
' xi - Double型一维数组，长度为n，返回n个根的虚部
' 返回值：Boolean型，若为False，则表示在QR方法中迭代已超过最大迭代次数但还未满足精度要求；
' 若为True，则表示求解成功。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLNdhRoot(n As Integer, a() As Double, xr() As Double, xi() As Double) As Boolean
Dim m As Integer, i As Integer, jt As Integer, k As Integer, nIs As Integer, it As Integer
Dim t As Double, x As Double, y As Double, x1 As Double, y1 As Double, dx As Double, dy As Double
Dim p As Double, q As Double, w As Double, dd As Double, dc As Double, c As Double
Dim g As Double, u As Double, v As Double, pq As Double, g1 As Double, u1 As Double, v1 As Double
' 系数个数
m = n + 1
' 用最高幂系数约化其余系数
k = a(1)
For i = 1 To m
a(i) = a(i) / k
Next i
' 迭代求解
k = m
nIs = 0
w = 1#
jt = 1
While (jt = 1)
pq = Abs(a(k))
While (pq < 0.000000000001)
xr(k - 1) = 0#
xi(k - 1) = 0#
k = k - 1
' 求解结束
If (k = 2) Then
xr(1) = -a(2) * w / a(1)
xi(1) = 0#
NLNdhRoot = True
Exit Function
End If
pq = Abs(a(k))
Wend
q = Log(pq)
q = q / (1# * k)
q = Exp(q)
p = q
w = w * p
For i = 2 To k
a(i) = a(i) / q
q = q * p
Next i
x = 0.0001
x1 = x
y = 0.2
y1 = y
dx = 1#
g = 1E+37
l40:
u = a(1)
v = 0#
For i = 2 To k
p = u * x1
q = v * y1
pq = (u + v) * (x1 + y1)
u = p - q + a(i)
v = pq - p - q
Next i
g1 = u * u + v * v
If (g1 >= g) Then
If (nIs <> 0) Then
it = 1
Call g65(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
If (it = 0) Then GoTo l40
Else
Call g60(t, x, y, x1, y1, dx, dy, p, q, k, it)
If (t >= 0.001) Then GoTo l40
If (g > 1E-18) Then
it = 0
Call g65(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
If (it = 0) Then GoTo l40
End If
End If
Call g90(xr, xi, a, x, y, p, q, w, k)
Else
g = g1
x = x1
y = y1
nIs = 0
If (g <= 1E-22) Then
Call g90(xr, xi, a, x, y, p, q, w, k)
Else
u1 = k * a(1)
v1 = 0#
For i = 3 To k
p = u1 * x
q = v1 * y
pq = (u1 + v1) * (x + y)
u1 = p - q + (k - i + 1) * a(i - 1)
v1 = pq - p - q
Next i
p = u1 * u1 + v1 * v1
If (p <= 1E-20) Then
it = 0
Call g65(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
If (it = 0) Then GoTo l40
Call g90(xr, xi, a, x, y, p, q, w, k)
Else
dx = (u * u1 + v * v1) / p
dy = (u1 * v - v1 * u) / p
t = 1# + 4# / k
Call g60(t, x, y, x1, y1, dx, dy, p, q, k, it)
If (t >= 0.001) Then GoTo l40
If (g > 1E-18) Then
it = 0
Call g65(x, y, x1, y1, dx, dy, dd, dc, c, k, nIs, it)
If (it = 0) Then GoTo l40
End If
Call g90(xr, xi, a, x, y, p, q, w, k)
End If
End If
End If
If (k = 2) Then
jt = 0
Else
jt = 1
End If
Wend
' 迭代结束
NLNdhRoot = True
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 模块名：NLModule.bas
' 函数名：g60
' 功能： 内部函数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub g60(t As Double, x As Double, y As Double, x1 As Double, y1 As Double, dx As Double, dy As Double, p As Double, q As Double, k As Integer, it As Integer)
it = 1
While (it = 1)
t = t / 1.67
it = 0
x1 = x - t * dx
y1 = y - t * dy
If (k > 50) Then
p = Sqr(x1 * x1 + y1 * y1)
q = Exp(85# / k)
If (p >= q) Then it = 1
End If
Wend
End Sub
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 模块名：NLModule.bas
' 函数名：g90
' 功能： 内部函数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub g90(xr() As Double, xi() As Double, a() As Double, x As Double, y As Double, p As Double, q As Double, w As Double, k As Integer)
Dim i As Integer
If (Abs(y) <= 0.000001) Then
p = -x
y = 0#
q = 0#
Else
p = -2# * x
q = x * x + y * y
xr(k - 1) = x * w
xi(k - 1) = -y * w
k = k - 1
End If
For i = 2 To k
a(i) = a(i) - a(i - 1) * p
a(i + 1) = a(i + 1) - a(i - 1) * q
Next i
xr(k - 1) = x * w
xi(k - 1) = y * w
k = k - 1
If (k = 2) Then
xr(1) = -a(2) * w / a(1)
xi(1) = 0#
End If
End Sub
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 模块名：NLModule.bas
' 函数名：g65
' 功能： 内部函数
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub g65(x As Double, y As Double, x1 As Double, y1 As Double, dx As Double, dy As Double, dd As Double, dc As Double, c As Double, k As Integer, nIs As Integer, it As Integer)
If (it = 0) Then
nIs = 1
dd = Sqr(dx * dx + dy * dy)
If (dd > 1#) Then dd = 1#
dc = 6.28 / (4.5 * k)
c = 0#
End If
While (True)
c = c + dc
dx = dd * Cos(c)
dy = dd * Sin(c)
x1 = x + dx
y1 = y + dy
If (c <= 6.29) Then
it = 0
Exit Sub
End If
dd = dd / 1.67
If (dd <= 0.0000001) Then
it = 1
Exit Sub
End If
c = 0#
Wend
End Sub
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 模块名：NLModule.bas
' 函数名：NLNdhcRoot
' 功能： 使用牛顿－下山法求解复系数代数方程的全部根
' 参数 n - 多项式方程的次数
' ar - Double型一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数的实部
' ai - Double型一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数的虚部
' xr - Double型一维数组，长度为n，返回n个根的实部
' xi - Double型一维数组，长度为n，返回n个根的虚部
' 返回值：Boolean型，若为False，则表示在QR方法中迭代已超过最大迭代次数但还未满足精度要求；
' 若为True，则表示求解成功。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function NLNdhcRoot(n As Integer, ar() As Double, ai() As Double, xr() As Double, xi() As Double) As Boolean
Dim m As Integer, i As Integer, jt As Integer, k As Integer, nIs As Integer, it As Integer
Dim t As Double, x As Double, y As Double, x1 As Double, y1 As Double, dx As Double, dy As Double
Dim p As Double, q As Double, w As Double, dd As Double, dc As Double, c As Double
Dim g As Double, u As Double, v As Double, pq As Double, g1 As Double, u1 As Double, v1 As Double
' 系数个数
m = n + 1
' 用最高幂系数约化其余系数
p = Sqr(ar(1) * ar(1) + ai(1) * ai(1))
For i = 1 To m
ar(i) = ar(i) / p
ai(i) = ai(i) / p
Next i
' 迭代求解
k = m
nIs = 0
w = 1#
jt = 1
While (jt = 1)
pq = Sqr(ar(k) * ar(k) + ai(k) * ai(k))
While (pq < 0.000000000001)
xr(k - 1) = 0#
xi(k - 1) = 0#
k = k - 1
' 求解结束
If (k = 2) Then
p = ar(1) * ar(1) + ai(1) * ai(1)
xr(1) = -w * (ar(1) * ar(2) + ai(1) * ai(2)) / p
xi(1) = w * (ar(2) * ai(1) - ar(1) * ai(2)) / p
NLNdhcRoot = True
Exit Function
End If
pq = Sqr(ar(k) * ar(k) + ai(k) * ai(k))
Wend
q = Log(pq)
q = q / (1# * k)
q = Exp(q)
p = q
w = w * p
For i = 2 To k
ar(i) = ar(i) / q
ai(i) = ai(i) / q
q = q * p
Next i
x = 0.0001
x1 = x
y = 0.2
y1 = y
dx = 1#
g = 1E+37
l40:
u = ar(1)
v = ai(1)
For i = 2 To k
p = u * x1
q = v * y1
pq = (u + v) * (x1 + y1)
u = p - q + ar(i)
v = pq - p - q + ai(i)
Next i
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -