815.sagews

a = ComplexField()#e.g.: a(4,2) = 4+2i
b1 = .5 
b2 = 1.5
b3 = 2.0
c1 = -.5
c2 = -1.5
c3 = -2.0
d1 = a(0,.5) #.5+i
d2 = a(0,2)#.75+i
d3 = a(0,1.3)
e1 = a(0,-.3)#.5-i
e2 = a(0,-1.3)#.75-i
e3 = a(0,-2.3)
TOL = 10^(-10)

f(x) = x^4 - 1
def Muller_plot(p_0,p_1,p_2,TOL,Max):
   h_1 = p_1 - p_0
   h_2 = p_2 - p_1
   delta_1 = (f(p_1)-f(p_0))/h_1
   delta_2 = (f(p_2)-f(p_1))/h_2
   d = (delta_2 - delta_1)/(h_2+h_1)
   i = 3
   E = 1
   temp = []
   data = []
   temp = (p_0.real(),p_0.imag())
   data.append(temp)
   temp = (p_1.real(),p_1.imag())
   data.append(temp)
   temp = (p_2.real(),p_2.imag())
   data.append(temp)
   while (i <= Max):
       b = delta_2 + h_2*d
       D = (b^2-4*f(p_2)*d)^(1/2)
       if (abs(b-D)<abs(b+D)):
           E = delta_2 +h_2*d + (b^2-4*f(p_2)*d)^(1/2)
       else:
           E = delta_2 +h_2*d - (b^2-4*f(p_2)*d)^(1/2)
       h = -2*f(p_2)/E
       p = p_2 + h
       temp = (p.real(),p.imag())
       data.append(temp)
       if abs(h)<TOL:
           print p
           break
       p_0 = p_1
       p_1 = p_2
       p_2 = p
       h_1 = p_1 - p_0
       h_2 = p_2 - p_1
       delta_1 = (f(p_1) - f(p_0))/h_1
       delta_2 = (f(p_2) - f(p_1))/h_2
       d = (delta_2 - delta_1)/(h_2 +h_1)
       i = i + 1
   return data

val1 = Muller_plot(b1,b2,b3,TOL,300)
val2 = Muller_plot(c1,c2,c3,TOL,300)
val3 = Muller_plot(d1,d2,d3,TOL,300)
val4 = Muller_plot(e1,e2,e3,TOL,300)

1.00000000000000 + 3.69323567632287e-26*I -1.00000000000000 + 3.69323567632287e-26*I -1.76122249269478e-24 + 1.00000000000000*I -8.66966994115636e-30 - 1.00000000000000*I

#plot the points for the roots in red,green,purple, blue
plot_b = point(val1,rgbcolor=(1,0,0))
plot_c = point(val2,rgbcolor=(0,1,0))
plot_d = point(val3,rgbcolor=(.5,0,.5))
plot_e = point(val4,rgbcolor=(0,0,1))
graph = plot_b+plot_c+plot_d+plot_e
graph.show(xmin = -2, xmax=2,ymin=-2,ymax=2)

def Muller_eval(p_0,p_1,p_2,TOL,Max,root):
   h_1 = p_1 - p_0
   h_2 = p_2 - p_1
   delta_1 = (f(p_1)-f(p_0))/h_1
   delta_2 = (f(p_2)-f(p_1))/h_2
   d = (delta_2 - delta_1)/(h_2+h_1)
   i = 3
   E = 1
   temp = []
   data = []
   temp = (0,p_0,f(p_0),abs(p_0-root),abs(p_0-root)/abs(root))
   data.append(temp)
   temp = (1,p_1,f(p_1),abs(p_1-root),abs(p_1-root)/abs(root))
   data.append(temp)
   temp = (2,p_2,f(p_2),abs(p_2-root),abs(p_2-root)/abs(root))
   data.append(temp)
   while (i <= Max):
       b = delta_2 + h_2*d
       D = (b^2-4*f(p_2)*d)^(1/2)
       if (abs(b-D)<abs(b+D)):
           E = delta_2 +h_2*d + (b^2-4*f(p_2)*d)^(1/2)
       else:
           E = delta_2 +h_2*d - (b^2-4*f(p_2)*d)^(1/2)
       h = -2*f(p_2)/E
       p = p_2 + h
       temp = (i,p,f(p),abs(p-root),abs(p-root)/abs(root))
       data.append(temp)
       if abs(h)<TOL:
           print p
           break
       p_0 = p_1
       p_1 = p_2
       p_2 = p
       h_1 = p_1 - p_0
       h_2 = p_2 - p_1
       delta_1 = (f(p_1) - f(p_0))/h_1
       delta_2 = (f(p_2) - f(p_1))/h_2
       d = (delta_2 - delta_1)/(h_2 +h_1)
       i = i + 1
   return data

