Müller's Method

Question

I have these question and I cannot solve it.

Can somebody help me?

Use Müller’s method to determine the roots of $$ f(x)=2x^5−2x^4+6x^3−6x^2+8x−8 $$ Choose $x_2=0.8 $, $x_0=0.808$ $x_1=0.792$.

Terminate your computation when the approximate relative error falls below $E_s = 10^{-4}$

Answer 1

The essential formulas in https://en.wikipedia.org/wiki/Muller%27s_method are $$ w = f[x_{k-1},x_{k-2}] + f[x_{k-1},x_{k-3}] - f[x_{k-2},x_{k-3}]. \, $$ and $$ x_{k} = x_{k-1} - \frac{2f(x_{k-1})}{w \pm \sqrt{w^2 - 4f(x_{k-1})f[x_{k-1}, x_{k-2}, x_{k-3}]}}. $$

The most difficult part is to find striking variable names... Set

• x1 for $x_{k-1}$, x2 for $x_{k-2}$ etc. and
• f1=f(x1) for f(x_{k-1}), f2=f(x2), etc. and
• for the divided differences set
  • f12=(f1-f2)/(x1-x2) for $f[x_{k-1},x_{k-2}]$,
  • f23=(f2-f3)/(x2-x3) and
  • f13=(f1-f3)/(x1-x3) and finally
  • f123=(f12-f23)/(x1-x3) for $f[x_{k-1},x_{k-2},x_{k-3}]$

---

#include<stdio.h>
#include<math.h>
#include<complex.h>

complex double f(complex double x) {
    complex double x2 = x*x;
    return (x-1)*((x2+2)*(x2+2)-x2);
}

complex double Muller(complex double (*f)(complex double), double a, double b, double c) {    
    // initialize
    unsigned int k=0;
    complex double x1 = c, f1 = f(x1);
    complex double x2 = b, f2 = f(x2), f12 = (f1-f2)/(x1-x2);
    complex double x3 = a, f3 = f(x3), f13 = (f1-f3)/(x1-x3), f23 = (f2-f3)/(x2-x3);
    complex double f123 = (f12-f23)/(x1-x3);
    double scale = fmax(cabs(f1),fmax(cabs(f2),cabs(f3)));

    printf("x[%2u]=%-15.12f%+15.12fi,\tf(x[%2u])=%-15.12f%+15.12fi,\n",
           k,creal(x3),cimag(x3),k,creal(f3),cimag(f3));
    k++;
    printf("x[%2u]=%-15.12f%+15.12fi,\tf(x[%2u])=%-15.12f%+15.12fi,\n",
           k,creal(x2),cimag(x2),k,creal(f2),cimag(f2));
    k++;
    printf("x[%2u]=%-15.12f%+15.12fi,\tf(x[%2u])=%-15.12f%+15.12fi,\n",
           k,creal(x1),cimag(x1),k,creal(f1),cimag(f1));
    k++;
    // loop
    double err = 1;
    while ( err > 1e-4 )
    {
        // iteration formulas, variables are type complex double, controls real
        complex double w = f12 - f23 + f13;
        complex double root = csqrt(w*w-4*f1*f123);
        if( cabs(w+root) < cabs(w-root) ) root = -root;
        complex double x0 = x1 - 2*f1/(w+root);
        complex double f0 = f(x0);
        printf("x[%2u]=%-15.12f%+15.12fi,\tf(x[%2u])=%-15.12f%+15.12fi,\n",
               k,creal(x0),cimag(x0),k,creal(f0),cimag(f0));
        k++;

        // shift the index for the next step, (1,2,3) <- (0,1,2)
        x3=x2; x2=x1; x1=x0;
        f3=f2; f2=f1; f1=f0;
        f23 = f12;

        // recompute remaining divided differences
        f12 = (f1-f2)/(x1-x2); f13 = (f1-f3)/(x1-x3);
        f123 = (f12-f23)/(x1-x3);

        // finish by computing control variables
        err = fmin( cabs(f1)/scale, cabs(x2-x1)/fabs(b-a) );
    }
    return x1;
}

