function brownian_animation ( )

%*****************************************************************************80
%
%% BROWNIAN_ANIMATION animates Brownian motion by drawing one step at a time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2018
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BROWNIAN_ANIMATION:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Simulate and plot Brownian motion in two dimensions.\n' );
  fprintf ( 1, '  Animate the process by drawing one step at a time.\n' );
%
%  Simulate.
%
  t = 1.0;
  d = 10.0;
  n = 100;
  x = brownian_2d_simulate ( n, d, t );
%
%  Display.
%
  brownian_animation_display ( n, x )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BROWNIAN_ANIMATION:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function x = brownian_2d_simulate ( n, d, t )

%*****************************************************************************80
%
%% BROWNIAN_2D_SIMULATE simulates Brownian motion in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of time steps to take. 
%
%    Input, real D, the diffusion coefficient.  
%
%    Input, real T, the total time.
%
%    Output, real X(2,N), the initial position at time 0.0, and 
%    the N-1 successive locations of the particle.
%

%
%  Set the time step.
%
  dt = t / ( n - 1 );
%
%  Make space X for the locations.
%
  x = zeros ( 2, n );
%
%  Choose stepsizes S that are normal random values.
%
  s = sqrt ( 2.0 * d * dt ) * randn ( 1, n - 1 );
%
%  The steps DX are a combination of a random stepsize and a random angle.
%
  a = randn ( 2, n - 1 );
  v = s ./ sqrt ( sum ( a.^2 ) );
  b = spdiags ( v', 0, n - 1, n - 1 );
  dx(1:2,1:n-1) = a * b;
%
%  Each position X is the sum of the previous steps.
%
  x(1:2,2:n) = cumsum ( dx(1:2,1:n-1), 2 );

  return
end
function brownian_animation_display ( n, x )

%*****************************************************************************80
%
%% BROWNIAN_ANIMATION_DISPLAY displays successive Brownian motion positions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of time steps. 
%
%    Input, real X(2,N), the particle positions.
%
  xmin = min ( x(1,:) );
  xmax = max ( x(1,:) );
  ymin = min ( x(2,:) );
  ymax = max ( x(2,:) );

  clf
  axis ( [ xmin, xmax, ymin, ymax ] )
  grid on
  xlabel ( '<--X-->' )
  ylabel ( '<--Y-->' )
  title ( 'Animation of Brownian motion' )

  hold on

    plot ( x(1,1), x(2,1), 'g.', 'MarkerSize', 35 );
    for i = 1 : n - 1
      plot ( [ x(1,i), x(1,i+1) ], [ x(2,i), x(2,i+1) ], 'linewidth', 2, 'Color', 'b' );
      pause ( 0.25 )
    end
    plot ( x(1,n), x(2,n), 'r.', 'MarkerSize', 35 );

  hold off

  filename = sprintf ( 'brownian_animation.png' );
  print ( '-dpng', filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved as "%s".\n', filename );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

