function biharmonic_fd1d ( n )

%*****************************************************************************80
%
%% biharmonic_fd1d() solves the 1D biharmonic equation using finite differences.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 May 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of equally spaced nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd1d():\n' );
  fprintf ( 1, '  Use finite differences to solve the biharmonic equation\n' );
  fprintf ( 1, '  over a 1D interval.\n' );

  if ( nargin < 1 )
    n = 17;
  end
%
%  Set the grid.
%
  xmin = -1.0;
  xmax = +1.0;
  x = linspace ( xmin, xmax, n );
  x = x';
  dx = ( xmax - xmin ) / ( n - 1 );
%
%  Get the basic matrix.
%
  A = biharmonic_matrix ( n );
%
%  Get the basic right hand side.
%
  b = exp ( x );
  b = b * dx^4;
%
%  Get the boundary conditions on the left and right endpoints for U and U'.
%
  ul = 0.0;
  upl = 0.0;
  ur = 0.0;
  upr = 0.0;
%
%  Modify A and b for the Dirichlet condition on U at left.
%
  A(1,1) = 1.0;
  A(1,2) = 0.0;
  A(1,3) = 0.0;
  b(1) = ul;
%
%  Modify A and b for the Dirichlet condition on U at right.
%
  A(n,n-2) = 0.0;
  A(n,n-1) = 0.0;
  A(n,n)   = 1.0;
  b(n) = ur;
%
%  Modify A and b for the Dirichlet condition on U' at left.
%
%  Dirichlet condition on u' at left is equivalent to 
%    ( u2 - u0 ) / 2dx = upl
%  so fictitious u0 satisfies
%    u0 = u2 - 2dx * upl
%  so equation at x2
%    ( 1 * u0 - 4 * u1 + 6 * u2 - 4 * u3 + 1 * u4 ) / dx^4 = f(x2)
%  becomes
%           ( - 4 * u1 + 7 * u2 - 4 * u3 + 1 * u4 ) / dx^4 = f(x2) + 2dx * upl / dx^4
%
  A(2,2) = 7.0;
  b(2) = b(2) + 2.0 * dx * upl;
%
%  Modify A and b for the Dirichlet condition on U' at right.
%
%  Dirichlet condition on u' at right is equivalent to 
%    ( un+1 - un-1 ) / 2dx = upr
%  so fictitious un+1 satisfies
%    un+1 = un-1 + 2dx * upr
%  so equation at xn-1
%    ( 1 * un-3 - 4 * un-2 + 6 * un-1 - 4 * un + 1 * un+1 ) / dx^4 = f(xn-1)
%  becomes
%    ( 1 * un-3 - 4 * un-2 + 7 * un-1 - 4 * un            ) / dx^4 = f(xn-1) - 2dx * upr / dx^4
%
  A(n-1,n-1) = 7.0;
  b(n-1) = b(n-1) - 2.0 * dx * upr;
%
%  Solve the system.
%
  u = A \ b;
%
%  Evaluate the exact solution.
%
  xe = linspace ( xmin, xmax, 51 );
  xe = xe';
  ue = biharmonic_exact ( xe );
%
%  Plot the solution.
%
  clf ( );
  hold ( 'on' );
    plot ( xe, ue, 'b-', 'linewidth', 3 );
    plot ( x, u, 'r.', 'markersize', 30 );
  hold ( 'off' );
  grid ( 'on' );
  title ( 'biharmonic 1d, Finite Differences' );
  legend ( 'Exact solution', 'Finite Difference solution' );
  filename = 'biharmonic_fd1d.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'biharmonic_fd1d():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function A = biharmonic_matrix ( n )

%*****************************************************************************80
%
%% biharmonic_matrix() returns a matrix for the 1D biharmonic equation.
%
%  Discussion:
%
%    To compute the fourth derivative, the matrix entries would normally
%    include a divisor of dx^4.  For simplicity, this divisor is omitted.
%    To solve A*x=b correctly, then, the right hand side b must be
%    multiplied by dx^4.
%
%    The first two and last two rows of A and b must be modified, to account
%    for boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 May 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of equally spaced nodes.
%
%  Output:
%
%    sparse A(N,N): the biharmonic matrix.
%
  e = ones ( n, 1 );
  B = [1,-4,6,-4,1];
  eB = e*B;

  A = spdiags ( eB, -2:+2, n, n );

  return
end
function c = biharmonic_coefficients ( )

%*****************************************************************************80
%
%% biharmonic_coefficients() returns the polynomial coefficients of the exact solution.
%
%  Discussion:
%
%    The exact solution has the form
%
%      u(x) = a + bx + cx^2 + dx^3 + exp(x)
%
%    Given u(-1) = u'(-1) = u(1) = u'(1) = 0, we can determine a, b, c and d.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 May 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of equally spaced nodes.
%
%  Output:
%
%    real C(4): the biharmonic matrix.
%
  xl = -1.0;
  xr = +1.0;
%
% ul  = 0 = a + b*xl +     c*xl^2 +     d*xl^3 + exp(xl);
% upl = 0 =     b    + 2.0*c*xl   + 3.0*d*xl^2 + exp(xl);
% ur  = 0 = a + b*xr +     c*xr^2 +     d*xr^3 + exp(xr);
% upr = 0 =     b    + 2.0*c*xr   + 3.0*d*xr^2 + exp(xr);
%
  A = [ 1.0,  xl,     xl^2,     xl^3; ...
        0.0, 1.0, 2.0*xl,   3.0*xl^2; ...
        1.0,  xr,     xr^2,     xr^3; ...
        0.0, 1.0, 2.0*xr,   3.0*xr^2 ];

  b = [ -exp(xl); -exp(xl); -exp(xr); -exp(xr) ];

  c = A \ b;
  
  return
end
function ue = biharmonic_exact ( xe )

%*****************************************************************************80
%
%% biharmonic_exact() evaluates the exact solution.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 May 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real XE(NE): evaluation points.
%
%  Output:
%
%    real UE(NE): the exact solution values.
%
  c = biharmonic_coefficients ( );

  ue = c(1) + c(2) * xe + c(3) * xe.^2 + c(4) * xe.^3 + exp(xe);
  
  return
end

