www.gusucode.com > funfun工具箱matlab源码程序 > funfun/bvpinit.m

    function solinit = bvpinit(x,v,parameters,varargin)
%BVPINIT  Form the initial guess for BVP solvers.
%   SOLINIT = BVPINIT(X,YINIT) forms the initial guess for BVP4C or BVP5C in 
%   common circumstances. The boundary value problem (BVP) is to be solved 
%   on [a,b]. The vector X specifies a and b as X(1) = a and X(end) = b. 
%   It is also a guess for an appropriate mesh. BVP solvers will adapt this 
%   mesh to the solution, so often a guess like X = linspace(a,b,10) will 
%   suffice, but in difficult cases, mesh points should be placed where the 
%   solution changes rapidly. 
%
%   The entries of X must be ordered. For two-point BVPs, the entries of X 
%   must be distinct, so if a < b, then X(1) < X(2) < ... < X(end), and 
%   similarly for a > b. For multipoint BVPs there are boundary conditions
%   at points in [a,b]. Generally, these points represent interfaces and 
%   provide a natural division of [a,b] into regions. BVPINIT enumerates 
%   the regions from left to right (from a to b), with indices starting 
%   from 1. You can specify interfaces by double entries in the initial 
%   mesh X. BVPINIT interprets one entry as the right end point of region k 
%   and the other as the left end point of region k+1. THREEBVP exemplifies 
%   this for a three-point BVP.
%
%   YINIT provides a guess for the solution. It must be possible to evaluate 
%   the differential equations and boundary conditions for this guess. 
%   YINIT can be either a vector or a function handle:
%
%   vector:  YINIT(i) is a constant guess for the i-th component Y(i,:) of 
%            the solution at all the mesh points in X.
%
%   function:  YINIT is a function of a scalar x. For example, use 
%              solinit = bvpinit(x,@yfun) if for any x in [a,b], yfun(x) 
%              returns a guess for the solution y(x). For multipoint BVPs, 
%              BVPINIT calls Y = YINIT(X,K) to get an initial guess for the 
%              solution at x in region k. 
%                       
%   SOLINIT = BVPINIT(X,YINIT,PARAMETERS) indicates that the BVP involves 
%   unknown parameters. A guess must be provided for all parameters in the 
%   vector PARAMETERS. 
%
%   See also BVPGET, BVPSET, BVP4C, BVP5C, BVPXTEND, DEVAL, FUNCTION_HANDLE.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2007 The MathWorks, Inc.

% Extend existing solution?  (backwards compatibility)
if isstruct(x)
  solinit = x;  
  if ~isfield(solinit,'solver')
    solinit.solver = 'bvp4c';  % backwards compatibility
  end
  if nargin < 2 || length(v) < 2
    error(message('MATLAB:bvpinit:NoSolInterval'))
  elseif nargin < 3
    solinit = bvpxtend(solinit,v(1),'solution');
    solinit = bvpxtend(solinit,v(2),'solution');      
  else % pass unknown parameters    
    solinit = bvpxtend(solinit,v(1),'solution',parameters);
    solinit = bvpxtend(solinit,v(2),'solution',parameters);
  end   
  % remove solver-specific information
  switch solinit.solver
    case 'bvp4c'
      solinit = rmfield(solinit,{'yp','solver'}); 
    case 'bvp5c'
      solinit = rmfield(solinit,{'idata','solver'});
  end
  return;
end

% Create a guess structure.
N = length(x);
if x(1) == x(N)
  error(message('MATLAB:bvpinit:XSameEndPts'))    
elseif x(1) < x(N)
  if any(diff(x) < 0)
    error(message('MATLAB:bvpinit:IncreasingXNotMonotonic'))
  end
else  % x(1) > x(N)
  if any(diff(x) > 0)
    error(message('MATLAB:bvpinit:DecreasingXNotMonotonic'))
  end
end

if nargin > 2
  params = parameters;
else
  params = [];
end

extraArgs = varargin;

mbcidx = find(diff(x) == 0);  % locate internal interfaces
ismbvp = ~isempty(mbcidx);  
if ismbvp
  Lidx = [1, mbcidx+1]; 
  Ridx = [mbcidx, length(x)];
end

if isnumeric(v) 
  w = v;
else
  if ismbvp
    w = feval(v,x(1),1,extraArgs{:});   % check region 1, only.
  else
    w = feval(v,x(1),extraArgs{:});
  end
end
[m,n] = size(w);
if m == 1
  L = n;
elseif n == 1
  L = m;
else
  error(message('MATLAB:bvpinit:SolGuessNotVector'))
end

yinit = zeros(L,N);
if isnumeric(v)
  yinit = repmat(v(:),1,N);
else 
  if ismbvp
    for region = 1:length(Lidx)
      for i = Lidx(region):Ridx(region)
        w = feval(v,x(i),region,extraArgs{:});
        yinit(:,i) = w(:);
      end  
    end        
  else  
    yinit(:,1) = w(:);
    for i = 2:N
      w = feval(v,x(i),extraArgs{:});
      yinit(:,i) = w(:);
    end
  end  
end

solinit.solver = 'bvpinit';
solinit.x = x(:)';  % row vector
solinit.y = yinit;
if ~isempty(params)  
  solinit.parameters = params;
end
solinit.yinit = v;