本文程序，使用有限元计算2D瞬态热传导方程。

本程序只使用2D三角形单元求解，没有四边形单元求解。不过，也非常容易就可做到。

主程序为:

Heat_Conduction_transient.m

大部分子程序和这篇文章相同，但是多了下面的初始温度场的子程序。

u_init.m

现将各个程序列于下：

%function Heat_Conduction_transient( )%*****************************************************************************80%%% Applies the finite element method to the heat equation.%%    The user supplies datafiles that specify the geometry of the region %    and its arrangement into triangular elements, and the location and %    type of the boundary conditions, which can be any mixture of Neumann and %    Dirichlet, and the initial condition for the solution.% %    Note that only uses triangular elements in this version.%%    The unknown state variable U(x,y,t) is assumed to satisfy%    the time dependent heat equation:%%      dU(x,y,t)/dt = Uxx(x,y,t) + Uyy(x,y,t) + F(x,y,t) in Omega x [0,T]%%    with initial conditions%%      U(x,y,0) = U_INIT(x,y,0)%%    with Dirichlet boundary conditions%%      U(x,y,t) = U_D(x,y,t) on Gamma_D%%    and Neumann boundary conditions on the outward normal derivative:%%      Un(x,y) = G(x,y,t) 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.%    %    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 the heat 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.%    %%  Local Parameters:%%    Local, real DT, the size of a single time step.%%    Local, integer NT, the number of time steps to take.%%    Local, real T, the current time.%%    Local, real T_FINAL, the final time.%%    Local, real T_START, the initial time.%%  Read the nodal coordinate data file.%  load coordinates.dat;%%  Read the triangular element data file.%  load elements3.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.%  There must be at least one Dirichlet boundary condition.%  load dirichlet.dat;%%  Determine the bound and free nodes.%  BoundNodes = unique ( dirichlet );  FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );  A = sparse ( size(coordinates,1), size(coordinates,1) );  B = sparse ( size(coordinates,1), size(coordinates,1) );  t_start = 0.0;  t_final = 1;  nt = 10;% nt means number of time intervals instead of timesteps  t = t_start;  dt = ( t_final - t_start ) / nt;  U = zeros ( size(coordinates,1), nt+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(elements3,1)    B(elements3(j,:),elements3(j,:)) = B(elements3(j,:),elements3(j,:)) ...      + det ( [1,1,1; coordinates(elements3(j,:),:)' ] )...      * [ 2, 1, 1; 1, 2, 1; 1, 1, 2 ] / 24;  end%%  Set the initial condition.%  U(:,1) = u_init ( coordinates, t );%%  Given the solution at step I-1, compute the solution at step I.%  for i = 1 : nt    t = ( ( nt - i ) * t_start   ...        + (      i ) * t_final ) ...        /   nt;    b = sparse ( size(coordinates,1), 1 );    u = sparse ( size(coordinates,1), 1 );%%  Account for the volume forces, evaluated at the new time T.%    for j = 1 : size(elements3,1)      b(elements3(j,:)) = b(elements3(j,:)) ...        + det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...        dt * f ( sum(coordinates(elements3(j,:),:))/3, t ) / 6;    end%%  Account for the Neumann conditions, evaluated at the new time T.%    for j = 1 : size(neumann,1)      b(neumann(j,:)) = b(neumann(j,:)) + ...        norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...        dt * g ( sum(coordinates(neumann(j,:),:))/2, t ) / 2;    end%%  Account for terms that involve the solution at the previous timestep.%    b = b + B * U(:,i);%%  Use the Dirichlet conditions, evaluated at the new time T, to eliminate the%  known state variables.%    u(BoundNodes) = u_d ( coordinates(BoundNodes,:), t );    b = b - ( dt * A  + B ) * u;%%  Compute the remaining unknowns by solving ( dt * A + B ) * U = b.%    u(FreeNodes) = ( dt * A(FreeNodes,FreeNodes) + B(FreeNodes,FreeNodes) ) ...      \ b(FreeNodes);    U(:,i+1) = u;    end  show ( elements3, coordinates, U, t_start, t_final );

 

u_init.m

function value = u_init ( xy, t )%*****************************************************************************80%%% U_INIT sets the initial condition for the state variable.%%  Discussion:%%    The user must supply the appropriate routine for a given problem%%    In many problems, the initial time is 0.  However, the value of%    T is passed, in case the user wishes to use this same routine to%    evaluate, for instance, the exact solution.%%%  Parameters:%%    Input, real XY(N,M), contains the M-dimensional coordinates of N points.%    N is probably the total number of points, and M is probably 2.%%    Input, real T, the initial time.  %%    Output, VALUE(N), contains the value of the solution U(X,Y,T).%  n = size ( xy, 1 );  value = zeros ( n, 1 );  middle = floor ( n / 2 );  value ( middle ) = 1.0;

 

u_d.m

function value = u_d ( u, t )%*****************************************************************************80%%% U_D evaluates the Dirichlet boundary conditions.%%  Discussion:%%    The user must supply the appropriate routine for a given problem%%%  Parameters:%%    Input, real U(N,M), contains the M-dimensional coordinates of N points.%%    Output, VALUE(N), contains the value of the Dirichlet boundary%    condition at each point.%  n = size ( u, 1 );  value =zeros(n,1);% 0.1 * ( u(:,1) + u(:,2) );value(1:13)=1;

 

stima3.m

function M = stima3 ( vertices )%*****************************************************************************80%%% STIMA3 determines the local stiffness matrix for a triangular element.%%  Discussion:%%    Although this routine is intended for 2D usage, the same formulas%    work for tetrahedral elements in 3D.  The spatial dimension intended%    is determined implicitly, from the spatial dimension of the vertices.%%  Modified:%%    23 February 2004%%  Author:%%    Jochen Alberty, Carsten Carstensen, Stefan Funken.%%  Reference:%%    Jochen Alberty, Carsten Carstensen, Stefan Funken,%    Remarks Around 50 Lines of MATLAB:%    Short Finite Element Implementation,%    Numerical Algorithms,%    Volume 20, pages 117-137, 1999.%%  Parameters:%%    Input, real VERTICES(1:(D+1),1:D), contains the D-dimensional %    coordinates of the vertices.%%    Output, real M(1:(D+1),1:(D+1)), the local stiffness matrix %    for the element.%  d = size ( vertices, 2 );  D_eta = [ ones(1,d+1); vertices' ] \ [ zeros(1,d); eye(d) ];  M = det ( [ ones(1,d+1); vertices' ] ) * D_eta * D_eta' / prod ( 1:d );  returnend

show.m

function show ( elements3, coordinates, U, t_start, t_final )%*****************************************************************************80%%% SHOW displays the solution of the finite element computation.%%%  Parameters:%%    Input, integer ELEMENTS3(N3,3), the nodes that make up each triangle.%%    Input, real COORDINATES(N,1:2), the coordinates of each node.%%    Input, real U(N,NT+1), the solution, for each time step.%  %    Input, integer NT, the number of time steps.%%    Input, real T_START, T_FINAL, the start and end times.%  umin = min ( min ( U ) );  umax = max ( max ( U ) );nt=size(U,2)-1;  for i = 0 : nt%%  Print the current time T to the command window, and in the plot title.%    t = ( ( nt - i     ) * t_start   ...        + (      i     ) * t_final ) ...        /   nt;    fprintf ( 1, 'T = %f\n', t );%    picture = trisurf ( elements3, coordinates(:,1), coordinates(:,2), ...      U(:,i+1)', 'EdgeColor', 'interp', 'FaceColor', 'interp' );%%  We want all the plots to use the same Z scale.%    axis ( [0 1 0 1 umin umax] );%%  Write some labels on the plot.%    xlabel ( 'X axis' );    ylabel ( 'Y axis' );    s = sprintf ( 'T = %8f', t );    title ( s );%    drawnowend

f.m

function value = f ( xy, t )%*****************************************************************************80%%% F evaluates the right hand side of the heat equation.%%%    This routine must be changed by the user to reflect a particular problem.%%%  Parameters:%%    Input, real XY(N,M), contains the M-dimensional coordinates of N points.%%    Output, VALUE(N), contains the value of the right hand side of Laplace's%    equation at each of the points.%  n = size ( xy, 1 );  value(1:n) = 0;

g.m

function value = g ( u, t )%*****************************************************************************80%%% G evaluates the outward normal values assigned at Neumann boundary conditions.%%    This routine must be changed by the user to reflect a particular problem.%%    For this particular problem, we want to set the value of G(X,Y,T) to 1%    if X is 1, and to 0 otherwise.%%  Parameters:%%    Input, real U(N,M), contains the M-dimensional coordinates of N points.%%    Input, real T, the current time.%%    Output, VALUE(N), contains the value of the outward normal at each point%    where a Neumann boundary condition is applied.%  n = size ( u, 1 );  value = zeros ( n, 1 );

coord3.m

% the coordinate index is from 1~Nx(from left to right) for the bottom line% and then from Nx+1~2*Nx  (from left to right) for the next line above % and then nextxmin=0;xmax=1;ymin=0;ymax=1;Nx=13;Ny=13;x=linspace(xmin,xmax,Nx);y=linspace(ymin,ymax,Ny);k=0;for i1=1:Ny    for i2=1:Nx        k=k+1;        Coord(k,:)=[x(i2),y(i1)];        endendsave coordinates.dat Coord -ascii% the element indexk=0;vertices=zeros((Nx-1)*(Ny-1)*2,3);for i1=1:Ny-1    for i2=1:Nx-1        k=k+1;        ijm1=i2+(i1-1)*Nx;        ijm2=i2+1+(i1-1)*Nx;        ijm3=i2+1+i1*Nx;        vertices(k,:)=[ijm1,ijm2,ijm3];                k=k+1;        ijm1=i2+1+i1*Nx;        ijm2=i2+i1*Nx;        ijm3=i2+(i1-1)*Nx;        vertices(k,:)=[ijm1,ijm2,ijm3];            endendsave elements3.dat vertices -ascii% The direchlet boundary condition (index of the two end nodes for each boundary line)boundary=zeros(2*(Nx+Ny-2),2);temp1=1:Nx-1;temp2=2:Nx;temp3=1:Nx-1;boundary(temp3',:)=[temp1',temp2'];temp1=Nx*(1:Ny-1);temp2=temp1+Nx;temp3=temp3+Nx-1;boundary(temp3',:)=[temp1',temp2'];temp1=Ny*Nx:-1:Ny*Nx-Nx+2;temp2=temp1-1;temp3=temp3+Nx-1;boundary(temp3',:)=[temp1',temp2'];temp1=Nx*(Ny-1)+1:-Nx:Nx+1;temp2=temp1-Nx;temp3=temp3+Nx-1;boundary(temp3',:)=[temp1',temp2'];save dirichlet.dat boundary -ascii