Construcción y utilización de las series de polinomios Sturm en un ordenador (página 2)
Enviado por Aladar Peter Santha
Else
cx = cx + " X^" + Mid$(Str$(gx), 2)
End If
End If
Else
If Left$(x0(i), 1) = "-" Then
ax = " – "
Else
ax = " + "
End If
If (x0(i) <> "1" And x0(i) <> "-1") Or i = gx Then
If Left$(x0(i), 1) = "-" Then
ax = ax + Mid$(x0(i), 2)
Else
ax = ax + x0(i)
End If
End If
If gx > 1 Then
If i < gx – 1 Then
ax = ax + " X^"
ax = ax + Mid$(Str$(gx – i), 2)
Else
If i = gx – 1 Then
ax = ax + " X"
End If
End If
End If
cx = cx + ax: ax = ""
End If
End If
Next i
FormatoPol = cx
End Function
Junto a estas funciones, en un módulo del programa, hay que incluir las funciones: Sumar, Restar, Multiplicar, DivisionEuclidea, MaxComDiv, MinComMult, SumarDec, ResrarDec, MultiplicarDec CalculoMCDPOL, DivisionEspecialPOL, expuestas en las monografías 3) y 4).
Ejemplo 4: Para el polinomio
, la función AlgoritmoSturm devuelve el resultado siguiente:
Polinomio:
P(X) = X^5 -2 X^4 + 4X^3 -13X^2 + 17 X -7
Polinomio derivado:
P '(X) = 5 X^4 – 8 X^3 + 12 X^2 -26X + 17
M.C.D.[ P( x ) ; P ' ( X ) ] = X – 1
P0(X) = P( X) I M.C.D.[ P( X) ; P ' ( X )
P0(X) = X^4 – X^3 + 3 X^2 – 10 X + 7
P0 ' (X) = 4X^3 – 3X^2 + 6X-10
Cota inferior de los ceros negativos = -11.01
Cota superior de los ceros positivoss = 11.01
Sucesión de polinomios Sturm:
Ro!!nomio (1) = X^4 – X^3 + 3 X^2 – 10 X + 7 |
Polinomio ( 2) = 4 X^3 – 3 X;^2 + 6 X – 10 |
Polinomio (3) = -7 X^2 + 38 X – 34 |
Polinomio (4) = -90 X + 103 |
Polinomio (5) = -1 |
Suceciones de signos:
+ | – | – | + | – |
+ | + | – | – | – |
Número total de los ceros reales distintos = 2
Ejemplo 5: Dado el polinomio
La función AlgoritmoSturm devuelve el resultado siguiente:
Polonomio:
P( X) = 3 X^9 -17 X^8 + 8 X^7 – 5 X ^6 + 9 X^5 + 23 X^4 + 12 X^3 -17 X^2 + 8 X + 6
Polinomio derivado:
P '( X) = 27 X^8 – 136 X^7 + 56 X^6 – 30 X^5 + 45 X^4 + 92 X^3 + 36 X^2 – 34 X + 8
M.C.D.[ P( X) ; P ' (X)] = 1
P0(X) = P(X) , P0"(X) = P"(X)
Cota inferior de las raíces negativas = -8.68
Cota superior de las raíces positivas = 8.67
Sucesión de polinomios Sturm:
Polinomio(1) = 3 X^9 – 17 X^8 + 8 X^7 – 5 X^6 + 9 X^5 + 23 X^4 + 12 X^3 -17 X^2 + 8 X + 6 |
P olinomio(2) = 27 X^8 – 136 X^7 + 56 X^6 – 30 X^5 + 45 X^4 + 92 X^3 + 36 X^2 – 34 X + 8 |
Polinomio(3) = 1880 X^7 – 547 X^6 – 462 X^5 – 3870 X^4 – 3508 X^3 + 2601 X^2 -1150 X -1594 |
Pollinomio(4) = -368721 X^6 + 85974 X^5 + 2449430 X^4 + 2683036 X^3 – 3342477 X^2 + 1301670 X + 1463938 |
Polinomio(5) = -76942737 X^5 – 58265407 X^4 + 136817242 X^3 – 65538156 X^2 – 38021673 X + 12984610 |
P olinomio(6) = -594615046409 X^4 – 1429075870234 X^3 + 1360436518848 X^2 – 415022669295 X – 597878770828 |
P olinomio(7) = 6840640104890716 X^3 – 5533174222235737 X^2 + 976628467414302 X + 2276743655767350 |
Polinomio ( 8) = 13187538709035480633 X^2 – 7355256216977515014 X – 5007134039647387622 |
Polinomio(9) = -3240548216208883581377 X – 2012479906499050426312 |
Polinomio (10) = -1 |
Sucesiones de signos:
– | + | – | – | + | – | – | + | + | – |
+ | + | + | – | – | – | + | + | – | – |
Número total de las ceros reales distintos = 3
Para hacer todo el cálculo correspondiente al ejemplo 5, un ordenador de 3GHz tarda 0.328 segundos en efectuarlo.
Para hallar el número de los ceros reales de un polinomio en un intervalo solo es necesario sustituir en el código anterior la función por la función que se expone a continuación
Public Function AlgoritmoSturmI(ByRef p0() As String, intervalo() As String, n As Integer) As String
Dim i As Integer, i0 As Integer, k1 As Integer, j0 As Integer, js As Integer
Dim g0 As Integer, g1 As Integer, g2 As Integer, gz As Integer, gx As Integer
Dim sw2 As Integer, gc As Integer, s1 As Integer, sr As Integer, ist As Integer
Dim si As Integer, sd As Integer, fp() As Integer, ai As String, bi As String
Dim m2 As String, m11 As String, rp As String, x(2) As String, p1() As String
Dim rr() As String, pd() As String, pi() As String, rc As String, x0() As String
Dim ps() As String, sci As String, scs As String, s() As String, pr As String
Dim mcd() As String, pd0() As String, cx1 As String, rrr As Integer, ra As String
Dim cx2 As String, c1() As String, c2() As String, q0() As String
Dim cxp0 As String, cxp0d As String, cx0 As String, cxpd As String
rc = Chr$(13) + Chr$(10)
q0() = p0(): cxp0 = FormatoPol(q0())
' Derivada del polinomio
g0 = UBound(p0())
ReDim pd0(g0 – 1)
For i = 0 To g0 – 1
x(1) = Mid$(Str$(g0 – i), 2): x(2) = p0(i)
pd0(i) = Multiplicar(x(), n)
Next i
cxp0d = FormatoPol(pd0())
mcd() = CalculoMCDPOL(q0(), pd0(), n)
cx0 = FormatoPol(mcd())
If UBound(mcd()) <> 0 Then
c1() = DivisionEspecialPOL(p0(), mcd(), n)
Else
c1() = p0()
End If
cx1 = FormatoPol(c1())
g1 = UBound(c1())
ai = intervalo(0): bi = intervalo(1)
ReDim pd(g1 – 1)
For i = 0 To g1 – 1
x(1) = Mid$(Str$(g1 – i), 2): x(2) = c1(i)
pd(i) = Multiplicar(x(), n)
Next i
cxpd = FormatoPol(pd())
c2() = pd(): g2 = g1 – 1
js = 0: ist = 1
sr = (g1 + 2) * (g1 + 3) / 2
ReDim s(sr), z(g1), fp(g1 + 2)
For i = 0 To g1
s(js) = c1(i): js = js + 1
Next i
fp(0) = -1: fp(ist) = js – 1
ist = ist + 1: gz = g2 + 1
c2() = RutinaSturm(c2(), n)
For i = 0 To gz – 1
s(js) = c2(i)
js = js + 1
Next i
fp(ist) = js – 1
ist = ist + 1
Do
gz = 0: s1 = 0
gc = g1 – g2
For i0 = 0 To gc
If c1(i0) <> "0" Then
x(1) = c1(i0): x(2) = c2(0)
x(1) = MinComMult(x(), n): x(2) = c1(i0)
rr() = DivisionEuclidea(x(), n): m11 = rr(1)
If Left$(m11, 1) = "-" Then m11 = Mid$(m11, 2)
If m11 <> "1" Then
For k1 = i0 To g1
x(1) = c1(k1): x(2) = m11
c1(k1) = Multiplicar(x(), n)
Next k1
End If
x(1) = c1(i0): x(2) = c2(0)
rr() = DivisionEuclidea(x(), n)
m2 = rr(1)
For k1 = i0 To i0 + g2
x(1) = c2(k1 – i0): x(2) = m2: rp = Multiplicar(x(), n)
x(1) = c1(k1): x(2) = rp: c1(k1) = Restar(x(), n)
Next k1
End If
Next i0
ReDim z(g1)
j0 = gc + 1: j = j0
For i = j0 To g1
If c1(i) <> "0" Then
z(i – j) = c1(i): gz = gz + 1
If s1 = 0 Then s1 = s1 + 1
Else
If s1 = 1 Then
z(i – j) = c1(i): gz = gz + 1
Else
j = j + 1
End If
End If
Next i
sw2 = 0
For i = 0 To gz – 1
If z(i) <> "0" Then
ReDim p1(gz – 1)
For j = 0 To gz – 1: p1(j) = z(j): Next j
c1() = c2(): g1 = UBound(c1())
For i0 = 0 To gz – 1
If Left$(p1(i0), 1) = "-" Then
p1(i0) = Mid$(p1(i0), 2)
Else
If p1(i0) <> "0" Then p1(i0) = "-" + p1(i0)
End If
Next i0
c2() = p1()
g2 = UBound(c2())
c2() = RutinaSturm(c2(), n)
For i0 = 0 To gz – 1
s(js) = c2(i0): js = js + 1
Next i0
fp(ist) = js – 1: ist = ist + 1
c2() = RutinaSturm(c2(), n)
sw2 = 1: Exit For
End If
Next i
If sw2 = 0 Then Exit Do
Loop
cx2 = ""
For i = 0 To ist – 2
cx2 = cx2 + "Polinomio"
cx2 = cx2 + " (" + Str$(i + 1) + ") = "
gx = fp(i + 1) – fp(i) – 1
ReDim x0(gx)
For j = 0 To gx
x0(j) = s(fp(i) + 1 + j)
Next j
cx2 = cx2 + FormatoPol(x0()) + rc
Next i
' =========== Número de los ceros reales en el intervalo ===========
' Valores de los polinomios Sturm en el punto ai.
ReDim pi(ist – 1), ps(ist – 1)
For i = 1 To ist – 2
pi(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = pi(i): x(2) = ai: pi(i) = MultiplicarDec(x(), n)
x(1) = pi(i): x(2) = s(j): pi(i) = SumarDec(x(), n)
Next j
If pi(i) <> "0" Then
If Left$(pi(i), 1) <> "-" Then pi(i) = "1" Else pi(i) = "-1"
End If
Next i
pi(i) = s(js – 1)
'Valores de los polinomios Sturm en el punto bi.
For i = 1 To ist – 2
ps(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = ps(i): x(2) = bi: ps(i) = MultiplicarDec(x(), n)
x(1) = ps(i): x(2) = s(j): ps(i) = SumarDec(x(), n)
Next j
If ps(i) <> "0" Then
If Left$(ps(i), 1) <> "-" Then ps(i) = "1" Else ps(i) = "-1"
End If
Next i
ps(i) = s(js – 1)
For i = 1 To ist – 2
If pi(i) <> "0" Then
If Val(pi(i)) * Val(pi(i + 1)) < 0 Then si = si + 1
Else
If Val(pi(i – 1)) * Val(pi(i + 1)) <= 0 Then si = si + 1
End If
Next i
For i = 1 To ist – 2
If ps(i) <> "0" Then
If Val(ps(i)) * Val(ps(i + 1)) < 0 Then sd = sd + 1
Else
If Val(ps(i – 1)) * Val(ps(i + 1)) <= 0 Then sd = sd + 1
End If
Next i
rrr = Abs(si – sd)
' Sucesiones de signos Sturm
sci = "": scs = ""
For i = 1 To ist – 1
If pi(i) <> "0" Then
If pi(i) = "1" Then sci = sci + " + " Else sci = sci + " – "
End If
Next i
For i = 1 To ist – 1
If ps(i) <> "0" Then
If ps(i) = "1" Then scs = scs + " + " Else scs = scs + " – "
End If
Next i
ra = rc + "Polonomio:" + rc
ra = ra + "P( X ) = " + cxp0 + rc + " " + rc
ra = ra + "Polinomio derivado:" + rc
ra = ra + "P'( X ) = " + cxp0d + rc + " " + rc
ra = ra + "M.C.D.[ P( X ) ; P'( X ) ) = " + cx0 + rc + " " + rc
If cx0 <> "1" Then
ra = ra + "P0(X) = P( X ) / M.C.D.[ P( X ) ; P'( X ) ] " + rc
ra = ra + rc + "P0(X)= " + cx1 + rc
ra = ra + "P0'(X) = " + cxpd + rc
End If
ra = ra + rc + "Intervalo: " + "(" + Str$(ai) + " , " + Str$(bi) + ")" + rc
ra = ra + rc + "Sucesión de polinomios Sturm:" + rc
ra = ra + rc + cx2 + rc
ra = ra + rc + "Suceciones de signos correspondientes a las extremidades del intervalo:" + rc
ra = ra + rc + " " + sci + rc
ra = ra + " " + scs + rc
ra = ra + rc + "Número de los ceros reales distintos en el intervalo = " + Str$(rrr) + rc
AlgoritmoSturmI = ra
End Function
Utilizando esta versión del código, se pueden verificar los resultados obtenidos en el ejemplo 2
Ejemplo 6: Si calcular el número de los ceros reales distintos en el intervalo y luego hallar el número de todos los ceros reales distintos.
Utilizando ela función AAlgoritmoSturmI resulta que:
Sucesióm de polinomios Sturm:
Sucesiones de signos en las extremidades del intervalo:
– | + | – | + | – | – |
+ | – | – | + | + | – |
Número de los ceros distintos en el intervalo = 1.
Para hallar el número de todos los ceros distintos hay que hallar primero una cota inferior y una superior de los ceros. Teniendo en cuenta que las cotas de los ceros reales del polinomio sirven también para el polinomio utilizando la función CotasCeros se obtienen las cotas siguientes: -11.84 y 11.85. Calculando ahora el número de los ceros reales en el intervalo (utilizando la función AlgoritmoSturmI) se obtiene que el número total de los ceros reales es 3.
Si se quiere saber únicamente el número total de los ceros reales distintos de un polinomio o el número de los ceros reales distintos en un intervalo (sin ver la serie de polinomios Sturm y otros resultados intermedios) se puden utilizar las siguientes funciones menos extensas:
Public Function AlgoritmoSturmB(ByRef p0() As String, ByVal n As Integer) As Integer
Dim i As Integer, i0 As Integer, k1 As Integer, j0 As Integer, js As Integer
Dim g0 As Integer, g1 As Integer, g2 As Integer, gz As Integer, gx As Integer
Dim sw2 As Integer, gc As Integer, s1 As Integer, sr As Integer, ist As Integer
Dim si As Integer, sd As Integer, fp() As Integer, ai As Double, bi As Double
Dim m2 As String, m11 As String, rp As String, p1() As String, x(2) As String
Dim rr() As String, pd() As String, pi() As String, rc As String, x0() As String
Dim ps() As String, s() As String, rrr As Integer, pr As String
Dim mcd() As String, pd0() As String, cotas() As Double
Dim cx2 As String, c1() As String, c2() As String, q0() As String
Dim ra As String
rc = Chr$(13) + Chr$(10): q0() = p0()
' Derivada del polinomio
g0 = UBound(p0())
ReDim pd0(g0 – 1)
For i = 0 To g0 – 1
x(1) = Mid$(Str$(g0 – i), 2): x(2) = p0(i): pd0(i) = Multiplicar(x(), n)
Next i
mcd() = CalculoMCDPOL(q0(), pd0(), n)
If UBound(mcd()) <> 0 Then
c1() = DivisionEspecialPOL(p0(), mcd(), n)
Else
c1() = p0()
End If
g1 = UBound(c1())
cotas() = CotasCeros(c1())
ai = cotas(2): bi = cotas(1)
ReDim pd(g1 – 1)
For i = 0 To g1 – 1
x(1) = Mid$(Str$(g1 – i), 2): x(2) = c1(i): pd(i) = Multiplicar(x(), n)
Next i
c2() = pd(): g2 = g1 – 1: js = 0: ist = 1
sr = (g1 + 2) * (g1 + 3) / 2
ReDim s(sr), z(g1), fp(g1 + 2)
For i = 0 To g1
s(js) = c1(i): js = js + 1
Next i
fp(0) = -1: fp(ist) = js – 1
ist = ist + 1: gz = g2 + 1
c2() = RutinaSturm(c2(), n)
For i = 0 To gz – 1
s(js) = c2(i): js = js + 1
Next i
fp(ist) = js – 1
ist = ist + 1
Do
gz = 0: s1 = 0: gc = g1 – g2
For i0 = 0 To gc
If c1(i0) <> "0" Then
x(1) = c1(i0): x(2) = c2(0): x(1) = MinComMult(x(), n): x(2) = c1(i0)
rr() = DivisionEuclidea(x(), n): m11 = rr(1)
If Left$(m11, 1) = "-" Then m11 = Mid$(m11, 2)
If m11 <> "1" Then
For k1 = i0 To g1
x(1) = c1(k1): x(2) = m11: c1(k1) = Multiplicar(x(), n)
Next k1
End If
x(1) = c1(i0): x(2) = c2(0): rr() = DivisionEuclidea(x(), n): m2 = rr(1)
For k1 = i0 To i0 + g2
x(1) = c2(k1 – i0): x(2) = m2: rp = Multiplicar(x(), n)
x(1) = c1(k1): x(2) = rp: c1(k1) = Restar(x(), n)
Next k1
End If
Next i0
ReDim z(g1)
j0 = gc + 1: j = j0
For i = j0 To g1
If c1(i) <> "0" Then
z(i – j) = c1(i): gz = gz + 1
If s1 = 0 Then s1 = s1 + 1
Else
If s1 = 1 Then
z(i – j) = c1(i): gz = gz + 1
Else
j = j + 1
End If
End If
Next i
sw2 = 0
For i = 0 To gz – 1
If z(i) <> "0" Then
ReDim p1(gz – 1)
For j = 0 To gz – 1: p1(j) = z(j): Next j
c1() = c2(): g1 = UBound(c1())
For i0 = 0 To gz – 1
If Left$(p1(i0), 1) = "-" Then
p1(i0) = Mid$(p1(i0), 2)
Else
If p1(i0) <> "0" Then p1(i0) = "-" + p1(i0)
End If
Next i0
c2() = p1()
g2 = UBound(c2())
c2() = RutinaSturm(c2(), n)
For i0 = 0 To gz – 1
s(js) = c2(i0): js = js + 1
Next i0
fp(ist) = js – 1: ist = ist + 1
c2() = RutinaSturm(c2(), n)
sw2 = 1: Exit For
End If
Next i
If sw2 = 0 Then Exit Do
Loop
cx2 = ""
' =========== Número de todos los ceros reales ===========
' Valores de los polinomios Sturm en el punto ai.
ReDim pi(ist – 1), ps(ist – 1)
For i = 1 To ist – 2
pi(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = pi(i): x(2) = ai: pi(i) = MultiplicarDec(x(), n)
x(1) = pi(i): x(2) = s(j): pi(i) = SumarDec(x(), n)
Next j
If pi(i) <> "0" Then
If Left$(pi(i), 1) <> "-" Then pi(i) = "1" Else pi(i) = "-1"
End If
Next i
pi(i) = s(js – 1)
'Valores de los polinomios Sturm en el punto bi.
For i = 1 To ist – 2
ps(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = ps(i): x(2) = bi: ps(i) = MultiplicarDec(x(), n)
x(1) = ps(i): x(2) = s(j): ps(i) = SumarDec(x(), n)
Next j
If ps(i) <> "0" Then
If Left$(ps(i), 1) <> "-" Then ps(i) = "1" Else ps(i) = "-1"
End If
Next i
ps(i) = s(js – 1)
For i = 1 To ist – 2
If pi(i) <> "0" Then
If Val(pi(i)) * Val(pi(i + 1)) < 0 Then si = si + 1
Else
If Val(pi(i – 1)) * Val(pi(i + 1)) <= 0 Then si = si + 1
End If
Next i
For i = 1 To ist – 2
If ps(i) <> "0" Then
If Val(ps(i)) * Val(ps(i + 1)) < 0 Then sd = sd + 1
Else
If Val(ps(i – 1)) * Val(ps(i + 1)) <= 0 Then sd = sd + 1
End If
Next i
rrr = Abs(si – sd)
AlgoritmoSturmB = rrr
End Function
" – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –
Public Function AlgoritmoSturmIB(ByRef p0() As String, intervalo() As String, n As Integer) As Integer
Dim i As Integer, i0 As Integer, k1 As Integer, j0 As Integer, js As Integer
Dim g0 As Integer, g1 As Integer, g2 As Integer, gz As Integer, gx As Integer
Dim sw2 As Integer, gc As Integer, s1 As Integer, sr As Integer, ist As Integer
Dim si As Integer, sd As Integer, fp() As Integer, ai As String, bi As String
Dim m2 As String, m11 As String, rp As String, x(2) As String, p1() As String
Dim rr() As String, pd() As String, pi() As String, rc As String, x0() As String
Dim ps() As String, s() As String, pr As String, sci As String, scs As String
Dim mcd() As String, pd0() As String, cx1 As String, rrr As Integer, ra As String
Dim cx2 As String, c1() As String, c2() As String, q0() As String
rc = Chr$(13) + Chr$(10)
q0() = p0()
' Derivada del polinomio
g0 = UBound(p0())
ReDim pd0(g0 – 1)
For i = 0 To g0 – 1
x(1) = Mid$(Str$(g0 – i), 2): x(2) = p0(i)
pd0(i) = Multiplicar(x(), n)
Next i
mcd() = CalculoMCDPOL(q0(), pd0(), n)
'cx0 = FormatoPol(mcd())
If UBound(mcd()) <> 0 Then
c1() = DivisionEspecialPOL(p0(), mcd(), n)
Else
c1() = p0()
End If
g1 = UBound(c1())
ai = intervalo(0): bi = intervalo(1)
ReDim pd(g1 – 1)
For i = 0 To g1 – 1
x(1) = Mid$(Str$(g1 – i), 2): x(2) = c1(i)
pd(i) = Multiplicar(x(), n)
Next i
c2() = pd(): g2 = g1 – 1
js = 0: ist = 1
sr = (g1 + 2) * (g1 + 3) / 2
ReDim s(sr), z(g1), fp(g1 + 2)
For i = 0 To g1
s(js) = c1(i): js = js + 1
Next i
fp(0) = -1: fp(ist) = js – 1
ist = ist + 1: gz = g2 + 1
c2() = RutinaSturm(c2(), n)
For i = 0 To gz – 1
s(js) = c2(i)
js = js + 1
Next i
fp(ist) = js – 1
ist = ist + 1
Do
gz = 0: s1 = 0
gc = g1 – g2
For i0 = 0 To gc
If c1(i0) <> "0" Then
x(1) = c1(i0): x(2) = c2(0)
x(1) = MinComMult(x(), n): x(2) = c1(i0)
rr() = DivisionEuclidea(x(), n): m11 = rr(1)
If Left$(m11, 1) = "-" Then m11 = Mid$(m11, 2)
If m11 <> "1" Then
For k1 = i0 To g1
x(1) = c1(k1): x(2) = m11
c1(k1) = Multiplicar(x(), n)
Next k1
End If
x(1) = c1(i0): x(2) = c2(0)
rr() = DivisionEuclidea(x(), n)
m2 = rr(1)
For k1 = i0 To i0 + g2
x(1) = c2(k1 – i0): x(2) = m2: rp = Multiplicar(x(), n)
x(1) = c1(k1): x(2) = rp: c1(k1) = Restar(x(), n)
Next k1
End If
Next i0
ReDim z(g1)
j0 = gc + 1: j = j0
For i = j0 To g1
If c1(i) <> "0" Then
z(i – j) = c1(i): gz = gz + 1
If s1 = 0 Then s1 = s1 + 1
Else
If s1 = 1 Then
z(i – j) = c1(i): gz = gz + 1
Else
j = j + 1
End If
End If
Next i
sw2 = 0
For i = 0 To gz – 1
If z(i) <> "0" Then
ReDim p1(gz – 1)
For j = 0 To gz – 1: p1(j) = z(j): Next j
c1() = c2(): g1 = UBound(c1())
For i0 = 0 To gz – 1
If Left$(p1(i0), 1) = "-" Then
p1(i0) = Mid$(p1(i0), 2)
Else
If p1(i0) <> "0" Then p1(i0) = "-" + p1(i0)
End If
Next i0
c2() = p1()
g2 = UBound(c2())
c2() = RutinaSturm(c2(), n)
For i0 = 0 To gz – 1
s(js) = c2(i0): js = js + 1
Next i0
fp(ist) = js – 1: ist = ist + 1
c2() = RutinaSturm(c2(), n)
sw2 = 1: Exit For
End If
Next i
If sw2 = 0 Then Exit Do
Loop
' =========== Número de los ceros reales en el intervalo ===========
' Valores de los polinomios Sturm en el punto ai.
ReDim pi(ist – 1), ps(ist – 1)
For i = 1 To ist – 2
pi(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = pi(i): x(2) = ai: pi(i) = MultiplicarDec(x(), n)
x(1) = pi(i): x(2) = s(j): pi(i) = SumarDec(x(), n)
Next j
If pi(i) <> "0" Then
If Left$(pi(i), 1) <> "-" Then pi(i) = "1" Else pi(i) = "-1"
End If
Next i
pi(i) = s(js – 1)
'Valores de los polinomios Sturm en el punto bi.
For i = 1 To ist – 2
ps(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = ps(i): x(2) = bi: ps(i) = MultiplicarDec(x(), n)
x(1) = ps(i): x(2) = s(j): ps(i) = SumarDec(x(), n)
Next j
If ps(i) <> "0" Then
If Left$(ps(i), 1) <> "-" Then ps(i) = "1" Else ps(i) = "-1"
End If
Next i
ps(i) = s(js – 1)
For i = 1 To ist – 2
If pi(i) <> "0" Then
If Val(pi(i)) * Val(pi(i + 1)) < 0 Then si = si + 1
Else
If Val(pi(i – 1)) * Val(pi(i + 1)) <= 0 Then si = si + 1
End If
Next i
For i = 1 To ist – 2
If ps(i) <> "0" Then
If Val(ps(i)) * Val(ps(i + 1)) < 0 Then sd = sd + 1
Else
If Val(ps(i – 1)) * Val(ps(i + 1)) <= 0 Then sd = sd + 1
End If
Next i
rrr = Abs(si – sd)
AlgoritmoSturmIB = rrr
End Function
La function AlgoritmoSturmIB es útil cuando se calcula el número de los ceros reales distintos de un polinomio en un solo intervalo. Tiene la desventaja de que para cada intervalo nuevo calcularía de nuevo la serie de polinomios Sturm. En el caso cuando los cálculos se deben efectuar para un polinomio y varios intervalos hay que separar el cálculo de la serie Sturm y el cálculo del número de los ceros distintos en un intervalo. Para esto, las variables s(), fp(), ist, js tienen que ser globales y declarados en el apartado Option Explicit del módulo 1. Con el código siguiente se puede construir un proyecto para hallar el número de los ceros reales distintos de un polinomio sin que se calcule la serie Sturm para cada intervalo:
Option Explict
Dim s() As String, fp() As Integer, ist As Integer, js as Integer
" = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
Public Function AlgoritmoSturmIB2(ByRef p0() As String, ByRef interv() As String, n As Integer) As Integer
Dim i As Integer, i0 As Integer, k1 As Integer, j0 As Integer, rrr As Integer
Dim g0 As Integer, g1 As Integer, g2 As Integer, gz As Integer, gx As Integer
Dim sw2 As Integer, gc As Integer, s1 As Integer, sr As Integer, si As Integer
Dim m2 As String, m11 As String, rp As String, x(2) As String, sd As Integer
Dim rr() As String, pd() As String, pi() As String, rc As String, x0() As String
Dim ps() As String, pr As String, sci As String, scs As String, p1() As String
Dim mcd() As String, pd0() As String, cx1 As String, ra As String, res As String
Dim cx2 As String, c1() As String, c2() As String, q0() As String
rc = Chr$(13) + Chr$(10)
q0() = p0()
' Derivada del polinomio
g0 = UBound(p0())
ReDim pd0(g0 – 1)
For i = 0 To g0 – 1
x(1) = Mid$(Str$(g0 – i), 2): x(2) = p0(i)
pd0(i) = Multiplicar(x(), n)
Next i
mcd() = CalculoMCDPOL(q0(), pd0(), n)
If UBound(mcd()) <> 0 Then
c1() = DivisionEspecialPOL(p0(), mcd(), n)
Else
c1() = p0()
End If
g1 = UBound(c1())
ReDim pd(g1 – 1)
For i = 0 To g1 – 1
x(1) = Mid$(Str$(g1 – i), 2): x(2) = c1(i)
pd(i) = Multiplicar(x(), n)
Next i
c2() = pd(): g2 = g1 – 1
js = 0: ist = 1
sr = (g1 + 2) * (g1 + 3) / 2
ReDim s(sr), z(g1), fp(g1 + 2)
For i = 0 To g1
s(js) = c1(i): js = js + 1
Next i
fp(0) = -1: fp(ist) = js – 1
ist = ist + 1: gz = g2 + 1
c2() = RutinaSturm(c2(), n)
For i = 0 To gz – 1
s(js) = c2(i)
js = js + 1
Next i
fp(ist) = js – 1
ist = ist + 1
Do
gz = 0: s1 = 0
gc = g1 – g2
For i0 = 0 To gc
If c1(i0) <> "0" Then
x(1) = c1(i0): x(2) = c2(0)
x(1) = MinComMult(x(), n): x(2) = c1(i0)
rr() = DivisionEuclidea(x(), n): m11 = rr(1)
If Left$(m11, 1) = "-" Then m11 = Mid$(m11, 2)
If m11 <> "1" Then
For k1 = i0 To g1
x(1) = c1(k1): x(2) = m11
c1(k1) = Multiplicar(x(), n)
Next k1
End If
x(1) = c1(i0): x(2) = c2(0)
rr() = DivisionEuclidea(x(), n)
m2 = rr(1)
For k1 = i0 To i0 + g2
x(1) = c2(k1 – i0): x(2) = m2: rp = Multiplicar(x(), n)
x(1) = c1(k1): x(2) = rp: c1(k1) = Restar(x(), n)
Next k1
End If
Next i0
ReDim z(g1)
j0 = gc + 1: j = j0
For i = j0 To g1
If c1(i) <> "0" Then
z(i – j) = c1(i): gz = gz + 1
If s1 = 0 Then s1 = s1 + 1
Else
If s1 = 1 Then
z(i – j) = c1(i): gz = gz + 1
Else
j = j + 1
End If
End If
Next i
sw2 = 0
For i = 0 To gz – 1
If z(i) <> "0" Then
ReDim p1(gz – 1)
For j = 0 To gz – 1: p1(j) = z(j): Next j
c1() = c2(): g1 = UBound(c1())
For i0 = 0 To gz – 1
If Left$(p1(i0), 1) = "-" Then
p1(i0) = Mid$(p1(i0), 2)
Else
If p1(i0) <> "0" Then p1(i0) = "-" + p1(i0)
End If
Next i0
c2() = p1()
g2 = UBound(c2())
c2() = RutinaSturm(c2(), n)
For i0 = 0 To gz – 1
s(js) = c2(i0): js = js + 1
Next i0
fp(ist) = js – 1: ist = ist + 1
c2() = RutinaSturm(c2(), n)
sw2 = 1: Exit For
End If
Next i
If sw2 = 0 Then Exit Do
Loop
rrr = NRCR(interv())
AlgoritmoSturmIB2 = rrr
End Function
" – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –
Public Function NRCR(ByRef interv() As String) As Integer
' =========== Número de los ceros reales en el intervalo ===========
' Valores de los polinomios Sturm en el punto ai.
Dim pi() As String, ps() As String, si As Integer, sd As Integer
Dim i As Integer, j As Integer, rrr As Integer, x(2) As String
Dim ai As String, bi As String
ReDim pi(ist – 1), ps(ist – 1)
ai = interv(0): bi = interv(1)
For i = 1 To ist – 2
pi(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = pi(i): x(2) = ai: pi(i) = MultiplicarDec(x(), n)
x(1) = pi(i): x(2) = s(j): pi(i) = SumarDec(x(), n)
Next j
If pi(i) <> "0" Then
If Left$(pi(i), 1) <> "-" Then pi(i) = "1" Else pi(i) = "-1"
End If
Next i
pi(i) = s(js – 1)
'======== Valores de los polinomios Sturm en el punto bi. ========
For i = 1 To ist – 2
ps(i) = s(fp(i – 1) + 1)
For j = fp(i – 1) + 2 To fp(i)
x(1) = ps(i): x(2) = bi: ps(i) = MultiplicarDec(x(), n)
x(1) = ps(i): x(2) = s(j): ps(i) = SumarDec(x(), n)
Next j
If ps(i) <> "0" Then
If Left$(ps(i), 1) <> "-" Then ps(i) = "1" Else ps(i) = "-1"
End If
Next i
ps(i) = s(js – 1)
For i = 1 To ist – 2
If pi(i) <> "0" Then
If Val(pi(i)) * Val(pi(i + 1)) < 0 Then
si = si + 1
End If
Else
If Val(pi(i – 1)) * Val(pi(i + 1)) <= 0 Then
si = si + 1
End If
End If
Next i
For i = 1 To ist – 2
If ps(i) <> "0" Then
If Val(ps(i)) * Val(ps(i + 1)) < 0 Then
sd = sd + 1
End If
Else
If Val(ps(i – 1)) * Val(ps(i + 1)) <= 0 Then
sd = sd + 1
End If
End If
Next i
NRCR = Abs(si – sd)
End Function
En algunos casos podría ocurrir que un polinomio tenga en un intervalo dos ceros simples tan cercanos que sería díficil de separarles en dos intervalos disjuntos y luego calcularles por uno de los métodos conocidos (método de la bipartición, método de Newton o de la cuerda). En este caso sea el intervalo que contiene los dos ceros simples muy cercanos del polinomio y sea el punto medio del intervalo A continuación dos casos son posibles:
– Si y según el teorema de Sturm en el intervalo el número de los ceros es
los dos ceros del polinomio quedarán separados en los intervalos y Si los dos ceros del polinomio pertenecerán al intervalo
– Si entonces es uno de los dos ceros del polinomio en el intervalo y se considera Si entonces entonces y son los dos ceros del polinomio en el intervalo y así estos dos ceros quedan separados en los intervalos y donde Sise considera el intervaloPuesto queel número de ceros del polinomio en el intervalo solo puede ser ó Con el teorema de Sturm se puede establecer cuál es el caso. Si los dos ceros del intervalo quedan separados en los intervalos y Si los dos ceros del polinomio se situarán en el intervaloLa longitud del intervaloes y el intervalotendrá la longitud
Por tanto, donde
Luego, con el intervalo se procederá de la misma manera que con el intervalo y es posible llegar a e.t.c. Se obtendrá finalmente una suceción de intervalos de longitudes, tales que
(11)
(12)
(13)
Puesto que los intervalos contienen los dos ceros distintos del polinomio resulta que la sucesión (11) tiene que ser finita y así en elgun momento el número de los ceros distintos del intervalo tendrá que ser lo que significa que los dos ceros quedarán separados por uno de los números
La separación de dos ceros del polinomio, pertenecientes a un intervalo se puede realizar por el código siguiente que combina el método de la bipartición con el teorema de Sturm:
Public Function SepCeros(ByRef p() As String, interv() As String) As Variant
Dim numce As Integer, a(2) As String, x(2) As String, su As String
Dim b(2) As String, pm As String, nr As Integer, n As Integer, vp As String
Dim rc As String, res As String, rr(2, 2) As String, b1 As String
Dim aa As String, bb As String, a0 As String, a1 As String
a0 = interv(0): a1 = interv(1)
rc = Chr$(13) + Chr$(10): n = 7
numce = AlgoritmoSturmIB(px(), interv(), n)
If numce <> 2 Then
MsgBox "El número de ceros tiene que ser 2"
Exit Function
End If
a(0) = interv(0): a(1) = interv(1)
Do
x(1) = a(0): x(2) = a(1): x(1) = SumarDec(x(), n): x(2) = "0.5"
aa = MultiplicarDec(x(), n)
vp = ValPolR(p(), aa)
If vp <> "0" Then
b(0) = a(0): b(1) = aa
nr = NRCR(b())
Else
x(1) = a(0): x(2) = aa: x(1) = SumarDec(x(), n)
x(2) = "0.5": bb = MultiplicarDec(x(), n)
vp = ValPolR(p(), bb)
If vp = "0" Then
x(1) = aa: x(2) = bb: x(1) = SumarDec(x(), n)
x(2) = "0.5": aa = MultiplicarDec(x(), n)
rr(1, 1) = a0: rr(2, 1) = aa: rr(1, 2) = aa: rr(2, 2) = a1
Exit Do
Else
b(0) = a(0): b(1) = aa
nr = NRCR(b())
End If
End If
If nr = "1" Then
rr(1, 1) = a0: rr(2, 1) = aa: rr(1, 2) = aa: rr(2, 2) = a1
Exit Do
Else
If nr = "0" Then a(0) = aa
If nr = "2" Then a(1) = aa
End If
Loop
res = "Los ceros se encuentran separados en los intervalos:"
res = res + "(" + rr(1, 1) + " , " + rr(2, 1) + " )" + " y "
res = res + "(" + rr(1, 2) + " , " + rr(2, 2) + " )"
SepCeros = res
End Function
" – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –
Public Function ValPolR(ByRef p() As String, ByVal a As String) As String
Dim i As Integer, g As Integer, q() As String, x(2) As String
g = UBound(p())
ReDim q(g)
q(0) = p(0)
For i = 0 To g – 1
x(1) = q(i): x(2) = a
x(1) = Multiplicar(x(), 7): x(2) = p(i + 1)
q(i + 1) = Sumar(x(), 7)
Next i
ValPolR = q(g)
End Function
Ejemplo 7: Si se considera el polinomio
, la función AlgoritmoSturmB indica que tiene 3 ceros reales distintos. Luego según la función AlgoritmoSturmIB, en el intervalo hay dos ceros simples. Para su cálculo, los ceros de este intervalo difícilmente podrían ser separados por el método gráfico (construyendo la gráfica del polinomio) o por el método del barrido. Sin embargo, utilizando la función SepCeros (que se basa en el teorema de Sturm) se obtiene enseguida que dichos ceros se pueden separar en los intervalos:
(1, 1.732050418853759765625) y (1.732050411885375976562, 2).
Bibliografía:
1) A. ?.??PO?, K?P? B?C?E? A????P?, HAYKA, MOCKBA, 1968
2) Prof. Dr. Gh. PIC, ALGEBRA SUPERIOARA, Editura Didactica si Pedagogica, Bucuresti, 1966
3) A. Peter Santha, Cálculos con números enteros grandes en ordenadores, Monografias.com, 2012
4) A. Peter Santha, Cálculos con números decimales largos en ordenadores, Monografías.com, 2012
5) B.Démidovitch, I. Maron: Éléments de Calcule Numérique, Éditions MIR, Moscou.
Autor :
Aladar Peter Santha
Página anterior | Volver al principio del trabajo | Página siguiente |