[Master Index] [Index for PublicToolbox/regutools]

# heat

## (PublicToolbox/regutools/heat.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    [A,b,x] = heat(n,kappa)

### Help Text

    HEAT Test problem: inverse heat equation.

     [A,b,x] = heat(n,kappa)

     A first kind Volterra integral equation with [0,1] as
     integration interval.  The kernel is K(s,t) = k(s-t) with
        k(t) = t^(-3/2)/(2*kappa*sqrt(pi))*exp(-1/(4*kappa^2*t^2)) .
     Here, kappa controls the ill-conditioning of the matrix:
        kappa = 5 gives a well-conditioned problem
        kappa = 1 gives an ill-conditioned problem.
     The default is kappa = 1.

     An exact soltuion is constructed, and then the right-hand side
     b is produced as b = A*x.

### Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\heat.m

    function [A,b,x] = heat(n,kappa)
    %HEAT Test problem: inverse heat equation.
    %
    % [A,b,x] = heat(n,kappa)
    %
    % A first kind Volterra integral equation with [0,1] as
    % integration interval.  The kernel is K(s,t) = k(s-t) with
    %    k(t) = t^(-3/2)/(2*kappa*sqrt(pi))*exp(-1/(4*kappa^2*t^2)) .
    % Here, kappa controls the ill-conditioning of the matrix:
    %    kappa = 5 gives a well-conditioned problem
    %    kappa = 1 gives an ill-conditioned problem.
    % The default is kappa = 1.
    %
    % An exact soltuion is constructed, and then the right-hand side
    % b is produced as b = A*x.

    % Reference: A. S. Carasso, "Determining surface temperatures
    % from interior observations", SIAM J. Appl. Math. 42 (1982),
    % 558-574.  See also L. Elden, "The numerical solution of a
    % non-characteristic Cauchy problem for a parabolic equation";
    % in P. Deuflhand & E. Hairer (Eds.), "Numerical Treatment of
    % Inverse Problems in Differential and Integral Equations",
    % Birkhauser, 1983.

    % Discretization by means of simple quadrature (midpoint rule).

    % Per Christian Hansen, UNI-C, 09/18/92.

    % Set default kappa.
    if (nargin==1), kappa = 1; end

    % Initialization.
    h = 1/n; t = h/2:h:1; e = ones(1,length(t));
    c = h/(2*kappa*sqrt(pi)); d = 1/(4*kappa^2);

    % Compute the matrix A.
    k = c*t.^(-1.5).*exp(-d*e./t);
    r = zeros(1,length(t)); r(1) = k(1); A = toeplitz(k,r);

    % Compute the vectors x and b.
    if (nargout>1)
      x = zeros(n,1);
      for i=1:n/2
        ti = i*20/n;
        if (ti < 2)
          x(i) = 0.75*ti^2/4;
        elseif (ti < 3)
          x(i) = 0.75 + (ti-2)*(3-ti);
        else
          x(i) = 0.75*exp(-(ti-3)*2);
        end
      end
      x(n/2+1:n) = zeros(1,n/2);
      b = A*x;
    end

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON