fork download 
  1. % ===== Helper function for Octave =====
  2. % This replaces MATLAB's "poly2str" (not available in Octave).
  3. function s = poly2str_oct(coeffs, varname)
  4.     n = length(coeffs);
  5.     terms = {};
  6.     for i = 1:n
  7.         c = coeffs(i);
  8.         p = n - i;
  9.         if abs(c) < 1e-12, continue; end
  10.         if p == 0
  11.             terms{end+1} = sprintf('%+g', c);
  12.         elseif p == 1
  13.             terms{end+1} = sprintf('%+g*%s', c, varname);
  14.         else
  15.             terms{end+1} = sprintf('%+g*%s^%d', c, varname, p);
  16.         end
  17.     end
  18.     % Join terms, remove leading "+"
  19.     s = strjoin(terms, ' ');
  20.     if s(1) == '+', s = s(2:end); end
  21. end
  22.  
  23. % ===== Main Script =====
  24. % Polynomial: a*x^3 + b*x^2 + c*x + d
  25. a = 1; b = -9; c = 0; d = 3.8197;
  26. coeffs = [a b c d];
  27. coeffs0 = coeffs;
  28.  
  29. % Initial guesses for each root
  30. x0_list = [rand()*10, rand()*10, rand()*10];
  31.  
  32. maxiter  = 50;
  33. tol      = 5e-4;   % 0.0005 => correct to 3 dp
  34.  
  35. roots_found = zeros(1,3);
  36.  
  37. % Show original polynomial
  38. pstr = poly2str_oct(coeffs0, 'x');   % <-- use helper function
  39. fprintf('f(x) = %s = 0\n', pstr);
  40.  
  41. for root_index = 1:3
  42.     fprintf('\n=== Root %d ===\n', root_index);
  43.  
  44.     dcoeffs = polyder(coeffs);
  45.     x = x0_list(root_index);
  46.  
  47.     % Table header
  48.     fprintf('iter        x_i             f(x_i)           f''(x_i)        |x_{i+1}-x_i|\n');
  49.  
  50.     for it = 1:maxiter
  51.         fval  = polyval(coeffs,  x);
  52.         fpval = polyval(dcoeffs, x);
  53.         xnew  = x - fval / fpval;
  54.         delta = abs(xnew - x);
  55.  
  56.         % Truncate to 5 dp
  57.         xt5 = fix(x*1e5)/1e5;
  58.         f5  = fix(fval*1e5)/1e5;
  59.         fp5 = fix(fpval*1e5)/1e5;
  60.         d5  = fix(delta*1e5)/1e5;
  61.  
  62.         % Print iteration
  63.         fprintf('%3d   %12.5f   %14.5f   %14.5f   %14.5f\n', it, xt5, f5, fp5, d5);
  64.  
  65.         if delta < tol
  66.             x = xnew;
  67.             break;
  68.         end
  69.         x = xnew;
  70.     end
  71.  
  72.     roots_found(root_index) = x;
  73.  
  74.     % Deflate by (x - root) unless this was the last root
  75.     if root_index < 3
  76.         [q, ~] = deconv(coeffs, [1, -x]);
  77.         coeffs = q;
  78.  
  79.         % Show remaining polynomial
  80.         q_show = poly2str_oct(fix(coeffs*1e5)/1e5, 'x');
  81.         fprintf('Remaining polynomial after dividing by (x - %.5f): %s = 0\n', fix(x*1e5)/1e5, q_show);
  82.     end
  83. end
  84.  
  85. % Final roots
  86. fprintf('\n=== All Roots (5 dp with 3 dp in parentheses) ===\n');
  87. for k = 1:length(roots_found)
  88.     r5 = fix(roots_found(k)*1e5)/1e5;
  89.     r3 = fix(roots_found(k)*1e3)/1e3;
  90.     fprintf('Root %d = %.5f (%.3f)\n', k, r5, r3);
  91. end
  92.  
Success #stdin #stdout 0.11s 46908KB
comments (?)
stdin 
Standard input is empty
stdout 
f(x) = 1*x^3 -9*x^2 +3.8197 = 0

=== Root 1 ===
iter        x_i             f(x_i)           f'(x_i)        |x_{i+1}-x_i|
  1        9.63571         62.84389        105.09812          0.59795
  2        9.03776          6.90399         82.36364          0.08382
  3        8.95393          0.12668         79.34809          0.00159
  4        8.95234          0.00004         79.29106          0.00000
Remaining polynomial after dividing by (x - 8.95233): 1*x^2 -0.04766*x -0.42667 = 0

=== Root 2 ===
iter        x_i             f(x_i)           f'(x_i)        |x_{i+1}-x_i|
  1        1.53550          1.85793          3.02335          0.61452
  2        0.92098          0.37764          1.79430          0.21046
  3        0.71051          0.04429          1.37336          0.03225
  4        0.67826          0.00104          1.30886          0.00079
  5        0.67746          0.00000          1.30727          0.00000
Remaining polynomial after dividing by (x - 0.67746): 1*x +0.6298 = 0

=== Root 3 ===
iter        x_i             f(x_i)           f'(x_i)        |x_{i+1}-x_i|
  1        5.31836          5.94816          1.00000          5.94816
  2       -0.62980          0.00000          1.00000          0.00000

=== All Roots (5 dp with 3 dp in parentheses) ===
Root 1 = 8.95233 (8.952)
Root 2 = 0.67746 (0.677)
Root 3 = -0.62980 (-0.629)
https://ideone.com/dz2whi
language:
Octave (octave 4.4.1)
created:
9 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!