How to change Matlab program for solving equation with finite element method?

Question

I don't know is this question more related to mathematics or programming and I'm absolute newbie in Matlab. Program FEM_50 applies the finite element method to Laplace's equation -Uxx(x, y) - Uyy(x, y) = F(x, y) in Omega. How to change it to apply FEM to equation -Uxx(x, y) - Uyy(x, y) + U(x, y) = F(x, y)? At this page: http://sc.fsu.edu/~burkardt/m_src/fem_50/fem_50.html additional code files in case you need them.

function fem_50 ( )

%% FEM_50 applies the finite element method to Laplace's equation.
%
%  Discussion:
%
%    FEM_50 is a set of MATLAB routines to apply the finite
%    element method to solving Laplace's equation in an arbitrary
%    region, using about 50 lines of MATLAB code.
%    
%    FEM_50 is partly a demonstration, to show how little it
%    takes to implement the finite element method (at least using
%    every possible MATLAB shortcut.)  The user supplies datafiles
%    that specify the geometry of the region and its arrangement
%    into triangular and quadrilateral elements, and the location
%    and type of the boundary conditions, which can be any mixture
%    of Neumann and Dirichlet.
%    
%    The unknown state variable U(x,y) is assumed to satisfy
%    Laplace's equation:
%      -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega
%    with Dirichlet boundary conditions
%      U(x,y) = U_D(x,y) on Gamma_D
%    and Neumann boundary conditions on the outward normal derivative:
%      Un(x,y) = G(x,y) on Gamma_N
%    If Gamma designates the boundary of the region Omega,
%    then we presume that
%      Gamma = Gamma_D + Gamma_N
%    but the user is free to determine which boundary conditions to
%    apply.  Note, however, that the problem will generally be singular
%    unless at least one Dirichlet boundary condition is specified.
%    
%    The code uses piecewise linear basis functions for triangular elements,
%    and piecewise isoparametric bilinear basis functions for quadrilateral
%    elements.
%    
%    The user is required to supply a number of data files and MATLAB
%    functions that specify the location of nodes, the grouping of nodes
%    into elements, the location and value of boundary conditions, and 
%    the right hand side function in Laplace's equation.  Note that the
%    fact that the geometry is completely up to the user means that
%    just about any two dimensional region can be handled, with arbitrary
%    shape, including holes and islands.
%    
  clear
%
%  Read the nodal coordinate data file.
%
  load coordinates.dat;
%
%  Read the triangular element data file.
%
  load elements3.dat;
%
%  Read the quadrilateral element data file.
%
  load elements4.dat;
%
%  Read the Neumann boundary condition data file.
%  I THINK the purpose of the EVAL command is to create an empty NEUMANN array
%  if no Neumann file is found.
%
  eval ( 'load neumann.dat;', 'neumann=[];' );
%
%  Read the Dirichlet boundary condition data file.
%
  load dirichlet.dat;

  A = sparse ( size(coordinates,1), size(coordinates,1) );
  b = sparse ( size(coordinates,1), 1 );
%
%  Assembly.
%
  for j = 1 : size(elements3,1)
    A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...
      + stima3(coordinates(elements3(j,:),:));
  end

  for j = 1 : size(elements4,1)
    A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ...
      + stima4(coordinates(elements4(j,:),:));
  end
%
%  Volume Forces.
%
  for j = 1 : size(elements3,1)
    b(elements3(j,:)) = b(elements3(j,:)) ...
      + det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...
      f(sum(coordinates(elements3(j,:),:))/3)/6;
  end

  for j = 1 : size(elements4,1)
    b(elements4(j,:)) = b(elements4(j,:)) ...
      + det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...
      f(sum(coordinates(elements4(j,:),:))/4)/4;
  end
%
%  Neumann conditions.
%
  if ( ~isempty(neumann) )
    for j = 1 : size(neumann,1)
      b(neumann(j,:)) = b(neumann(j,:)) + ...
        norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...
        g(sum(coordinates(neumann(j,:),:))/2)/2;
    end
  end
%
%  Determine which nodes are associated with Dirichlet conditions.
%  Assign the corresponding entries of U, and adjust the right hand side.
%
  u = sparse ( size(coordinates,1), 1 );
  BoundNodes = unique ( dirichlet );
  u(BoundNodes) = u_d ( coordinates(BoundNodes,:) );
  b = b - A * u;
%
%  Compute the solution by solving A * U = B for the remaining unknown values of U.
%
  FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );

  u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
%
%  Graphic representation.
%
  show ( elements3, elements4, coordinates, full ( u ) );

  return
end

Answer 1

I won't do your homework,... but here are some hints:

1. Go back to formulas and work out the weak form of the new equation to solve, you'll notice a new contribution to the A matrix, which you have to consider as well in the "Assembly" section.

2. Find the local (element-wise) matrix expression of this new contribution, so you can add it as well in the "Assembly" section. You'll see it's actually the mass matrix of the element.

3. Modify the "Assembly section" for all element types (triangles and quads)

4. et voilà,...

This was the "clean" way of doing it. There is actually another way that would require almost no modification of the provided program. Just turn it into a Matlab function to allow you to repeatedly solve the original problem with different body forces (F(x,y)).

Then, repeatedly call this function where you modify the body force to:

New_body_Force = Original_body_force - Solution_from_previous_call

This should converge to the desired result.

This second option is not as elegant as the "clean" way presented first but I just love fixed point iterations.

Hope this helps.