bvpinit_shi.m

用matlab编写的梁几何大变形程序包括静力和动力程序

function solinit = bvpinit(x,v,parameters,varargin)

if isstruct(x)
  if nargin < 2 || length(v) < 2
    error('MATLAB:bvpinit:NoSolInterval',...
          'Did not specify [ANEW BNEW] in BVPINIT(SOL,[ANEW BNEW]).')
  elseif nargin < 3
    parameters = [];
  end   
  solinit = bvpxtrp(x,v,parameters);
  return;
end

% Create a guess structure.
N = length(x);
if x(1) == x(N)
  error('MATLAB:bvpinit:XSameEndPts',...
        'The entries of x must satisfy a = x(1) ~= x(end) = b.')    
elseif x(1) < x(N)
  if any(diff(x) < 0)
    error('MATLAB:bvpinit:IncreasingXNotMonotonic',...
          'The entries of x must satisfy a = x(1) < x(2) < ... < x(end) = b.')
  end
else  % x(1) > x(N)
  if any(diff(x) > 0)
    error('MATLAB:bvpinit:DecreasingXNotMonotonic',...
          'The entries of x must satisfy a = x(1) > x(2) > ... > x(end) = b.')
  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('MATLAB:bvpinit:SolGuessNotVector',...
        'The guess for the solution must return a vector.')
end

yinit = zeros(L,N);
if isnumeric(v)
  yinit = repmat(v(:),1,N);
else 
  if ismbvp
    for region = 1:nregions
      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.x = x(:)';  % row vector
solinit.y = yinit;
if ~isempty(params)  
  solinit.parameters = params;
end

%---------------------------------------------------------------------------

function solxtrp = bvpxtrp(sol,v,parameters)
% Extend a solution SOL on [sol.x(1),sol.x(end)]
% to a guess for [v(1),v(end)] by extrapolation.

a = sol.x(1);
b = sol.x(end);
  
anew = v(1);
bnew = v(end);

if (a < b && (anew > a || bnew < b)) || ...
   (a > b && (anew < a || bnew > b))
    
  msg = sprintf(['The call BVPINIT(SOL,[ANEW BNEW]) must have\n',...
                'ANEW <= SOL.x(1) < SOL.x(end) <= BNEW  or \n',...
                'ANEW >= SOL.x(1) > SOL.x(end) >= BNEW. \n']);
  error('MATLAB:bvpinit:bvpxtrp:SolInterval', msg);            
   
end   
   
solxtrp.x = sol.x;
solxtrp.y = sol.y;
if abs(anew - a) <= 100*eps*max(abs(anew),abs(a))
    solxtrp.x(1) = anew;
else
    S = Hermite(sol.x(1),sol.y(:,1),sol.yp(:,1),...
                sol.x(2),sol.y(:,2),sol.yp(:,2),anew);
    solxtrp.x = [anew solxtrp.x];
    solxtrp.y = [S solxtrp.y];
end
if abs(bnew - b) <= 100*eps*max(abs(bnew),abs(b))
    solxtrp.x(end) = bnew;
else  
    S = Hermite(sol.x(end-1),sol.y(:,end-1),sol.yp(:,end-1),...
                sol.x(end),sol.y(:,end),sol.yp(:,end),bnew);
    solxtrp.x = [solxtrp.x bnew];
    solxtrp.y = [solxtrp.y S];
end

if ~isempty(parameters)
  solxtrp.parameters = parameters;
elseif isfield(sol,'parameters')
  solxtrp.parameters = sol.parameters;
end

%---------------------------------------------------------------------------

function S = Hermite(x1,y1,yp1,x2,y2,yp2,xi)
% Evaluate cubic Hermite interpolant at xi.
h = x2 - x1;
s = (xi - x1)/h;
A1 = (s - 1)^2 * (1 + 2*s);
A2 = (3 - 2*s)*s^2;
B1 = h*s*(s - 1)^2;
B2 = h*s^2 *(s - 1);
S = A1*y1 + A2*y2 + B1*yp1 + B2*yp2;