[Master Index] [Index for PublicToolbox/regutools]

# shaw

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

### Function Synopsis

    [A,b,x] = shaw(n)

### Help Text

    SHAW Test problem: one-dimensional iimage restoration model.

     [A,b,x] = shaw(n)

     Discretization of a first kind Fredholm integral equation with
     [-pi/2,pi/2] as both integration intervals.  The kernel K and
     the solution f, which are given by
        K(s,t) = (cos(s) + cos(t))*(sin(u)/u)^2
        u = pi*(sin(s) + sin(t))
        f(t) = a1*exp(-c1*(t - t1)^2) + a2*exp(-c2*(t - t2)^2) ,
     are discretized by simple quadrature to produce A and x.
     Then the discrete right-hand b side is produced as b = A*x.

     The order n must be even.

### Cross-Reference Information

This function is called by

    regudemo                           C:\BrainStorm_2001\PublicToolbox\regutools\regudemo.m

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

    function [A,b,x] = shaw(n)
    %SHAW Test problem: one-dimensional iimage restoration model.
    %
    % [A,b,x] = shaw(n)
    %
    % Discretization of a first kind Fredholm integral equation with
    % [-pi/2,pi/2] as both integration intervals.  The kernel K and
    % the solution f, which are given by
    %    K(s,t) = (cos(s) + cos(t))*(sin(u)/u)^2
    %    u = pi*(sin(s) + sin(t))
    %    f(t) = a1*exp(-c1*(t - t1)^2) + a2*exp(-c2*(t - t2)^2) ,
    % are discretized by simple quadrature to produce A and x.
    % Then the discrete right-hand b side is produced as b = A*x.
    %
    % The order n must be even.

    % Reference: C. B. Shaw, Jr., "Improvements of the resolution of
    % an instrument by numerical solution of an integral equation",
    % J. Math. Anal. Appl. 37 (1972), 83-112.

    % Per Christian Hansen, UNI-C, 08/20/91.

    % Check input.
    if (rem(n,2)~=0), error('The order n must be even'), end

    % Initialization.
    h = pi/n; A = zeros(n,n);

    % Compute the matrix A.
    co = cos(-pi/2 + [.5:n-.5]*h);
    psi = pi*sin(-pi/2 + [.5:n-.5]*h);
    for i=1:n/2
      for j=i:n-i
        ss = psi(i) + psi(j);
        A(i,j) = ((co(i) + co(j))*sin(ss)/ss)^2;
        A(n-j+1,n-i+1) = A(i,j);
      end
      A(i,n-i+1) = (2*co(i))^2;
    end
    A = A + triu(A,1)'; A = A*h;

    % Compute the vectors x and b.
    a1 = 2; c1 = 6; t1 =  .8;
    a2 = 1; c2 = 2; t2 = -.5;
    if (nargout>1)
      x =   a1*exp(-c1*(-pi/2 + [.5:n-.5]'*h - t1).^2) ...
          + a2*exp(-c2*(-pi/2 + [.5:n-.5]'*h - t2).^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