Detecting an error in Muller method.

Hi! I am programming Muller method as I wan to find the complex roots of an ecuation.

However, it does not work properly and I cannot detect why.

If someone could help me, I would really appreciate it.

Thanks in advance.

Here is the code:

%Muller, mejora del método de la secante, nos permite sacar tanto reales
%como complejas.
f = @(x) x.^3 + 2*x.^2 + 10*x -20;
x0 = -1;
x1 = 0;
x2 = 1;
eps_x = 10e-10;
eps_f = 10e-10;
maxits = 100;
[x0,x1,x2,n,error] = muller(f,x0,x1,x2,eps_x, eps_f,maxits)
function [x0,x1,x2, n, error] = muller(f, x0, x1, x2, eps_x, eps_f, maxits)
x0 = -1;
x1 = 0;
x2 = 1;
n = 0;
error(1) = eps_x + 1;
while n < maxits && (error(n+1) > eps_x || abs(f(x0)) > eps_f)
c = f(x2);
b = ((x0 - x2)^2*(f(x1)-f(x2))-(x1-x2)^2*(f(x0)-f(x2)))/((x0-x2)*(x1-x2)*(x0-x1));
a = ((x1-x2)*(f(x0)-f(x2))-(x0-x2)*(f(x1)-f(x2)))/((x0-x2)*(x1-x2)*(x0-x1));
x3 = x2 - (2*c)/(b+sign(b)*sqrt(b*b-4*a*c));
x0 = x1;
x1 = x2;
x2 = x3;
n = n+1;
error(n+1) = abs(x2-x1);
end
fprintf ('Las raíces son %.6f %.6f %.6f', x0, x1,x2)

end

[llucia]

I have modified some things and the result is 1 0 -1 which are the initial values, so I don´t know why it does not actualize.

[VBBV]

while n<m.axits && (error(n+1) ...

Change the above line to

while n<maxits && (error(n+1) ...

[llucia]

I have written an fprintf and I just get the real root but not the complex ones.

[VBBV]

That is due to usage of abs in this line

error(n+1) = abs(x2-x1);
             %abs 

Answer 1

[Torsten]

The argument of the root in the calculation of x2 must become negative.

Or simply start with complex values for x1,x2 and/or x3, e.g.

x0 = -1;
x1 = 0;
x2 = 3*1i;

[llucia]

I just get the real ones.

[Torsten]

%Muller, mejora del método de la secante, nos permite sacar tanto reales
%como complejas.
f = @(x) x.^3 + 2*x.^2 + 10*x -20;
x0 = -1;
x1 = 0;
x2 = 3*1i;
eps_x = 10e-10;
eps_f = 10e-10;
maxits = 100;
[x0,x1,x2,n,error] = muller(f,x0,x1,x2,eps_x, eps_f,maxits)

x0 = -1.6844 + 3.4313i

x1 = -1.6844 + 3.4313i

x2 = -1.6844 + 3.4313i

n = 29

error = 1×30

1.0000 2.5739 1.8489 1.2743 0.1398 0.0584 0.0235 0.0088 0.0033 0.0013 0.0005 0.0002 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

f(x2)

ans = -4.6850e-11 - 3.3729e-11i

function [x0,x1,x2, n, error] = muller(f, x0, x1, x2, eps_x, eps_f, maxits)
n = 0;
error(1) = eps_x + 1;
while n < maxits && (error(n+1) > eps_x || abs(f(x0)) > eps_f)
    c = f(x2);
    b = ((x0 - x2)^2*(f(x1)-f(x2))-(x1-x2)^2*(f(x0)-f(x2)))/((x0-x2)*(x1-x2)*(x0-x1));
    a = ((x1-x2)*(f(x0)-f(x2))-(x0-x2)*(f(x1)-f(x2)))/((x0-x2)*(x1-x2)*(x0-x1));
    x3 = x2 - (2*c)/(b+sign(b)*sqrt(b*b-4*a*c));
    x0 = x1;
    x1 = x2;
    x2 = x3;
    n = n+1;
    error(n+1) = abs(x2-x1);
end
%fprintf ('Las raíces son %.6f %.6f %.6f', x0, x1,x2)
end

[llucia]

perfect. Thanks a lot. It works perfectly now.

[Torsten]

In numerical computations, you will never automatically get complex numbers if you don't start with complex numbers and if there are no expressions that can generate complex numbers (like the sqrt here).