eval1 = Muller_eval(b1,b2,b3,TOL,300,1)
eval2 = Muller_eval(c1,c2,c3,TOL,300,-1)
eval3 = Muller_eval(d1,d2,d3,TOL,300,a(0,1))
eval4 = Muller_eval(e1,e2,e3,TOL,300,a(0,-1))

1.00000000000000 + 3.69323567632287e-26*I -1.00000000000000 + 3.69323567632287e-26*I -1.76122249269478e-24 + 1.00000000000000*I -8.66966994115636e-30 - 1.00000000000000*I

mat1 = matrix(eval1)
mat2 = matrix(eval2)
mat3 = matrix(eval3)
mat4 = matrix(eval4)

latex(mat1)

\left(\begin{array}{rrrrr} 0.000000000000000 & 0.500000000000000 & -0.937500000000000 & 0.500000000000000 & 0.500000000000000 \\ 1 & 1.50000000000000 & 4.06250000000000 & 0.500000000000000 & 0.500000000000000 \\ 2 & 2.00000000000000 & 15.0000000000000 & 1.00000000000000 & 1.00000000000000 \\ 3 & 1.17839458616267 & 0.928248287487500 & 0.178394586162665 & 0.178394586162665 \\ 4 & 1.00912729026519 - 0.185000743308827i & -0.170933942200123 - 0.734893174564513i & 0.185225760768861 & 0.185225760768861 \\ 5 & 1.01960098407226 + 0.0419358521920085i & 0.0697730954470841 + 0.177501041567696i & 0.0462905419677801 & 0.0462905419677801 \\ 6 & 0.998532750572800 + 0.000716778415954127i & -0.00585916700196543 + 0.00285451038838434i & 0.00163297035465910 & 0.00163297035465910 \\ 7 & 1.00000904811418 - 0.0000109678673520536i & 0.0000361922261520320 - 0.0000438726602759054i & 0.0000142183854376186 & 0.0000142183854376186 \\ 8 & 1.00000000105996 - 2.05290683612274 \times 10^{-10}i & 4.23983692598995 \times 10^{-9} - 8.21162737060294 \times 10^{-10}i & 1.07965635144463 \times 10^{-9} & 1.07965635144463 \times 10^{-9} \\ 9 & 1.00000000000000 - 2.50380277496256 \times 10^{-17}i & -1.00152110998503 \times 10^{-16}i & 2.50380277496256 \times 10^{-17} & 2.50380277496256 \times 10^{-17} \\ 10 & 1.00000000000000 + 3.69323567632287 \times 10^{-26}i & 1.47729427052915 \times 10^{-25}i & 3.69323567632287 \times 10^{-26} & 3.69323567632287 \times 10^{-26} \end{array}\right)

latex(mat2)

\left(\begin{array}{rrrrr} 0.000000000000000 & -0.500000000000000 & -0.937500000000000 & 0.500000000000000 & 0.500000000000000 \\ 1 & -1.50000000000000 & 4.06250000000000 & 0.500000000000000 & 0.500000000000000 \\ 2 & -2.00000000000000 & 15.0000000000000 & 1.00000000000000 & 1.00000000000000 \\ 3 & -1.17839458616267 & 0.928248287487500 & 0.178394586162665 & 0.178394586162665 \\ 4 & -1.00912729026519 - 0.185000743308827i & -0.170933942200123 + 0.734893174564513i & 0.185225760768861 & 0.185225760768861 \\ 5 & -1.01960098407226 + 0.0419358521920085i & 0.0697730954470841 - 0.177501041567696i & 0.0462905419677801 & 0.0462905419677801 \\ 6 & -0.998532750572800 + 0.000716778415954127i & -0.00585916700196543 - 0.00285451038838434i & 0.00163297035465910 & 0.00163297035465910 \\ 7 & -1.00000904811418 - 0.0000109678673520536i & 0.0000361922261520320 + 0.0000438726602759054i & 0.0000142183854376186 & 0.0000142183854376186 \\ 8 & -1.00000000105996 - 2.05290683612274 \times 10^{-10}i & 4.23983692598995 \times 10^{-9} + 8.21162737060294 \times 10^{-10}i & 1.07965635144463 \times 10^{-9} & 1.07965635144463 \times 10^{-9} \\ 9 & -1.00000000000000 - 2.50380277496256 \times 10^{-17}i & 1.00152110998503 \times 10^{-16}i & 2.50380277496256 \times 10^{-17} & 2.50380277496256 \times 10^{-17} \\ 10 & -1.00000000000000 + 3.69323567632287 \times 10^{-26}i & -1.47729427052915 \times 10^{-25}i & 3.69323567632287 \times 10^{-26} & 3.69323567632287 \times 10^{-26} \end{array}\right)

