[Master Index] [Index for PublicToolbox/regutools]

# phillips

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

### Function Synopsis

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

### Help Text

    PHILLIPS Phillips' "famous" test problem.

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

     Discretization of the `famous' first-kind Fredholm integral
     equation deviced by D. L. Phillips.  Define the function
        phi(x) = | 1 + cos(x*pi/3) ,  |x| <  3 .
                 | 0               ,  |x| >= 3
     Then the kernel K, the solution f, and the right-hand side
     g are given by:
        K(s,t) = phi(s-t) ,
        f(t)   = phi(t) ,
        g(s)   = (6-|s|)*(1+.5*cos(s*pi/3)) + 9/(2*pi)*sin(|s|*pi/3) .
     Both integration intervals are [-6,6].

     The order n must be a multiple of 4.

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

    function [A,b,x] = phillips(n)
    %PHILLIPS Phillips' "famous" test problem.
    %
    % [A,b,x] = phillips(n)
    %
    % Discretization of the `famous' first-kind Fredholm integral
    % equation deviced by D. L. Phillips.  Define the function
    %    phi(x) = | 1 + cos(x*pi/3) ,  |x| <  3 .
    %             | 0               ,  |x| >= 3
    % Then the kernel K, the solution f, and the right-hand side
    % g are given by:
    %    K(s,t) = phi(s-t) ,
    %    f(t)   = phi(t) ,
    %    g(s)   = (6-|s|)*(1+.5*cos(s*pi/3)) + 9/(2*pi)*sin(|s|*pi/3) .
    % Both integration intervals are [-6,6].
    %
    % The order n must be a multiple of 4.

    % Reference: D. L. Phillips, "A technique for the numerical solution
    % of certain integral equations of the first kind", J. ACM 9
    % (1962), 84-97.

    % Discretized by Galerkin method with orthonormal box functions.

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

    % Check input.
    if (rem(n,4)~=0), error('The order n must be a multiple of 4'), end

    % Compute the matrix A.
    h = 12/n; n4 = n/4; r1 = zeros(1,n);
    c = cos([-1:n4]*4*pi/n);
    r1(1:n4) = h + 9/(h*pi^2)*(2*c(2:n4+1) - c(1:n4) - c(3:n4+2));
    r1(n4+1) = h/2 + 9/(h*pi^2)*(cos(4*pi/n)-1);
    A = toeplitz(r1);

    % Compute the right-hand side b.
    if (nargout>1),
      b = zeros(n,1); c = pi/3;
      for i=n/2+1:n
        t1 = -6 + i*h; t2 = t1 - h;
        b(i) =   t1*(6-abs(t1)/2) ...
               + ((3-abs(t1)/2)*sin(c*t1) - 2/c*(cos(c*t1) - 1))/c ...
               - t2*(6-abs(t2)/2) ...
               - ((3-abs(t2)/2)*sin(c*t2) - 2/c*(cos(c*t2) - 1))/c;
        b(n-i+1) = b(i);
      end
      b = b/sqrt(h);
    end

    % Compute the solution x.
    if (nargout==3),
      x = zeros(n,1);
      x(2*n4+1:3*n4) = (h + diff(sin([0:h:(3+10*eps)]'*c))/c)/sqrt(h);
      x(n4+1:2*n4) = x(3*n4:-1:2*n4+1);
    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