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