function p1_ex4 ( )

%*****************************************************************************80
%
%% P1_EX4 carries out project 1 exercise 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2019
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'P1_EX4\n' );
  fprintf ( 1, '  Create datasets S1 and S2.\n' );
  fprintf ( 1, '  Create polynomial models for S1.\n' );
  fprintf ( 1, '  Test them against the data in S2.\n' );
%
%  Create the data sets.
%
  m1 = 10;
  data1 = sinus ( m1 );
  x1 = data1(:,1);
  y1 = data1(:,2);

  m2 = 20;
  data2 = sinus ( m2 );
  x2 = data2(:,1);
  y2 = data2(:,2);
%
%  Plot the data.
%
  x = linspace ( 0.0, 2.0 * pi, 101 );
  y = sin ( x );
  clf
  hold on
  plot ( x, y, 'k-', 'linewidth', 3 );
  plot ( x1, y1, 'b.', 'markersize', 30 );
  plot ( x2, y2, 'r.', 'markersize', 30 );
  xlabel ( '<-- X -->' );
  ylabel ( '<-- Y -->' );
  title ( 'The S1 (blue) and S2 (red) Datasets', 'fontsize', 16 );
  grid ( 'on' )
  hold off
  filename = 'p1_ex4_data.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Fill in the table for n = 1 : 10
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  TABLE 4 (using SVD/pinv):\n' );
  fprintf ( 1, '   N    rms(r1)    rms(r1)\n' );
  fprintf ( 1, '\n' );
%
%  Set up the linear system A * c = b for [1,..,x^(n-1)] and solve
%  with SVD (pinv).
%
  A = [];

  for n = 1 : 10
    A = [ A, x1.^(n-1) ];
    b = y1;
    c = pinv ( A ) * b;
    r1 = A * c - b;
    r2 = y2 - polyval ( c(n:-1:1), x2 );
    fprintf ( 1, '  %2d  %8.5f  %8.5f\n', n, rms ( r1 ), rms ( r2 ) );
  end
%
%  Set up the linear system A * c = b for [1,..,x^(n-1)] and solve
%  with QR.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  TABLE 4 (using QR):\n' );
  fprintf ( 1, '   N    rms(r1)    rms(r1)\n' );
  fprintf ( 1, '\n' );

  A = [];

  for n = 1 : 10
    A = [ A, x1.^(n-1) ];
    b = y1;
    [ Q, R ] = qr ( A );
    c = R \ ( Q' * b );
    r1 = A * c - b;
    r2 = y2 - polyval ( c(n:-1:1), x2 );
    fprintf ( 1, '  %2d  %8.5f  %8.5f\n', n, rms ( r1 ), rms ( r2 ) );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'P1_EX4\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end