Consider the cubic equation ax3 bx2 xx d 0 where a b c and

Consider the cubic equation ax^3 + bx^2 -xx + d= 0 where a, b, c and d are real input coefficients. Write a MATLAB function function [r1, r2, r3] = cubic (a, b, c, d) % a: Coefficient of x^3 % d: Coefficient of x^2 % c: Coefficient of x % d: Coefficient of 1 % r1: The first root of the cubic % r2: The second root of the cubic %, r^3: The third root of the cubic to find the largest real root of this equation accurate to within a relative error 10^-6

Solution

MATLAB code

function [x,nroot]=cubic(a, b, c, d)

    % Usage:       [x,nroot]=cubicfcn(a, b, c, d)
    % Solve a cubic equation where a, b, c, and d are real.
    %   a*x^3 + b*x^2 + c*x + d = 0

      %   a, b, c, d ... coefficients (input)
    %   x   three complex solutions (output)
    %   nroot == number of roots (output)
    

      % Local Variables:
    %   y1, y2, y3 ... three transformed solutions

    % I) If a is 0. use the quadratic formula to avoid dividing by 0.
    % II)    Calculate p and q
    %           p = ( 3*c/a - (b/a)*2 ) / 3
    %           q = ( 2*(b/a)*3 - 9*b*c/a/a + 27*d/a ) / 27

    %   III) Calculate discriminant D
    %           D = (p/3)*3 + (q/2)*2

    % IV) Depending on the sign of D,
    %           If D<0, thre distinct real roots.
    %           If D=0, three real roots of which at least two are equal.
    %           If D>0, one real and two complex roots.
  
         V) For D>0 and D=0,
    %           Calculate u and v
    %           u = cubic_root(-q/2 + sqrt(D)) // sqrt means Square root of D
    %           v = cubic_root(-q/2 - sqrt(D))

    %      VI)   Find the 3 transformed roots
    %           y1 = u + v
    %           y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
    %           y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
    %   VI) Alternately, for D<0, a trigonometric formulation is more convenient
    %           y1 = 2 * sqrt(|p|/3) * cos(phi/3)
    %           y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3)
    %           y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3)
    %           where phi = acos(-q/2/sqrt(|p|**3/27))
    %                 pi = 3.141592654...
    %   VII) Finally, find the three roots
    %           x = y - b/a/3
    %     pi = 3.141592654

    % VIII) If a is 0 use the quadratic formula.
          if (a == 0)
            [x, nroot] = quadfcn(b, c, d);
            return
          end

    % Cubic equation with 3 roots
          nroot = 3;

    % Calculate p and q
          p = c/a - b*b/a/a/3. ;
          q = (2*b*b*b/a/a/a - 9.*b*c/a/a + 27.*d/a) / 27. ;

    % Calculate DD (discriminant)
          DD = p*p*p/27. + q*q/4. ;

    % Branch to different algorithms based on DD
          if (DD < 0)
    %       3 real unequal roots -- use the trigonometric formulation
            phi = acos(-q/2./sqrt(abs(p*p*p)/27.));
            temp1=2.*sqrt(abs(p)/3.);
            y1 = temp1*cos(phi/3.);
            y2 = -temp1*cos((phi+pi)/3.);
            y3 = -temp1*cos((phi-pi)/3.);
          else
    %       1 real root & 2 conjugate complex roots OR 3 real roots
            temp1 = -q/2 + sqrt(DD);
            temp2 = -q/2 - sqrt(DD);
            u = abs(temp1)^(1./3);
            v = abs(temp2)^(1./3);
            if (temp1 < 0) u=-u;

          end
            if (temp2 < 0) v=-v;

   end
            y1 = u + v;
            y2r = -(u+v)/2;
            y2i = (u-v)*sqrt(3)/2;
end

    Final transformation
          temp1 = b/a/3;
          y1 = y1-temp1;
          if (DD < 0)
            y2 = y2-temp1;
            y3 = y3-temp1;
          else
            y2r=y2r-temp1;
          end

    % Assign answers
          if (DD < 0)
            x(1) = y1;
            x(2) = y2;
            x(3) = y3;
          elseif (DD == 0)
            x(1) = y1;
            x(2) = y2r;
            x(3) = y2r;
          else
            x(1) = y1;
            x(2) = y2r + y2i*i;
            x(3) = y2r - y2i*i;
          end

-----------------------------------------------------------------------------------------------------------------

According to above matlab code the largest root of the equation.....,

    x(2) = y2r + y2i*i;

  

Relative Error F = (error of F/ Value of F)

10^-6   = (error of F/ Value of F)