int main() {
    Muller(*f, 0.808, 0.792, 0.8);
    printf("\n\nRestart\n\n");
    Muller(*f, -0.808, -0.792, -0.8);
    printf("\n\nRestart\n\n");
    Muller(*f, -1.808, -1.792, -1.8);
    return 0;
}

---

Log of the computation

x[ 0]=0.808000000000 +0.000000000000i,  f(x[ 0])=-1.225886093279+0.000000000000i,
x[ 1]=0.792000000000 +0.000000000000i,  f(x[ 1])=-1.305252442145+0.000000000000i,
x[ 2]=0.800000000000 +0.000000000000i,  f(x[ 2])=-1.265920000000+0.000000000000i,
x[ 3]=1.007594143502 +0.000000000000i,  f(x[ 3])=0.061333813157 +0.000000000000i,
x[ 4]=0.999691840144 +0.000000000000i,  f(x[ 4])=-0.002464329484+0.000000000000i,
x[ 5]=0.999999517022 +0.000000000000i,  f(x[ 5])=-0.000003863821+0.000000000000i,


Restart

x[ 0]=-0.808000000000+0.000000000000i,  f(x[ 0])=-11.543760711713+0.000000000000i,
x[ 1]=-0.792000000000+0.000000000000i,  f(x[ 1])=-11.245251809247+0.000000000000i,
x[ 2]=-0.800000000000+0.000000000000i,  f(x[ 2])=-11.393280000000+0.000000000000i,
x[ 3]=-0.313139488169-0.597993095065i,  f(x[ 3])=-3.578785212716-3.102850582564i,
x[ 4]=-0.152423434021-0.896568995507i,  f(x[ 4])=-2.174496945609-2.419209170588i,
x[ 5]=-0.113762245198-1.283180404737i,  f(x[ 5])=-1.974140842175-2.072205685530i,
x[ 6]=-0.590242127121-1.562744397202i,  f(x[ 6])=-1.355281768619+5.516580476110i,
x[ 7]=-0.451736511236-1.277417823451i,  f(x[ 7])=-0.401351555157-0.782834069193i,
x[ 8]=-0.493961290160-1.317349685898i,  f(x[ 8])=-0.054707398842-0.107598547954i,
x[ 9]=-0.499942023498-1.322680432217i,  f(x[ 9])=0.000220582660 -0.003039079565i,
x[10]=-0.500000033538-1.322875803330i,  f(x[10])=-0.000000312546+0.000002246640i,


Restart

x[ 0]=-1.808000000000+0.000000000000i,  f(x[ 0])=-68.773715292193+0.000000000000i,
x[ 1]=-1.792000000000+0.000000000000i,  f(x[ 1])=-66.857255669727+0.000000000000i,
x[ 2]=-1.800000000000+0.000000000000i,  f(x[ 2])=-67.809280000000+0.000000000000i,
x[ 3]=-1.182333251904-0.563769729100i,  f(x[ 3])=-10.589662224283-18.750119777722i,
x[ 4]=-0.933263172300-0.753306592649i,  f(x[ 4])=-2.028773338306-12.084168632225i,
x[ 5]=-0.707684635587-0.951498316047i,  f(x[ 5])=0.866873722369 -6.119938934357i,
x[ 6]=-0.547056458429-1.156597078482i,  f(x[ 6])=0.791543009771 -2.225497123247i,
x[ 7]=-0.493379197649-1.278355700046i,  f(x[ 7])=0.097313082739 -0.629785563704i,
x[ 8]=-0.496959319137-1.320479889960i,  f(x[ 8])=-0.029703785932-0.049268561946i,
x[ 9]=-0.499969027836-1.322900955102i,  f(x[ 9])=-0.000567484947+0.000190261047i,

Check: $$ 2x^5−2x^4+6x^3−6x^2+8x−8=2(x-1)(x^4+3x^2+4)=2(x-1)((x^2+2)^2-x^2)\\ =2(x-1)(x^2-x+2)(x^2+x+2) $$ with roots $1$ and $\pm\frac12\pm i·\frac{\sqrt{7}}2$, independent signs for 4 roots.