郭先抢 发表于 2013-1-15 19:01:33

如何求解多项式的所有实数根呢?用什么样的算法呢?

任意给一个多项式,如何求解这个多项式的所有实数根呢?其中系数都是实数!

郭先抢 发表于 2013-1-15 19:03:15

http://blog.163.com/shikang999@126/blog/static/172624896201131702358523/
这儿有个答案,但是我不知道原理,原理没有说

郭先抢 发表于 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

      '在区间求解

      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

郭先抢 发表于 2013-1-15 19:23:50

Polynomial Root-finder (Real Coefficients)
http://www.akiti.ca/PolyRootRe.html

郭先抢 发表于 2013-1-15 21:46:53

http://www.akiti.ca/rpoly_ak1_cpp.html

郭先抢 发表于 2013-1-15 21:50:23

http://www.netlib.org/toms/493
读那么长的代码,也是一种折磨呀!!!!!!!!!!!!!!!!!!!!1

郭先抢 发表于 2013-1-15 22:03:13

General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses

http://arxiv.org/abs/1203.1034
页: [1]
查看完整版本: 如何求解多项式的所有实数根呢?用什么样的算法呢?