| [df,neq,t,y,tout,info,rtol,atol,idid,h,tolfac,yp,f1,f2,f3,f4,f5,ys,told,dtsign,u26,rer,init,ksteps,kop,iquit,stiff,nonstf,ntstep,nstifs,rpar,ipar]=drkfs(df,neq,t,y,tout,info,rtol,atol,idid,h,tolfac,yp,f1,f2,f3,f4,f5,ys,told,dtsign,u26,rer,init,ksteps,kop, |
function [df,neq,t,y,tout,info,rtol,atol,idid,h,tolfac,yp,f1,f2,f3,f4,f5,ys,told,dtsign,u26,rer,init,ksteps,kop,iquit,stiff,nonstf,ntstep,nstifs,rpar,ipar]=drkfs(df,neq,t,y,tout,info,rtol,atol,idid,h,tolfac,yp,f1,f2,f3,f4,f5,ys,told,dtsign,u26,rer,init,ksteps,kop,iquit,stiff,nonstf,ntstep,nstifs,rpar,ipar);
%***BEGIN PROLOGUE DRKFS
%***SUBSIDIARY
%***PURPOSE Subsidiary to DDERKF
%***LIBRARY SLATEC
%***TYPE doubleprecision (DERKFS-S, DRKFS-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% Fehlberg Fourth-Fifth Order Runge-Kutta Method
% **********************************************************************
%
% DRKFS integrates a system of first order ordinary differential
% equations as described in the comments for DDERKF .
%
% The arrays YP,F1,F2,F3,F4,F5,and YS (of length at least NEQ)
% appear in the call list for variable dimensioning purposes.
%
% The variables H,TOLFAC,TOLD,DTSIGN,U26,RER,INIT,KSTEPS,KOP,IQUIT,
% STIFF,NONSTF,NTSTEP, and NSTIFS are used internally by the code
% and appear in the call list to eliminate local retention of
% variables between calls. Accordingly, these variables and the
% array YP should not be altered.
% Items of possible interest are
% H - An appropriate step size to be used for the next step
% TOLFAC - Factor of change in the tolerances
% YP - Derivative of solution vector at T
% KSTEPS - Counter on the number of steps attempted
%
% **********************************************************************
%
%***SEE ALSO DDERKF
%***ROUTINES CALLED D1MACH, DFEHL, DHSTRT, DHVNRM, XERMSG
%***REVISION HISTORY (YYMMDD)
% 820301 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 891024 Changed references from DVNORM to DHVNRM. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 900510 Convert XERRWV calls to XERMSG calls, change GO TOs to
% IF-THEN-ELSEs. (RWC)
% 910722 Updated AUTHOR section. (ALS)
%***end PROLOGUE DRKFS
%
persistent a big dt dy ee eeoet es estiff esttol et firstCall gt hfaild hmin k ktol mxkop mxstep natolp nrtolp output remin s tol u ute xern1 xern3 xern4 yavg ; if isempty(firstCall),firstCall=1;end;
if isempty(k), k=0; end;
if isempty(ktol), ktol=0; end;
if isempty(mxkop), mxkop=0; end;
if isempty(mxstep), mxstep=0; end;
if isempty(natolp), natolp=0; end;
if isempty(nrtolp), nrtolp=0; end;
if isempty(gt), gt=zeros(1,2); end;
if isempty(a), a=0; end;
if isempty(big), big=0; end;
if isempty(dt), dt=0; end;
if isempty(dy), dy=0; end;
if isempty(ee), ee=0; end;
if isempty(eeoet), eeoet=0; end;
if isempty(es), es=0; end;
if isempty(estiff), estiff=0; end;
if isempty(esttol), esttol=0; end;
if isempty(et), et=0; end;
if isempty(hmin), hmin=0; end;
if isempty(remin), remin=0; end;
if isempty(s), s=0; end;
if isempty(tol), tol=0; end;
if isempty(u), u=0; end;
if isempty(ute), ute=0; end;
if isempty(yavg), yavg=0; end;
if isempty(hfaild), hfaild=false; end;
if isempty(output), output=false; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,16); end;
if isempty(xern4), xern4=repmat(' ',1,16); end;
%
y_shape=size(y);y=reshape(y,1,[]);
yp_shape=size(yp);yp=reshape(yp,1,[]);
f1_shape=size(f1);f1=reshape(f1,1,[]);
f2_shape=size(f2);f2=reshape(f2,1,[]);
f3_shape=size(f3);f3=reshape(f3,1,[]);
f4_shape=size(f4);f4=reshape(f4,1,[]);
f5_shape=size(f5);f5=reshape(f5,1,[]);
ys_shape=size(ys);ys=reshape(ys,1,[]);
rtol_shape=size(rtol);rtol=reshape(rtol,1,[]);
atol_shape=size(atol);atol=reshape(atol,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
%
%
% ..................................................................
%
% A FIFTH ORDER METHOD WILL GENERALLY NOT BE CAPABLE OF DELIVERING
% ACCURACIES NEAR LIMITING PRECISION ON COMPUTERS WITH LONG
% WORDLENGTHS. TO PROTECT AGAINST LIMITING PRECISION DIFFICULTIES
% ARISING FROM UNREASONABLE ACCURACY REQUESTS, AN APPROPRIATE
% TOLERANCE THRESHOLD REMIN IS ASSIGNED FOR THIS METHOD. THIS
% VALUE SHOULD NOT BE CHANGED ACROSS DIFFERENT MACHINES.
%
if firstCall, remin=[1.0d-12]; end;
%
% ..................................................................
%
% THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE
% NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MXSTEP, THE
% COUNTER IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE
% EXCESSIVE WORK.
%
if firstCall, mxstep=[500]; end;
%
% ..................................................................
%
% INEFFICIENCY CAUSED BY TOO FREQUENT OUTPUT IS MONITORED BY
% COUNTING THE NUMBER OF STEP SIZES WHICH ARE SEVERELY SHORTENED
% DUE SOLELY TO THE CHOICE OF OUTPUT POINTS. WHEN THE NUMBER OF
% ABUSES EXCEED MXKOP, THE COUNTER IS RESET TO ZERO AND THE USER
% IS INFORMED ABOUT POSSIBLE MISUSE OF THE CODE.
%
if firstCall, mxkop=[100]; end;
firstCall=0;
%
% ..................................................................
%
%***FIRST EXECUTABLE STATEMENT DRKFS
gt(:)=0;
if( info(1)==0 )
%
% ON THE FIRST CALL , PERFORM INITIALIZATION --
% DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING THE
% function ROUTINE D1MACH. THE USER MUST MAKE SURE THAT THE
% VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
%
[u ]=d1mach(4);
% -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS
u26 = 26.0d0.*u;
rer = 2.0d0.*u + remin;
% -- SET TERMINATION FLAG
iquit = 0;
% -- SET INITIALIZATION INDICATOR
init = 0;
% -- SET COUNTER FOR IMPACT OF OUTPUT POINTS
kop = 0;
% -- SET COUNTER FOR ATTEMPTED STEPS
ksteps = 0;
% -- SET INDICATORS FOR STIFFNESS DETECTION
stiff = false;
nonstf = false;
% -- SET STEP COUNTERS FOR STIFFNESS DETECTION
ntstep = 0;
nstifs = 0;
% -- RESET INFO(1) FOR SUBSEQUENT CALLS
info(1) = 1;
end;
%
%.......................................................................
%
% CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY
%
if( info(1)~=0 && info(1)~=1 )
xern1=sprintf(['%8i'], info(1));
xermsg('SLATEC','DRKFS',['IN DDERKF, INFO(1) MUST BE SET TO 0 ',['FOR THE START OF A NEW PROBLEM, AND MUST BE SET TO 1 ',['FOLLOWING AN INTERRUPTED TASK. YOU ARE ATTEMPTING TO ',['CONTINUE THE INTEGRATION ILLEGALLY BY CALLING THE CODE ',['WITH INFO(1) = ',xern1]]]]],3,1);
idid = -33;
end;
%
if( info(2)~=0 && info(2)~=1 )
xern1=sprintf(['%8i'], info(2));
xermsg('SLATEC','DRKFS',['IN DDERKF, INFO(2) MUST BE 0 OR 1 ',['INDICATING SCALAR AND VECTOR ERROR TOLERANCES, ',['RESPECTIVELY. YOU HAVE CALLED THE CODE WITH INFO(2) = ',xern1]]],4,1);
idid = -33;
end;
%
if( info(3)~=0 && info(3)~=1 )
xern1=sprintf(['%8i'], info(3));
xermsg('SLATEC','DRKFS',['IN DDERKF, INFO(3) MUST BE 0 OR 1 ',['INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT MODE OF ',['INTEGRATION, RESPECTIVELY. YOU HAVE CALLED THE CODE ',['WITH INFO(3) = ',xern1]]]],5,1);
idid = -33;
end;
%
if( neq<1 )
xern1=sprintf(['%8i'], neq);
xermsg('SLATEC','DRKFS',['IN DDERKF, THE NUMBER OF EQUATIONS ',['NEQ MUST BE A POSITIVE INTEGER. YOU HAVE CALLED THE ',['CODE WITH NEQ = ',xern1]]],6,1);
idid = -33;
end;
%
nrtolp = 0;
natolp = 0;
for k = 1 : neq;
if( nrtolp==0 && rtol(k)<0.0d0 )
xern1=sprintf(['%8i'], k);
xern3=sprintf([repmat('%15.6f',1,1)], rtol(k));
xermsg('SLATEC','DRKFS',['IN DDERKF, THE RELATIVE ERROR ',['TOLERANCES RTOL MUST BE NON-NEGATIVE. YOU HAVE ',['CALLED THE CODE WITH RTOL(',[xern1,[') = ',[xern3,['. IN THE CASE OF VECTOR ERROR TOLERANCES, ','NO FURTHER CHECKING OF RTOL COMPONENTS IS DONE.']]]]]]],7,1);
idid = -33;
nrtolp = 1;
end;
%
if( natolp==0 && atol(k)<0.0d0 )
xern1=sprintf(['%8i'], k);
xern3=sprintf([repmat('%15.6f',1,1)], atol(k));
xermsg('SLATEC','DRKFS',['IN DDERKF, THE ABSOLUTE ERROR ',['TOLERANCES ATOL MUST BE NON-NEGATIVE. YOU HAVE ',['CALLED THE CODE WITH ATOL(',[xern1,[') = ',[xern3,['. IN THE CASE OF VECTOR ERROR TOLERANCES, ','NO FURTHER CHECKING OF ATOL COMPONENTS IS DONE.']]]]]]],8,1);
idid = -33;
natolp = 1;
end;
%
if( info(2)==0 )
break;
end;
if( natolp>0 && nrtolp>0 )
break;
end;
end;
%
%
% CHECK SOME CONTINUATION POSSIBILITIES
%
if( init~=0 )
if( t==tout )
xern3=sprintf([repmat('%15.6f',1,1)], t);
xermsg('SLATEC','DRKFS',['IN DDERKF, YOU HAVE CALLED THE ',['CODE WITH T = TOUT = ',[xern3,['$$THIS IS NOT ','ALLOWED ON CONTINUATION CALLS.']]]],9,1);
idid = -33;
end;
%
if( t~=told )
xern3=sprintf([repmat('%15.6f',1,1)], told);
xern4=sprintf([repmat('%15.6f',1,1)], t);
xermsg('SLATEC','DRKFS',['IN DDERKF, YOU HAVE CHANGED THE ',['VALUE OF T FROM ',[xern3,[' TO ',[xern4,'$$THIS IS NOT ALLOWED ON CONTINUATION CALLS.']]]]],10,1);
idid = -33;
end;
%
if( init~=1 )
if( dtsign.*(tout-t)<0.0d0 )
xern3=sprintf([repmat('%15.6f',1,1)], tout);
xermsg('SLATEC','DRKFS',['IN DDERKF, BY CALLING THE CODE WITH TOUT = ',[xern3,[' YOU ARE ATTEMPTING TO CHANGE THE ',['DIRECTION OF INTEGRATION.$$THIS IS NOT ALLOWED ','WITHOUT RESTARTING.']]]],11,1);
idid = -33;
end;
end;
end;
%
% INVALID INPUT DETECTED
%
while (1);
if( idid==(-33) )
if( iquit~=(-33) )
iquit = -33;
gt(2)=1;
break;
else;
xermsg('SLATEC','DRKFS',['IN DDERKF, INVALID INPUT WAS ',['DETECTED ON SUCCESSIVE ENTRIES. IT IS IMPOSSIBLE ',['TO PROCEED BECAUSE YOU HAVE NOT CORRECTED THE ','PROBLEM, SO EXECUTION IS BEING TERMINATED.']]],12,2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
f1_shape=zeros(f1_shape);f1_shape(:)=f1(1:numel(f1_shape));f1=f1_shape;
f2_shape=zeros(f2_shape);f2_shape(:)=f2(1:numel(f2_shape));f2=f2_shape;
f3_shape=zeros(f3_shape);f3_shape(:)=f3(1:numel(f3_shape));f3=f3_shape;
f4_shape=zeros(f4_shape);f4_shape(:)=f4(1:numel(f4_shape));f4=f4_shape;
f5_shape=zeros(f5_shape);f5_shape(:)=f5(1:numel(f5_shape));f5=f5_shape;
ys_shape=zeros(ys_shape);ys_shape(:)=ys(1:numel(ys_shape));ys=ys_shape;
rtol_shape=zeros(rtol_shape);rtol_shape(:)=rtol(1:numel(rtol_shape));rtol=rtol_shape;
atol_shape=zeros(atol_shape);atol_shape(:)=atol(1:numel(atol_shape));atol=atol_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;
end;
%
% ............................................................
%
% RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND
% INTERPRETED AS ASKING FOR THE MOST ACCURATE SOLUTION
% POSSIBLE. IN THIS CASE, THE RELATIVE ERROR TOLERANCE
% RTOL IS RESET TO THE SMALLEST VALUE RER WHICH IS LIKELY
% TO BE REASONABLE FOR THIS METHOD AND MACHINE.
%
for k = 1 : neq;
if( rtol(k)+atol(k)<=0.0d0 )
rtol(k) = rer;
idid = -2;
end;
% ...EXIT
if( info(2)==0 )
break;
end;
end;
%
if( idid~=(-2) )
%
% BRANCH ON STATUS OF INITIALIZATION INDICATOR
% INIT=0 MEANS INITIAL DERIVATIVES AND
% STARTING STEP SIZE
% NOT YET COMPUTED
% INIT=1 MEANS STARTING STEP SIZE NOT YET
% COMPUTED INIT=2 MEANS NO FURTHER
% INITIALIZATION REQUIRED
%
if( init==0 )
%
% ................................................
%
% MORE INITIALIZATION --
% -- EVALUATE INITIAL
% DERIVATIVES
%
init = 1;
a = t;
[a,y,yp,rpar,ipar]=df(a,y,yp,rpar,ipar);
if( t==tout )
%
% INTERVAL MODE
idid = 2;
t = tout;
told = t;
% .....................EXIT
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
f1_shape=zeros(f1_shape);f1_shape(:)=f1(1:numel(f1_shape));f1=f1_shape;
f2_shape=zeros(f2_shape);f2_shape(:)=f2(1:numel(f2_shape));f2=f2_shape;
f3_shape=zeros(f3_shape);f3_shape(:)=f3(1:numel(f3_shape));f3=f3_shape;
f4_shape=zeros(f4_shape);f4_shape(:)=f4(1:numel(f4_shape));f4=f4_shape;
f5_shape=zeros(f5_shape);f5_shape(:)=f5(1:numel(f5_shape));f5=f5_shape;
ys_shape=zeros(ys_shape);ys_shape(:)=ys(1:numel(ys_shape));ys=ys_shape;
rtol_shape=zeros(rtol_shape);rtol_shape(:)=rtol(1:numel(rtol_shape));rtol=rtol_shape;
atol_shape=zeros(atol_shape);atol_shape(:)=atol(1:numel(atol_shape));atol=atol_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;
% ......EXIT
elseif( init~=1 ) ;
break;
% .........EXIT
end;
%
% -- SET SIGN OF INTEGRATION DIRECTION AND
% -- ESTIMATE STARTING STEP SIZE
%
init = 2;
dtsign = (abs(1.0d0).*sign(tout-t));
[u ]=d1mach(4);
big = sqrt(d1mach(2));
ute = u.^0.375d0;
dy = ute.*dhvnrm(y,neq);
if( dy==0.0d0 )
dy = ute;
end;
ktol = 1;
for k = 1 : neq;
if( info(2)==1 )
ktol = fix(k);
end;
tol = rtol(ktol).*abs(y(k)) + atol(ktol);
if( tol==0.0d0 )
tol = dy.*rtol(ktol);
end;
f1(k) = tol;
end; k = fix(neq+1);
%
[df,neq,t,tout,y,yp,f1,dumvar8,u,big,f2,f3,f4,f5,rpar,ipar,h]=dhstrt(df,neq,t,tout,y,yp,f1,4,u,big,f2,f3,f4,f5,rpar,ipar,h);
else;
%
% RTOL=ATOL=0 ON INPUT, SO RTOL WAS CHANGED TO A
% SMALL POSITIVE VALUE
tolfac = 1.0d0;
gt(2)=1;
break;
end;
break;
end;
%
% ......................................................
%
% SET STEP SIZE FOR INTEGRATION IN THE DIRECTION
% FROM T TO TOUT AND SET OUTPUT POINT INDICATOR
%
while(gt(2)==0);
gt(:)=0;
dt = tout - t;
h = (abs(h).*sign(dt));
output = false;
%
% TEST TO SEE IF DDERKF IS BEING SEVERELY IMPACTED BY
% TOO MANY OUTPUT POINTS
%
if( abs(h)>=2.0d0.*abs(dt) )
kop = fix(kop + 1);
end;
if( kop<=mxkop )
%
if( abs(dt)>u26.*abs(t) )
% BEGIN BLOCK PERMITTING ...EXITS TO 490
%
% *********************************************
% *********************************************
% STEP BY STEP INTEGRATION
%
% BEGIN BLOCK PERMITTING ...EXITS TO 480
while (1);
gt(:)=0;
hfaild = false;
%
% TO PROTECT AGAINST IMPOSSIBLE ACCURACY
% REQUESTS, COMPUTE A TOLERANCE FACTOR
% BASED ON THE REQUESTED ERROR TOLERANCE
% AND A LEVEL OF ACCURACY ACHIEVABLE AT
% LIMITING PRECISION
%
tolfac = 0.0d0;
ktol = 1;
for k = 1 : neq;
if( info(2)==1 )
ktol = fix(k);
end;
et = rtol(ktol).*abs(y(k)) + atol(ktol);
if( et>0.0d0 )
tolfac = max(tolfac,abs(y(k)).*(rer./et));
else;
tolfac = max(tolfac,rer./rtol(ktol));
end;
end; k = fix(neq+1);
if( tolfac<=1.0d0 )
%
% SET SMALLEST ALLOWABLE STEP SIZE
%
hmin = u26.*abs(t);
%
% ADJUST STEP SIZE IF NECESSARY TO HIT
% THE OUTPUT POINT -- LOOK AHEAD TWO
% STEPS TO AVOID DRASTIC CHANGES IN THE
% STEP SIZE AND THUS LESSEN THE IMPACT OF
% OUTPUT POINTS ON THE CODE. STRETCH THE
% STEP SIZE BY, AT MOST, AN AMOUNT EQUAL
% TO THE SAFETY FACTOR OF 9/10.
%
dt = tout - t;
if( abs(dt)<2.0d0.*abs(h) )
if( abs(dt)>abs(h)./0.9d0 )
%
h = 0.5d0.*dt;
else;
%
% THE NEXT STEP, IF SUCCESSFUL,
% WILL COMPLETE THE INTEGRATION TO
% THE OUTPUT POINT
%
output = true;
h = dt;
end;
end;
%
%
% ***************************************
% CORE INTEGRATOR FOR TAKING A
% SINGLE STEP
% ***************************************
% TO AVOID PROBLEMS WITH ZERO
% CROSSINGS, RELATIVE ERROR IS
% MEASURED USING THE AVERAGE OF THE
% MAGNITUDES OF THE SOLUTION AT THE
% BEGINNING AND END OF A STEP.
% THE ERROR ESTIMATE FORMULA HAS
% BEEN GROUPED TO CONTROL LOSS OF
% SIGNIFICANCE.
% LOCAL ERROR ESTIMATES FOR A FIRST
% ORDER METHOD USING THE SAME
% STEP SIZE AS THE FEHLBERG METHOD
% ARE CALCULATED AS PART OF THE
% TEST FOR STIFFNESS.
% TO DISTINGUISH THE VARIOUS
% ARGUMENTS, H IS NOT PERMITTED
% TO BECOME SMALLER THAN 26 UNITS OF
% ROUNDOFF IN T. PRACTICAL LIMITS
% ON THE CHANGE IN THE STEP SIZE ARE
% ENFORCED TO SMOOTH THE STEP SIZE
% SELECTION PROCESS AND TO AVOID
% EXCESSIVE CHATTERING ON PROBLEMS
% HAVING DISCONTINUITIES. TO
% PREVENT UNNECESSARY FAILURES, THE
% CODE USES 9/10 THE STEP SIZE
% IT ESTIMATES WILL SUCCEED.
% AFTER A STEP FAILURE, THE STEP
% SIZE IS NOT ALLOWED TO INCREASE
% FOR THE NEXT ATTEMPTED STEP. THIS
% MAKES THE CODE MORE EFFICIENT ON
% PROBLEMS HAVING DISCONTINUITIES
% AND MORE EFFECTIVE IN GENERAL
% SINCE LOCAL EXTRAPOLATION IS BEING
% USED AND EXTRA CAUTION SEEMS
% WARRANTED.
% .......................................
%
% MONITOR NUMBER OF STEPS ATTEMPTED
%
while (1);
gt(:)=0;
if( ksteps<=mxstep )
%
% ADVANCE AN APPROXIMATE SOLUTION OVER
% ONE STEP OF LENGTH H
%
[df,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar]=dfehl(df,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar);
ksteps = fix(ksteps + 1);
%
% ....................................
%
% COMPUTE AND TEST ALLOWABLE
% TOLERANCES VERSUS LOCAL ERROR
% ESTIMATES. NOTE THAT RELATIVE
% ERROR IS MEASURED WITH RESPECT
% TO THE AVERAGE OF THE
% MAGNITUDES OF THE SOLUTION AT
% THE BEGINNING AND END OF THE
% STEP. LOCAL ERROR ESTIMATES
% FOR A SPECIAL FIRST ORDER
% METHOD ARE CALCULATED ONLY WHEN
% THE STIFFNESS DETECTION IS
% TURNED ON.
%
eeoet = 0.0d0;
estiff = 0.0d0;
ktol = 1;
for k = 1 : neq;
yavg = 0.5d0.*(abs(y(k))+abs(ys(k)));
if( info(2)==1 )
ktol = fix(k);
end;
et = rtol(ktol).*yavg + atol(ktol);
if( et>0.0d0 )
%
ee = abs((-2090.0d0.*yp(k)+(21970.0d0.*f3(k)-15048.0d0.*f4(k)))+(22528.0d0.*f2(k)-27360.0d0.*f5(k)));
if( ~(stiff || nonstf) )
es = abs(h.*(0.055455d0.*yp(k)-0.035493d0.*f1(k)-0.036571d0.*f2(k)+0.023107d0.*f3(k)-0.009515d0.*f4(k)+0.003017d0.*f5(k)));
estiff = max(estiff,es./et);
end;
eeoet = max(eeoet,ee./et);
else;
%
% PURE RELATIVE ERROR INAPPROPRIATE WHEN SOLUTION
% VANISHES
idid = -3;
% ...........................EXIT
gt(2)=1;
break;
end;
end;
if(gt(2)==1)
break;
end;
%
esttol = abs(h).*eeoet./752400.0d0;
%
% ...EXIT
if( esttol<=1.0d0 )
%
% .......................................
%
% SUCCESSFUL STEP
% STORE SOLUTION AT T+H
% AND EVALUATE
% DERIVATIVES THERE
%
t = t + h;
for k = 1 : neq;
y(k) = ys(k);
end; k = fix(neq+1);
a = t;
[a,y,yp,rpar,ipar]=df(a,y,yp,rpar,ipar);
%
% CHOOSE NEXT STEP SIZE
% THE INCREASE IS LIMITED TO A FACTOR OF
% 5 IF STEP FAILURE HAS JUST OCCURRED,
% NEXT
% STEP SIZE IS NOT ALLOWED TO INCREASE
%
s = 5.0d0;
if( esttol>1.889568d-4 )
s = 0.9d0./esttol.^0.2d0;
end;
if( hfaild )
s = min(s,1.0d0);
end;
h = (abs(max(s.*abs(h),hmin)).*sign(h));
%
% .......................................
%
% CHECK FOR STIFFNESS (IF NOT
% ALREADY DETECTED)
%
% IN A SEQUENCE OF 50 SUCCESSFUL
% STEPS BY THE FEHLBERG METHOD, 25
% SUCCESSFUL STEPS BY THE FIRST
% ORDER METHOD INDICATES STIFFNESS
% AND TURNS THE TEST OFF. IF 26
% FAILURES BY THE FIRST ORDER METHOD
% OCCUR, THE TEST IS TURNED OFF
% UNTIL THIS SEQUENCE OF 50 STEPS BY
% THE FEHLBERG METHOD IS COMPLETED.
%
% ...EXIT
if( ~(stiff) )
ntstep = fix(rem(ntstep+1,50));
if( ntstep==1 )
nonstf = false;
end;
% ...EXIT
if( ~(nonstf) )
if( estiff<=1.0d0 )
%
% SUCCESSFUL STEP WITH FIRST ORDER
% METHOD
nstifs = fix(nstifs + 1);
% TURN TEST OFF AFTER 25 INDICATIONS
% OF STIFFNESS
if( nstifs==25 )
stiff = true;
end;
%
% UNSUCCESSFUL STEP WITH FIRST ORDER
% METHOD
elseif( ntstep-nstifs>25 ) ;
% TURN STIFFNESS DETECTION OFF FOR THIS BLOCK OF
% FIFTY STEPS
nonstf = true;
% RESET STIFF STEP COUNTER
nstifs = 0;
end;
end;
end;
%
% ******************************************
% end OF CORE INTEGRATOR
% ******************************************
%
%
% SHOULD WE TAKE ANOTHER STEP
%
% ......EXIT
if( ~(output) )
if( info(3)==0 )
gt(1)=1;
%to120
break;
end;
%
% *********************************************
% *********************************************
%
% INTEGRATION SUCCESSFULLY COMPLETED
%
% ONE-STEP MODE
idid = 1;
told = t;
% .....................EXIT
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
f1_shape=zeros(f1_shape);f1_shape(:)=f1(1:numel(f1_shape));f1=f1_shape;
f2_shape=zeros(f2_shape);f2_shape(:)=f2(1:numel(f2_shape));f2=f2_shape;
f3_shape=zeros(f3_shape);f3_shape(:)=f3(1:numel(f3_shape));f3=f3_shape;
f4_shape=zeros(f4_shape);f4_shape(:)=f4(1:numel(f4_shape));f4=f4_shape;
f5_shape=zeros(f5_shape);f5_shape(:)=f5(1:numel(f5_shape));f5=f5_shape;
ys_shape=zeros(ys_shape);ys_shape(:)=ys(1:numel(ys_shape));ys=ys_shape;
rtol_shape=zeros(rtol_shape);rtol_shape(:)=rtol(1:numel(rtol_shape));rtol=rtol_shape;
atol_shape=zeros(atol_shape);atol_shape(:)=atol(1:numel(atol_shape));atol=atol_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;
%
% ....................................
%
% UNSUCCESSFUL STEP
%
elseif( abs(h)>hmin ) ;
%
% REDUCE THE STEP SIZE , TRY AGAIN
% THE DECREASE IS LIMITED TO A FACTOR
% OF 1/10
%
hfaild = true;
output = false;
s = 0.1d0;
if( esttol<59049.0d0 )
s = 0.9d0./esttol.^0.2d0;
end;
h = (abs(max(s.*abs(h),hmin)).*sign(h));
continue;
else;
%
% REQUESTED ERROR UNATTAINABLE AT SMALLEST
% ALLOWABLE STEP SIZE
tolfac = 1.69d0.*esttol;
idid = -2;
% ........................EXIT
gt(2)=1;
break;
end;
else;
%
% A SIGNIFICANT AMOUNT OF WORK HAS
% BEEN EXPENDED
idid = -1;
ksteps = 0;
% ........................EXIT
if( stiff )
%
% PROBLEM APPEARS TO BE STIFF
idid = -4;
stiff = false;
nonstf = false;
ntstep = 0;
% ........................EXIT
nstifs = 0;
end;
gt(2)=1;
break;
end;
break;
end;
if(gt(1)==1)
continue;
end;
if(gt(2)==1)
break;
end;
else;
%
% REQUESTED ERROR UNATTAINABLE DUE TO LIMITED
% PRECISION AVAILABLE
tolfac = 2.0d0.*tolfac;
idid = -2;
% .....................EXIT
gt(2)=1;
break;
end;
break;
end;
if(gt(2)==1)
break;
end;
else;
%
% IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND
% RETURN
%
for k = 1 : neq;
y(k) = y(k) + dt.*yp(k);
end; k = fix(neq+1);
a = tout;
[a,y,yp,rpar,ipar]=df(a,y,yp,rpar,ipar);
ksteps = fix(ksteps + 1);
end;
%
% INTERVAL MODE
idid = 2;
t = tout;
told = t;
% ...............EXIT
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
f1_shape=zeros(f1_shape);f1_shape(:)=f1(1:numel(f1_shape));f1=f1_shape;
f2_shape=zeros(f2_shape);f2_shape(:)=f2(1:numel(f2_shape));f2=f2_shape;
f3_shape=zeros(f3_shape);f3_shape(:)=f3(1:numel(f3_shape));f3=f3_shape;
f4_shape=zeros(f4_shape);f4_shape(:)=f4(1:numel(f4_shape));f4=f4_shape;
f5_shape=zeros(f5_shape);f5_shape(:)=f5(1:numel(f5_shape));f5=f5_shape;
ys_shape=zeros(ys_shape);ys_shape(:)=ys(1:numel(ys_shape));ys=ys_shape;
rtol_shape=zeros(rtol_shape);rtol_shape(:)=rtol(1:numel(rtol_shape));rtol=rtol_shape;
atol_shape=zeros(atol_shape);atol_shape(:)=atol(1:numel(atol_shape));atol=atol_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;
else;
%
% UNNECESSARY FREQUENCY OF OUTPUT IS RESTRICTING
% THE STEP SIZE CHOICE
idid = -5;
kop = 0;
end;
break;
end;
%
% INTEGRATION TASK INTERRUPTED
%
info(1) = -1;
told = t;
% ...EXIT
if( idid==(-2) )
%
% THE ERROR TOLERANCES ARE INCREASED TO VALUES
% WHICH ARE APPROPRIATE FOR CONTINUING
rtol(1) = tolfac.*rtol(1);
atol(1) = tolfac.*atol(1);
% ...EXIT
if( info(2)~=0 )
for k = 2 : neq;
rtol(k) = tolfac.*rtol(k);
atol(k) = tolfac.*atol(k);
end; k = fix(neq+1);
end;
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
f1_shape=zeros(f1_shape);f1_shape(:)=f1(1:numel(f1_shape));f1=f1_shape;
f2_shape=zeros(f2_shape);f2_shape(:)=f2(1:numel(f2_shape));f2=f2_shape;
f3_shape=zeros(f3_shape);f3_shape(:)=f3(1:numel(f3_shape));f3=f3_shape;
f4_shape=zeros(f4_shape);f4_shape(:)=f4(1:numel(f4_shape));f4=f4_shape;
f5_shape=zeros(f5_shape);f5_shape(:)=f5(1:numel(f5_shape));f5=f5_shape;
ys_shape=zeros(ys_shape);ys_shape(:)=ys(1:numel(ys_shape));ys=ys_shape;
rtol_shape=zeros(rtol_shape);rtol_shape(:)=rtol(1:numel(rtol_shape));rtol=rtol_shape;
atol_shape=zeros(atol_shape);atol_shape(:)=atol(1:numel(atol_shape));atol=atol_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 drkfs
%DECK DRLCAL
|
|