latex(mat3)

\left(\begin{array}{rrrrr} 0 & 0.500000000000000i & -0.937500000000000 & 0.500000000000000 & 0.500000000000000 \\ 1 & 2.00000000000000i & 15.0000000000000 & 1.00000000000000 & 1.00000000000000 \\ 2 & 1.30000000000000i & 1.85610000000000 & 0.300000000000000 & 0.300000000000000 \\ 3 & 4.68253832726578 \times 10^{-18} + 1.10850462100757i & 0.509906419514941 - 2.55125492342538 \times 10^{-17}i & 0.108504621007567 & 0.108504621007567 \\ 4 & 0.0972351339863748 + 0.937500304904909i & -0.277292058482261 - 0.317029964384680i & 0.115589286563797 & 0.115589286563797 \\ 5 & -0.00327811632213316 + 1.00245158317929i & 0.00977766056269713 + 0.0132089995551326i & 0.00409344680024076 & 0.00409344680024076 \\ 6 & -0.0000489497544265845 + 0.999983584219219i & -0.0000656758822733527 + 0.000195789374814230i & 0.0000516290259163795 & 0.0000516290259163795 \\ 7 & 1.46962469681775 \times 10^{-8} + 0.999999980934139i & -7.62634427831088 \times 10^{-8} - 5.87849845103508 \times 10^{-8}i & 2.40725305364647 \times 10^{-8} & 2.40725305364647 \times 10^{-8} \\ 8 & 1.69264559561835 \times 10^{-15} + 0.999999999999995i & -1.90958360235527 \times 10^{-14} - 6.77058238247328 \times 10^{-15}i & 5.06514893189400 \times 10^{-15} & 5.06514893189400 \times 10^{-15} \\ 9 & -1.76122249269478 \times 10^{-24} + 1.00000000000000i & 7.04488997077913 \times 10^{-24}i & 1.76122249269478 \times 10^{-24} & 1.76122249269478 \times 10^{-24} \end{array}\right)

latex(mat4)

\left(\begin{array}{rrrrr} 0 & -0.300000000000000i & -0.991900000000000 & 0.700000000000000 & 0.700000000000000 \\ 1 & -1.30000000000000i & 1.85610000000000 & 0.300000000000000 & 0.300000000000000 \\ 2 & -2.30000000000000i & 26.9841000000000 & 1.30000000000000 & 1.30000000000000 \\ 3 & -1.59698477644970 \times 10^{-17} - 1.14919620539212i & 0.744121489805099 - 9.69489963307747 \times 10^{-17}i & 0.149196205392119 & 0.149196205392119 \\ 4 & 0.146549270698739 - 0.985623744553619i & -0.180997268017854 + 0.548868262542281i & 0.147252726504441 & 0.147252726504441 \\ 5 & -0.0347634855006931 - 1.02201484849641i & 0.0834379158904821 - 0.148269612808408i & 0.0411479462242871 & 0.0411479462242871 \\ 6 & -0.000617650758879630 - 0.999283843784179i & -0.00286383473116769 - 0.00246529788097410i & 0.000945712527887977 & 0.000945712527887977 \\ 7 & 5.69107014374308 \times 10^{-6} - 1.00000094924505i & 3.79679125916965 \times 10^{-6} + 0.0000227643454009382i & 5.76969197945164 \times 10^{-6} & 5.76969197945164 \times 10^{-6} \\ 8 & 2.01381392153645 \times 10^{-10} - 1.00000000010195i & 4.07810674118991 \times 10^{-10} + 8.05525568860955 \times 10^{-10}i & 2.25718434617283 \times 10^{-10} & 2.25718434617283 \times 10^{-10} \\ 9 & -9.57638368886063 \times 10^{-20} - 1.00000000000000i & -3.83055347554425 \times 10^{-19}i & 9.57638368886063 \times 10^{-20} & 9.57638368886063 \times 10^{-20} \\ 10 & -8.66966994115636 \times 10^{-30} - 1.00000000000000i & -3.46786797646254 \times 10^{-29}i & 8.66966994115636 \times 10^{-30} & 8.66966994115636 \times 10^{-30} \end{array}\right)