- 注册时间
- 2012-11-2
- 最后登录
- 1970-1-1
- 威望
- 星
- 金币
- 枚
- 贡献
- 分
- 经验
- 点
- 鲜花
- 朵
- 魅力
- 点
- 上传
- 次
- 下载
- 次
- 积分
- 1843
- 在线时间
- 小时
|
楼主 |
发表于 2013-1-15 19:03:45
|
显示全部楼层
函数格式:Math_Matrix_PolyRoots1(A,A_n,Ret,Erro)
A:1*A_n的多项式系数矩阵
A_n:A的列数或大小
Ret:返回的实数根
Erro:误差控制
附加说明:
Math_Matrix_PolySturm:求多项式在某一点的变号数
Math_Matrix_PolyRoots1_Call:获取求解区间,气源代码在下面
Math_Matrix_PolyDiv:多项式除法
特别说明:执行本函数前,最好先处理重根。然后执行本函数后,根据函数返回值的不同以确定是否继续求解复数根。如果返回0表示多项式的根都为实数根,这是Ret里面就是全部的根,否则(不管返回1还是-1)都得再求复数根,这时A里面的多项式就只存在复数解,则可以调用【求多项式复数根贝尔斯托(Bairstow)算法Math_Matrix_PolyRoots2】函数求解复数根。
源代码:
Public Function Math_Matrix_PolyRoots1(ByRef A(,) As Double, ByVal A_n As Integer, ByRef Ret() As Double, ByVal Erro As Integer) As Int16
'本函数求解多项式,求解其实根到Ret里面
'如果A里面存在虚根,返回False并且返回求虚根多项式A;A里面全是实根,返回True
'使用本函数前,请确定多项式没有重根
'返回1表示求得有根,但还有虚根没有求解
'返回0表示求得全部根,且无虚根
'返回-1,没有求得根,即多项式只存在虚根
Dim min As Double = Math.Abs(A(0, 1))
Dim Max As Double = Math.Abs(A(0, 0))
Dim i As Integer
Dim j As Integer
A_n -= 1
For i = 2 To A_n
If Math.Abs(A(0, i)) > min Then
min = Math.Abs(A(0, i))
End If
Next
min /= Max
Max = 1 + min
min = -Max
'在区间[min,max]求解
i = Math_Matrix_PolySturm(A, min)
j = Math_Matrix_PolySturm(A, Max)
i = Math.Abs(i - j)
If i > 0 Then
Dim N As Integer = i
Dim temp(2) As Double
Dim mtemp As Double
Dim erro1 As Double = Math.Pow(0.1, Erro)
N -= 1
ReDim Ret(N)
Dim Qujian(N, 1) As Double
Dim Loopnumber As Int16
If N > 0 Then
Math_Matrix_PolyRoots1_Call(A, min, Max, N + 1, Qujian)
Else
Qujian(0, 0) = min
Qujian(0, 1) = Max
End If
For i = 0 To N
temp(0) = A(0, A_n)
For j = 0 To A_n - 1
temp(0) += A(0, j) * Math.Pow(Qujian(i, 0), A_n - j)
Next
If Math.Abs(temp(0)) <= erro1 Then
Ret(i) = Qujian(i, 0)
Else
temp(2) = A(0, A_n)
For j = 0 To A_n - 1
temp(2) += A(0, j) * Math.Pow(Qujian(i, 1), A_n - j)
Next
If Math.Abs(temp(2)) <= erro1 Then
Ret(i) = Qujian(i, 1)
Else
temp(1) = 1
Loopnumber = 1000
While temp(1) > erro1
mtemp = (Qujian(i, 0) + Qujian(i, 1)) / 2
temp(1) = A(0, A_n)
For j = 0 To A_n - 1
temp(1) += A(0, j) * Math.Pow(mtemp, A_n - j)
Next
If Math.Sign(temp(1)) = Math.Sign(temp(0)) Then
Qujian(i, 0) = mtemp
Else
Qujian(i, 1) = mtemp
End If
temp(1) = Math.Abs(temp(1))
Loopnumber -= 1
If Loopnumber = 0 Then
Exit While
End If
End While
Ret(i) = mtemp
End If
End If
Next
If N = A.Length - 2 Then
Return 0
Else
Dim mod1(0, 0) As Double
Dim Atemp(0, 1) As Double
For i = 0 To N
Atemp(0, 0) = 1
Atemp(0, 1) = -Ret(i)
Math_Matrix_PolyDiv(A, Atemp, Nothing, mod1, 11)
j = mod1.Length - 1
ReDim A(0, j)
While j >= 0
A(0, j) = mod1(0, j)
j -= 1
End While
Next
Return 1
End If
Else
Return -1
End If
End Function
Public Sub Math_Matrix_PolyRoots1_Call(ByVal A(,) As Double, ByVal min As Double, ByVal max As Double, ByVal N As Integer, ByRef Ret(,) As Double)
'返回求解区域
Dim i As Integer = Math_Matrix_PolySturm(A, min)
Dim j As Integer = Math_Matrix_PolySturm(A, max)
Dim temp As Double = (max + min) / 2
Dim k As Integer = Math_Matrix_PolySturm(A, temp)
If N > 1 Then
i -= k
i = Math.Abs(i)
j -= k
j = Math.Abs(j)
If i = 0 Then
Math_Matrix_PolyRoots1_Call(A, temp, max, N, Ret)
ElseIf j = 0 Then
Math_Matrix_PolyRoots1_Call(A, min, temp, N, Ret)
Else
If i = 1 Then
Ret(0, 0) = min
Ret(0, 1) = temp
N -= 1
If N = 1 Then
Ret(1, 0) = temp
Ret(1, 1) = max
Else
Dim ret1(N - 1, 1) As Double
Math_Matrix_PolyRoots1_Call(A, temp, max, N, ret1)
For i = 0 To N - 1
Ret(i + 1, 0) = ret1(i, 0)
Ret(i + 1, 1) = ret1(i, 1)
Next
End If
ElseIf j = 1 Then
Ret(0, 0) = temp
Ret(0, 1) = max
N -= 1
Dim ret1(N - 1, 1) As Double
Math_Matrix_PolyRoots1_Call(A, min, temp, N, ret1)
For i = 0 To N - 1
Ret(i + 1, 0) = ret1(i, 0)
Ret(i + 1, 1) = ret1(i, 1)
Next
Else
Dim ret1(i - 1, 1) As Double
N = i - 1
Math_Matrix_PolyRoots1_Call(A, min, temp, i, ret1)
For i = 0 To N
Ret(i, 0) = ret1(i, 0)
Ret(i, 1) = ret1(i, 1)
Next
ReDim ret1(j - 1, 1)
i = j - 1
Math_Matrix_PolyRoots1_Call(A, temp, max, j, ret1)
N += 1
For j = 0 To i
Ret(j + N, 0) = ret1(j, 0)
Ret(j + N, 1) = ret1(j, 1)
Next
End If
End If
End If
End Sub |
|