deriv2.m

这是在网上下的一个东东

function [A,b,x] = deriv2(n,example) 
%DERIV2 Test problem: computation of the second derivative. 
% 
% [A,b,x] = deriv2(n,example) 
% 
% This is a mildly ill-posed problem.  It is a discretization of a 
% first kind Fredholm integral equation whose kernel K is the 
% Green's function for the second derivative: 
%    K(s,t) = | s(t-1)  ,  s <  t . 
%             | t(s-1)  ,  s >= t 
% Both integration intervals are [0,1], and as right-hand side g 
% and correspond solution f one can choose between the following: 
%    example = 1 : g(s) = (s^3 - s)/6          ,  f(t) = t 
%    example = 2 : g(s) = exp(s) + (1-e)s - 1  ,  f(t) = exp(t) 
%    example = 3 : g(s) = | (4s^3 - 3s)/24               ,  s <  0.5 
%                         | (-4s^3 + 12s^2 - 9s + 1)/24  ,  s >= 0.5 
%                  f(t) = | t    ,  t <  0.5 
%                         | 1-t  ,  t >= 0.5 
 
% References.  The first two examples are from L. M. Delves & J. L. 
% Mohamed, "Computational Methods for Integral Equations", Cambridge 
% University Press, 1985; p. 310.  The third example is from A. K. 
% Louis & P. Maass, "A mollifier method for linear operator equations 
% of the first kind", Inverse Problems 6 (1990), 427-440. 
 
% Discretized by the Galerkin method with orthonormal box functions. 
 
% Per Christian Hansen, IMM, 04/21/97. 
 
% Initialization. 
if (nargin==1), example = 1; end 
h = 1/n; sqh = sqrt(h); h32 = h*sqh; h2 = h^2; sqhi = 1/sqh; 
t = 2/3; A = zeros(n,n); 
 
% Compute the matrix A. 
for i=1:n 
  A(i,i) = h2*((i^2 - i + 0.25)*h - (i - t)); 
  for j=1:i-1 
    A(i,j) = h2*(j-0.5)*((i-0.5)*h-1); 
  end 
end 
A = A + tril(A,-1)'; 
 
% Compute the right-hand side vector b. 
if (nargout>1) 
  b = zeros(n,1); 
  if (example==1) 
    for i=1:n 
      b(i) = h32*(i-0.5)*((i^2 + (i-1)^2)*h2/2 - 1)/6; 
    end 
  elseif (example==2) 
    ee = 1 - exp(1); 
    for i=1:n 
      b(i) = sqhi*(exp(i*h) - exp((i-1)*h) + ee*(i-0.5)*h2 - h); 
    end 
  elseif (example==3) 
    if (rem(n,2)~=0), error('Order n must be even'), else 
      for i=1:n/2 
        s12 = (i*h)^2; s22 = ((i-1)*h)^2; 
        b(i) = sqhi*(s12 + s22 - 1.5)*(s12 - s22)/24; 
      end 
      for i=n/2+1:n 
        s1 = i*h; s12 = s1^2; s2 = (i-1)*h; s22 = s2^2; 
        b(i) = sqhi*(-(s12+s22)*(s12-s22) + 4*(s1^3 - s2^3) - ... 
                    4.5*(s12 - s22) + h)/24; 
      end 
    end 
  else 
    error('Illegal value of example') 
  end 
end 
 
% Compute the solution vector x. 
if (nargout==3) 
  x = zeros(n,1); 
  if (example==1) 
    for i=1:n, x(i) = h32*(i-0.5); end 
  elseif(example==2) 
    for i=1:n, x(i) = sqhi*(exp(i*h) - exp((i-1)*h)); end 
  else 
    for i=1:n/2,   x(i) = sqhi*((i*h)^2 - ((i-1)*h)^2)/2; end 
    for i=n/2+1:n, x(i) = sqhi*(h - ((i*h)^2 - ((i-1)*h)^2)/2); end 
  end 
end