| [df,neqn,y,x,h,eps,wt,start,hold,k,kold,crash,phi,p,yp,psi,alpha,beta,sig,v,w,g,phase1,ns,nornd,ksteps,twou,fouru,xold,kprev,ivc,iv,kgi,gi,rpar,ipar]=dsteps(df,neqn,y,x,h,eps,wt,start,hold,k,kold,crash,phi,p,yp,psi,alpha,beta,sig,v,w,g,phase1,ns,nornd,kst |
function [df,neqn,y,x,h,eps,wt,start,hold,k,kold,crash,phi,p,yp,psi,alpha,beta,sig,v,w,g,phase1,ns,nornd,ksteps,twou,fouru,xold,kprev,ivc,iv,kgi,gi,rpar,ipar]=dsteps(df,neqn,y,x,h,eps,wt,start,hold,k,kold,crash,phi,p,yp,psi,alpha,beta,sig,v,w,g,phase1,ns,nornd,ksteps,twou,fouru,xold,kprev,ivc,iv,kgi,gi,rpar,ipar);
persistent absh big erk erkm1 erkm2 erkp1 err firstCall gstr gt hnew i ifail im1 ip1 iq j jv km1 km2 knew kp1 kp2 l limit1 limit2 nsm2 nsp1 nsp2 p5eps r reali realns rho round tau temp1 temp2 temp3 temp4 temp5 temp6 two u ; if isempty(firstCall),firstCall=1;end;
if isempty(jv), jv=0; end;
if isempty(gt), gt=0; end;
%***BEGIN PROLOGUE DSTEPS
%***PURPOSE Integrate a system of first order ordinary differential
% equations one step.
%***LIBRARY SLATEC (DEPAC)
%***CATEGORY I1A1B
%***TYPE doubleprecision (STEPS-S, DSTEPS-D)
%***KEYWORDS ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
% ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
%***AUTHOR Shampine, L. F., (SNLA)
% Gordon, M. K., (SNLA)
% MODIFIED BY H.A. WATTS
%***DESCRIPTION
%
% Written by L. F. Shampine and M. K. Gordon
%
% Abstract
%
% subroutine DSTEPS is normally used indirectly through subroutine
% DDEABM . Because DDEABM suffices for most problems and is much
% easier to use, using it should be considered before using DSTEPS
% alone.
%
% subroutine DSTEPS integrates a system of NEQN first order ordinary
% differential equations one step, normally from X to X+H, using a
% modified divided difference form of the Adams Pece formulas. Local
% extrapolation is used to improve absolute stability and accuracy.
% The code adjusts its order and step size to control the local error
% per unit step in a generalized sense. Special devices are included
% to control roundoff error and to detect when the user is requesting
% too much accuracy.
%
% This code is completely explained and documented in the text,
% Computer Solution of Ordinary Differential Equations, The Initial
% Value Problem by L. F. Shampine and M. K. Gordon.
% Further details on use of this code are available in 'Solving
% Ordinary Differential Equations with ODE, STEP, and INTRP',
% by L. F. Shampine and M. K. Gordon, SLA-73-1060.
%
%
% The parameters represent --
% DF -- subroutine to evaluate derivatives
% NEQN -- number of equations to be integrated
% Y(*) -- solution vector at X
% X -- independent variable
% H -- appropriate step size for next step. Normally determined by
% code
% EPS -- local error tolerance
% WT(*) -- vector of weights for error criterion
% START -- logical variable set true for first step, false
% otherwise
% HOLD -- step size used for last successful step
% K -- appropriate order for next step (determined by code)
% KOLD -- order used for last successful step
% CRASH -- logical variable set true when no step can be taken,
% false otherwise.
% YP(*) -- derivative of solution vector at X after successful
% step
% KSTEPS -- counter on attempted steps
% TWOU -- 2.*U where U is machine unit roundoff quantity
% FOURU -- 4.*U where U is machine unit roundoff quantity
% RPAR,IPAR -- parameter arrays which you may choose to use
% for communication between your program and subroutine F.
% They are not altered or used by DSTEPS.
% The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G,
% W,P,IV and GI are required for the interpolation subroutine SINTRP.
% The remaining variables and arrays are included in the call list
% only to eliminate local retention of variables between calls.
%
% Input to DSTEPS
%
% First call --
%
% The user must provide storage in his calling program for all arrays
% in the call list, namely
%
% DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
% 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
% 2 RPAR(*),IPAR(*)
%
% **Note**
%
% The user must also declare START , CRASH , PHASE1 and NORND
% logical variables and DF an EXTERNAL subroutine, supply the
% subroutine DF(X,Y,YP) to evaluate
% DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN))
% and initialize only the following parameters.
% NEQN -- number of equations to be integrated
% Y(*) -- vector of initial values of dependent variables
% X -- initial value of the independent variable
% H -- nominal step size indicating direction of integration
% and maximum size of step. Must be variable
% EPS -- local error tolerance per step. Must be variable
% WT(*) -- vector of non-zero weights for error criterion
% START -- true
% YP(*) -- vector of initial derivative values
% KSTEPS -- set KSTEPS to zero
% TWOU -- 2.*U where U is machine unit roundoff quantity
% FOURU -- 4.*U where U is machine unit roundoff quantity
% Define U to be the machine unit roundoff quantity by calling
% the function routine D1MACH, U = D1MACH(4), or by
% computing U so that U is the smallest positive number such
% that 1.0+U .GT. 1.0.
%
% DSTEPS requires that the L2 norm of the vector with components
% LOCAL ERROR(L)/WT(L) be less than EPS for a successful step. The
% array WT allows the user to specify an error test appropriate
% for his problem. For example,
% WT(L) = 1.0 specifies absolute error,
% = ABS(Y(L)) error relative to the most recent value of the
% L-th component of the solution,
% = ABS(YP(L)) error relative to the most recent value of
% the L-th component of the derivative,
% = MAX(WT(L),ABS(Y(L))) error relative to the largest
% magnitude of L-th component obtained so far,
% = ABS(Y(L))*RELERR/EPS + ABSERR/EPS specifies a mixed
% relative-absolute test where RELERR is relative
% error, ABSERR is absolute error and EPS =
% MAX(RELERR,ABSERR) .
%
% Subsequent calls --
%
% subroutine DSTEPS is designed so that all information needed to
% continue the integration, including the step size H and the order
% K , is returned with each step. With the exception of the step
% size, the error tolerance, and the weights, none of the parameters
% should be altered. The array WT must be updated after each step
% to maintain relative error tests like those above. Normally the
% integration is continued just beyond the desired endpoint and the
% solution interpolated there with subroutine SINTRP . If it is
% impossible to integrate beyond the endpoint, the step size may be
% reduced to hit the endpoint since the code will not take a step
% larger than the H input. Changing the direction of integration,
% i.e., the sign of H , requires the user set START = true before
% calling DSTEPS again. This is the only situation in which START
% should be altered.
%
% Output from DSTEPS
%
% Successful Step --
%
% The subroutine returns after each successful step with START and
% CRASH set false . X represents the independent variable
% advanced one step of length HOLD from its value on input and Y
% the solution vector at the new value of X . All other parameters
% represent information corresponding to the new X needed to
% continue the integration.
%
% Unsuccessful Step --
%
% When the error tolerance is too small for the machine precision,
% the subroutine returns without taking a step and CRASH = true .
% An appropriate step size and error tolerance for continuing are
% estimated and all other information is restored as upon input
% before returning. To continue with the larger tolerance, the user
% just calls the code again. A restart is neither required nor
% desirable.
%
%***REFERENCES L. F. Shampine and M. K. Gordon, Solving ordinary
% differential equations with ODE, STEP, and INTRP,
% Report SLA-73-1060, Sandia Laboratories, 1973.
%***ROUTINES CALLED D1MACH, DHSTRT
%***REVISION HISTORY (YYMMDD)
% 740101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE DSTEPS
%
if isempty(i), i=0; end;
if isempty(ifail), ifail=0; end;
if isempty(im1), im1=0; end;
if isempty(ip1), ip1=0; end;
if isempty(iq), iq=0; end;
if isempty(j), j=0; end;
if isempty(km1), km1=0; end;
if isempty(km2), km2=0; end;
if isempty(knew), knew=0; end;
if isempty(kp1), kp1=0; end;
if isempty(kp2), kp2=0; end;
if isempty(l), l=0; end;
if isempty(limit1), limit1=0; end;
if isempty(limit2), limit2=0; end;
if isempty(nsm2), nsm2=0; end;
if isempty(nsp1), nsp1=0; end;
if isempty(nsp2), nsp2=0; end;
if isempty(absh), absh=0; end;
if isempty(big), big=0; end;
if isempty(erk), erk=0; end;
if isempty(erkm1), erkm1=0; end;
if isempty(erkm2), erkm2=0; end;
if isempty(erkp1), erkp1=0; end;
if isempty(err), err=0; end;
if isempty(gstr), gstr=zeros(1,13); end;
if isempty(hnew), hnew=0; end;
if isempty(p5eps), p5eps=0; end;
if isempty(r), r=0; end;
if isempty(reali), reali=0; end;
if isempty(realns), realns=0; end;
if isempty(rho), rho=0; end;
if isempty(round), round=0; end;
if isempty(tau), tau=0; end;
if isempty(temp1), temp1=0; end;
if isempty(temp2), temp2=0; end;
if isempty(temp3), temp3=0; end;
if isempty(temp4), temp4=0; end;
if isempty(temp5), temp5=0; end;
if isempty(temp6), temp6=0; end;
if isempty(two), two=zeros(1,13); end;
if isempty(u), u=0; end;
y_shape=size(y);y=reshape(y,1,[]);
wt_shape=size(wt);wt=reshape(wt,1,[]);
phi_orig=phi;phi_shape=[neqn,16];phi=reshape([phi_orig(1:min(prod(phi_shape),numel(phi_orig))),zeros(1,max(0,prod(phi_shape)-numel(phi_orig)))],phi_shape);
p_shape=size(p);p=reshape(p,1,[]);
yp_shape=size(yp);yp=reshape(yp,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
%
if firstCall, two(1) =[2.0d0]; end;
if firstCall, two(2) =[4.0d0]; end;
if firstCall, two(3) =[8.0d0]; end;
if firstCall, two(4) =[16.0d0]; end;
if firstCall, two(5) =[32.0d0]; end;
if firstCall, two(6) =[64.0d0]; end;
if firstCall, two(7)=[128.0d0]; end;
if firstCall, two(8) =[256.0d0]; end;
if firstCall, two(9) =[512.0d0]; end;
if firstCall, two(10) =[1024.0d0]; end;
if firstCall, two(11) =[2048.0d0]; end;
if firstCall, two(12) =[4096.0d0]; end;
if firstCall, two(13)=[8192.0d0]; end;
if firstCall, gstr(1) =[0.5d0]; end;
if firstCall, gstr(2) =[0.0833d0]; end;
if firstCall, gstr(3) =[0.0417d0]; end;
if firstCall, gstr(4) =[0.0264d0]; end;
if firstCall, gstr(5) =[0.0188d0]; end;
if firstCall, gstr(6) =[0.0143d0]; end;
if firstCall, gstr(7) =[0.0114d0]; end;
if firstCall, gstr(8) =[0.00936d0]; end;
if firstCall, gstr(9) =[0.00789d0]; end;
if firstCall, gstr(10) =[0.00679d0]; end;
if firstCall, gstr(11) =[0.00592d0]; end;
if firstCall, gstr(12)=[0.00524d0]; end;
if firstCall, gstr(13)=[0.00468d0]; end;
firstCall=0;
%
% *** BEGIN BLOCK 0 ***
% CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
% PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
% STARTING STEP SIZE.
% ***
%
% IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
%
%***FIRST EXECUTABLE STATEMENT DSTEPS
crash = true;
gt=0;
if( abs(h)>=fouru.*abs(x) )
p5eps = 0.5d0.*eps;
%
% IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
%
round = 0.0d0;
for l = 1 : neqn;
round = round +(y(l)./wt(l)).^2;
end; l = fix(neqn+1);
round = twou.*sqrt(round);
if( p5eps>=round )
crash = false;
g(1) = 1.0d0;
g(2) = 0.5d0;
sig(1) = 1.0d0;
if( start )
%
% INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
%
% CALL DF(X,Y,YP,RPAR,IPAR)
% SUM = 0.0
for l = 1 : neqn;
phi(l,1) = yp(l);
phi(l,2) = 0.0d0;
end; l = fix(neqn+1);
%20 SUM = SUM + (YP(L)/WT(L))**2
% SUM = SQRT(SUM)
% ABSH = ABS(H)
% IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
% H = SIGN(MAX(ABSH,FOURU*ABS(X)),H)
%
[u ]=d1mach(4);
big = sqrt(d1mach(2));
[df,neqn,x,dumvar4,y,yp,wt,dumvar8,u,big,dumvar11,dumvar12,dumvar13,dumvar14,rpar,ipar,h]=dhstrt(df,neqn,x,x+h,y,yp,wt,1,u,big,phi(sub2ind(size(phi),1,3):end),phi(sub2ind(size(phi),1,4):end),phi(sub2ind(size(phi),1,5):end),phi(sub2ind(size(phi),1,6):end),rpar,ipar,h); phi(sub2ind(size(phi),1,3):end)=dumvar11; phi(sub2ind(size(phi),1,4):end)=dumvar12; phi(sub2ind(size(phi),1,5):end)=dumvar13; phi(sub2ind(size(phi),1,6):end)=dumvar14;
%
hold = 0.0d0;
k = 1;
kold = 0;
kprev = 0;
start = false;
phase1 = true;
nornd = true;
if( p5eps<=100.0d0.*round )
nornd = false;
for l = 1 : neqn;
phi(l,15) = 0.0d0;
end; l = fix(neqn+1);
end;
end;
ifail = 0;
% *** end BLOCK 0 ***
%
% *** BEGIN BLOCK 1 ***
% COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING
% THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
% ***
%
gt=0;
while( true );
kp1 = fix(k + 1);
kp2 = fix(k + 2);
km1 = fix(k - 1);
km2 = fix(k - 2);
%
% NS IS THE NUMBER OF DSTEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
% ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE
%
if( h~=hold )
ns = 0;
end;
if( ns<=kold )
ns = fix(ns + 1);
end;
nsp1 = fix(ns + 1);
if( k>=ns )
%
% COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
% ARE CHANGED
%
beta(ns) = 1.0d0;
realns = ns;
alpha(ns) = 1.0d0./realns;
temp1 = h.*realns;
sig(nsp1) = 1.0d0;
if( k>=nsp1 )
for i = nsp1 : k;
im1 = fix(i - 1);
temp2 = psi(im1);
psi(im1) = temp1;
beta(i) = beta(im1).*psi(im1)./temp2;
temp1 = temp2 + h;
alpha(i) = h./temp1;
reali = i;
sig(i+1) = reali.*alpha(i).*sig(i);
end; i = fix(k+1);
end;
psi(k) = temp1;
%
% COMPUTE COEFFICIENTS G(*)
%
% INITIALIZE V(*) AND SET W(*).
%
if( ns>1 )
%
% IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
%
if( k>kprev )
if( ivc==0 )
jv = 1;
temp4 = k.*kp1;
v(k) = 1.0d0./temp4;
w(k) = v(k);
if( k==2 )
kgi = 1;
gi(1) = w(2);
end;
else;
jv = fix(kp1 - iv(ivc));
ivc = fix(ivc - 1);
end;
nsm2 = fix(ns - 2);
if( nsm2>=jv )
for j = jv : nsm2;
i = fix(k - j);
v(i) = v(i) - alpha(j+1).*v(i+1);
w(i) = v(i);
end; j = fix(nsm2+1);
if( i==2 )
kgi = fix(ns - 1);
gi(kgi) = w(2);
end;
end;
end;
%
% UPDATE V(*) AND SET W(*)
%
limit1 = fix(kp1 - ns);
temp5 = alpha(ns);
for iq = 1 : limit1;
v(iq) = v(iq) - temp5.*v(iq+1);
w(iq) = v(iq);
end; iq = fix(limit1+1);
g(nsp1) = w(1);
if( limit1~=1 )
kgi = fix(ns);
gi(kgi) = w(2);
end;
w(limit1+1) = v(limit1+1);
if( k<kold )
ivc = fix(ivc + 1);
iv(ivc) = fix(limit1 + 2);
end;
else;
for iq = 1 : k;
temp3 = iq.*(iq+1);
v(iq) = 1.0d0./temp3;
w(iq) = v(iq);
end; iq = fix(k+1);
ivc = 0;
kgi = 0;
if( k~=1 )
kgi = 1;
gi(1) = w(2);
end;
end;
%
% COMPUTE THE G(*) IN THE WORK VECTOR W(*)
%
nsp2 = fix(ns + 2);
kprev = fix(k);
if( kp1>=nsp2 )
for i = nsp2 : kp1;
limit2 = fix(kp2 - i);
temp6 = alpha(i-1);
for iq = 1 : limit2;
w(iq) = w(iq) - temp6.*w(iq+1);
end; iq = fix(limit2+1);
g(i) = w(1);
end; i = fix(kp1+1);
end;
end;
% *** end BLOCK 1 ***
%
% *** BEGIN BLOCK 2 ***
% PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
% SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
% K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
% ***
%
% INCREMENT COUNTER ON ATTEMPTED DSTEPS
%
ksteps = fix(ksteps + 1);
%
% CHANGE PHI TO PHI STAR
%
if( k>=nsp1 )
for i = nsp1 : k;
temp1 = beta(i);
for l = 1 : neqn;
phi(l,i) = temp1.*phi(l,i);
end; l = fix(neqn+1);
end; i = fix(k+1);
end;
%
% PREDICT SOLUTION AND DIFFERENCES
%
for l = 1 : neqn;
phi(l,kp2) = phi(l,kp1);
phi(l,kp1) = 0.0d0;
p(l) = 0.0d0;
end; l = fix(neqn+1);
for j = 1 : k;
i = fix(kp1 - j);
ip1 = fix(i + 1);
temp2 = g(i);
for l = 1 : neqn;
p(l) = p(l) + temp2.*phi(l,i);
phi(l,i) = phi(l,i) + phi(l,ip1);
end; l = fix(neqn+1);
end; j = fix(k+1);
if( nornd )
for l = 1 : neqn;
p(l) = y(l) + h.*p(l);
end; l = fix(neqn+1);
else;
for l = 1 : neqn;
tau = h.*p(l) - phi(l,15);
p(l) = y(l) + tau;
phi(l,16) =(p(l)-y(l)) - tau;
end; l = fix(neqn+1);
end;
xold = x;
x = x + h;
absh = abs(h);
[x,p,yp,rpar,ipar]=df(x,p,yp,rpar,ipar);
%
% ESTIMATE ERRORS AT ORDERS K,K-1,K-2
%
erkm2 = 0.0d0;
erkm1 = 0.0d0;
erk = 0.0d0;
for l = 1 : neqn;
temp3 = 1.0d0./wt(l);
temp4 = yp(l) - phi(l,1);
if( km2>=0 )
if( km2~=0 )
erkm2 = erkm2 +((phi(l,km1)+temp4).*temp3).^2;
end;
erkm1 = erkm1 +((phi(l,k)+temp4).*temp3).^2;
end;
erk = erk +(temp4.*temp3).^2;
end; l = fix(neqn+1);
if( km2>=0 )
if( km2~=0 )
erkm2 = absh.*sig(km1).*gstr(km2).*sqrt(erkm2);
end;
erkm1 = absh.*sig(k).*gstr(km1).*sqrt(erkm1);
end;
temp5 = absh.*sqrt(erk);
err = temp5.*(g(k)-g(kp1));
erk = temp5.*sig(kp1).*gstr(k);
knew = fix(k);
%
% TEST IF ORDER SHOULD BE LOWERED
%
if( km2>=0 )
if( km2==0 )
if( erkm1<=0.5d0.*erk )
knew = fix(km1);
end;
else;
if( max(erkm1,erkm2)<=erk )
knew = fix(km1);
end;
end;
end;
%
% TEST IF STEP SUCCESSFUL
%
if( err<=eps )
% *** end BLOCK 3 ***
%
% *** BEGIN BLOCK 4 ***
% THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION, EVALUATE
% THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
% DIFFERENCES. DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
% ***
kold = fix(k);
hold = h;
%
% CORRECT AND EVALUATE
%
temp1 = h.*g(kp1);
if( nornd )
for l = 1 : neqn;
temp3 = y(l);
y(l) = p(l) + temp1.*(yp(l)-phi(l,1));
p(l) = temp3;
end; l = fix(neqn+1);
else;
for l = 1 : neqn;
temp3 = y(l);
rho = temp1.*(yp(l)-phi(l,1)) - phi(l,16);
y(l) = p(l) + rho;
phi(l,15) =(y(l)-p(l)) - rho;
p(l) = temp3;
end; l = fix(neqn+1);
end;
[x,y,yp,rpar,ipar]=df(x,y,yp,rpar,ipar);
%
% UPDATE DIFFERENCES FOR NEXT STEP
%
for l = 1 : neqn;
phi(l,kp1) = yp(l) - phi(l,1);
phi(l,kp2) = phi(l,kp1) - phi(l,kp2);
end; l = fix(neqn+1);
for i = 1 : k;
for l = 1 : neqn;
phi(l,i) = phi(l,i) + phi(l,kp1);
end; l = fix(neqn+1);
end; i = fix(k+1);
%
% ESTIMATE ERROR AT ORDER K+1 UNLESS:
% IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
% ALREADY DECIDED TO LOWER ORDER,
% STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
%
erkp1 = 0.0d0;
if( knew==km1 || k==12 )
phase1 = false;
end;
while (1);
if( ~(phase1) )
if( knew~=km1 )
if( kp1>ns )
gt=1;
break;
end;
for l = 1 : neqn;
erkp1 = erkp1 +(phi(l,kp2)./wt(l)).^2;
end; l = fix(neqn+1);
erkp1 = absh.*gstr(kp1).*sqrt(erkp1);
%
% USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
% FOR NEXT STEP
%
if( k>1 )
if( erkm1>min(erk,erkp1) )
if( erkp1<erk && k~=12 )
break;
end;
gt=1;
break;
end;
elseif( erkp1>=0.5d0.*erk ) ;
gt=1;
break;
else;
break;
end;
end;
%
% LOWER ORDER
%
k = fix(km1);
erk = erkm1;
gt=1;
break;
end;
break;
end;
if(gt==1)
break;
end;
%
% HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
% BEEN LOWERED IN BLOCK 2. THUS ORDER IS TO BE RAISED
%
% RAISE ORDER
%
k = fix(kp1);
erk = erkp1;
gt=1;
break;
else;
% *** end BLOCK 2 ***
%
% *** BEGIN BLOCK 3 ***
% THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) .
% IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE
% THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR
% TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
% PRECISION.
% ***
%
% RESTORE X, PHI(*,*) AND PSI(*)
%
phase1 = false;
x = xold;
for i = 1 : k;
temp1 = 1.0d0./beta(i);
ip1 = fix(i + 1);
for l = 1 : neqn;
phi(l,i) = temp1.*(phi(l,i)-phi(l,ip1));
end; l = fix(neqn+1);
end; i = fix(k+1);
if( k>=2 )
for i = 2 : k;
psi(i-1) = psi(i) - h;
end; i = fix(k+1);
end;
%
% ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, USE OPTIMAL STEP
% SIZE
%
ifail = fix(ifail + 1);
temp2 = 0.5d0;
if( ifail>=3 )
if( ifail~=3 )
if( p5eps<0.25d0.*erk )
temp2 = sqrt(p5eps./erk);
end;
end;
knew = 1;
end;
h = temp2.*h;
k = fix(knew);
ns = 0;
if( abs(h)<fouru.*abs(x) )
break;
end;
end;
end;
if(gt==0)
crash = true;
h = (abs(fouru.*abs(x)).*sign(h));
eps = eps + eps;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
return;
end;
%
% WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
%
hnew = h + h;
if( ~(phase1) )
if( p5eps<erk.*two(k+1) )
hnew = h;
if( p5eps<erk )
temp2 = k + 1;
r =(p5eps./erk).^(1.0d0./temp2);
hnew = absh.*max(0.5d0,min(0.9d0,r));
hnew = (abs(max(hnew,fouru.*abs(x))).*sign(h));
end;
end;
end;
h = hnew;
else;
eps = 2.0d0.*round.*(1.0d0+fouru);
end;
else;
h = (abs(fouru.*abs(x)).*sign(h));
end;
% *** end BLOCK 4 ***
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
end %subroutine dsteps
%DECK DSTOD
|
|
