Support Forum / Форум поддержки

Code: Select all

' test for poly_roots() function
' L1 Lagrange point for Sun - Earth system is used

n=5 ! n1=n+1
dim eq(n1),roots(n1),diff(n1), xx(n1)

' basic data
M_sun=1.989E+30        ' kg  ( mass of Sun )
M_earth=5.972E+24      ' kg  ( mass of Earth )
alfa = M_sun / M_earth
d = 1.496E+11              ' m   ( distance Sun <-> Earth )

eq(0) = -(d^5)
eq(1) = +2*d^4
eq(2) = -(d^3)
eq(3) = 3*alfa*d^2
eq(4) = -3*alfa*d
eq(5) = alfa

poly_roots(n,eq,roots,d,.001,50)

print "L1 = ";roots(1)/1000;" km"
end


' finding the (real) roots of a n-degree polynomial
' option base 0 is assumed
' n is the degree of the polynomial
' Poly() contains the n+1 coefficients p(i) in the following order:
'    Poly(n) = p(0)+p(1).x+.....+p(n).x^n
' roots() contains the roots, returned by the function, starting with
'    roots(1). If there are less than n real roots, the series of
'    roots is closed with -999
' xo is a starting value for the proces, is not critical generally
' precision is the required absolute +/- precision
' max_iterations is the allowed max # of iterations
' 
def poly_roots(n,Poly(),roots(),xo,precision,max_iterations)
dim eq(n+1),out(n+1)
for i=0 to n ! eq(i)=Poly(i) ! next i
for i=n to 1 step -1
  x=poly_root(i,eq,xo,precision,max_iterations)
  k=n-i+1 ! roots(k)=x ! if x=-999 then return x
  poly_divide(i,eq,out,x)
  for j=0 to i-1 ! eq(j)=out(j+1) ! next j ! eq(i)=0
  next i
return 1
end def

' find one (real) root of p(n)
'
def poly_root(n,cin(),xo,eps,max_iter)
dim diff(n)
poly_diff(n,cin,diff) ! x=xo
for iter=1 to max_iter
  dx=poly_value(n,cin,x)/poly_value(n-1,diff,x) ! x-=dx
  if abs(dx)<eps then return x
  next iter
return -999
end def

def poly_diff(n,cin(),cout())
for i=0 to n-1 ! cout(i)=(i+1)*cin(i+1) ! next i
end def

def poly_value(n,cin(),x)
val=cin(n)
for i=n-1 to 0 step -1 ! val*=x ! val+=cin(i) ! next i
return val
end def

def poly_divide(n,cin(),cout(),a)
for i=0 to n ! cout(i)=cin(i) ! next i
for i=n-1 to 0 step -1 ! cout(i)+=a*cout(i+1) ! next i
return cout(0)
end def

Code: Select all

' test for poly_roots() function
' 6th degree polynomial with 4 real roots

n=6 ! n1=n+1
dim eq(n1),roots(n1),diff(n1), xx(n1)

print "          (2x^2-3x+5)(x-1)(x+3)(3x-1)(x-2) = 0"
print "equals : 6x^6-11x^5-24x^4+108x^3-192x^2+143x-30 = 0" ! print

for i=0 to n ! read eq(i) ! next i
data -30, 143, -192, 108, -24, -11, 6

poly_roots(n,eq,roots,d,.001,50)

for i=1 to n ! print "x(";i;") = ";roots(i) ! next i
end


' finding the (real) roots of a n-degree polynomial
' option base 0 is assumed
' n is the degree of the polynomial
' Poly() contains the n+1 coefficients p(i) in the following order:
'    Poly(n) = p(0)+p(1).x+.....+p(n).x^n
' roots() contains the roots, returned by the function, starting with
'    roots(1). If there are less than n real roots, the series of
'    roots is closed with -999
' xo is a starting value for the proces, is not critical generally
' precision is the required absolute +/- precision
' max_iterations is the allowed max # of iterations
' 
def poly_roots(n,Poly(),roots(),xo,precision,max_iterations)
dim eq(n+1),out(n+1)
for i=0 to n ! eq(i)=Poly(i) ! next i
for i=n to 1 step -1
  x=poly_root(i,eq,xo,precision,max_iterations)
  k=n-i+1 ! roots(k)=x ! if x=-999 then return x
  poly_divide(i,eq,out,x)
  for j=0 to i-1 ! eq(j)=out(j+1) ! next j ! eq(i)=0
  next i
return 1
end def

' find one (real) root of p(n)
'
def poly_root(n,cin(),xo,eps,max_iter)
dim diff(n)
poly_diff(n,cin,diff) ! x=xo
for iter=1 to max_iter
  dx=poly_value(n,cin,x)/poly_value(n-1,diff,x) ! x-=dx
  if abs(dx)<eps then return x
  next iter
return -999
end def

def poly_diff(n,cin(),cout())
for i=0 to n-1 ! cout(i)=(i+1)*cin(i+1) ! next i
end def

def poly_value(n,cin(),x)
val=cin(n)
for i=n-1 to 0 step -1 ! val*=x ! val+=cin(i) ! next i
return val
end def

def poly_divide(n,cin(),cout(),a)
for i=0 to n ! cout(i)=cin(i) ! next i
for i=n-1 to 0 step -1 ! cout(i)+=a*cout(i+1) ! next i
return cout(0)
end def