| [f,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]=derkfs(f,neq,t,y,tout,info,rtol,atol,idid,h,tolfac,yp,f1,f2,f3,f4,f5,ys,told,dtsign,u26,rer,init,ksteps,kop,i |
function [f,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]=derkfs(f,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);
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(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(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,4); 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;
%***BEGIN PROLOGUE DERKFS
%***SUBSIDIARY
%***PURPOSE Subsidiary to DERKF
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (DERKFS-S, DRKFS-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% Fehlberg Fourth-Fifth order Runge-Kutta Method
% **********************************************************************
%
% DERKFS integrates a system of first order ordinary differential
% equations as described in the comments for DERKF .
%
% 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 DERKF
%***ROUTINES CALLED DEFEHL, HSTART, HVNRM, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800501 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 891024 Changed references from VNORM to HVNRM. (WRB)
% 891024 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 900510 Convert XERRWV calls to XERMSG calls, replace GO TOs with
% IF-THEN-ELSEs. (RWC)
% 910722 Updated AUTHOR section. (ALS)
%***end PROLOGUE DERKFS
%
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.0e-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 DERKFS
gt(:)=0;
if( info(1)==0 )
%
% ON THE FIRST CALL , PERFORM INITIALIZATION --
% DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING THE
% function ROUTINE R1MACH. THE USER MUST MAKE SURE THAT THE
% VALUES SET IN R1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
%
[u ]=r1mach(4);
% -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS
u26 = 26..*u;
rer = 2..*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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, INFO(3) MUST BE 0 OR 1 INDICATING THE ',['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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, 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','DERKFS',['IN DERKF, 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);
gt(:)=0;
if( idid~=(-33) )
%
%.......................................................................
%
% 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. )
rtol(k) = rer;
idid = -2;
end;
if( info(2)==0 )
break;
end;
end;
%
if( idid==(-2) )
%
% RTOL=ATOL=0 ON INPUT, SO RTOL WAS CHANGED TO A
% SMALL POSITIVE VALUE
tolfac = 1.;
else;
%
% 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]=f(a,y,yp,rpar,ipar);
if( t==tout )
idid = 2;
t = tout;
told = t;
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;
elseif( init~=1 ) ;
gt(1)=1;
end;
%
% -- SET SIGN OF INTEGRATION DIRECTION AND
% -- ESTIMATE STARTING STEP SIZE
%
if(gt(1)==0)
init = 2;
dtsign = (abs(1.).*sign(tout-t));
[u ]=r1mach(4);
big = sqrt(r1mach(2));
ute = u.^0.375;
dy = ute.*hvnrm(y,neq);
if( dy==0. )
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. )
tol = dy.*rtol(ktol);
end;
f1(k) = tol;
end; k = fix(neq+1);
%
[f,neq,t,tout,y,yp,f1,dumvar8,u,big,f2,f3,f4,f5,rpar,ipar,h]=hstart(f,neq,t,tout,y,yp,f1,4,u,big,f2,f3,f4,f5,rpar,ipar,h);
end;
gt(1)=0;
%
%.......................................................................
%
% SET STEP SIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT
% AND SET OUTPUT POINT INDICATOR
%
dt = tout - t;
h = (abs(h).*sign(dt));
output = false;
%
% TEST TO SEE IF DERKF IS BEING SEVERELY IMPACTED BY TOO MANY
% OUTPUT POINTS
%
if( abs(h)>=2..*abs(dt) )
kop = fix(kop + 1);
end;
if( kop>mxkop )
%
% UNNECESSARY FREQUENCY OF OUTPUT IS RESTRICTING
% THE STEP SIZE CHOICE
idid = -5;
kop = 0;
gt(4)=1;
break;
elseif( abs(dt)>u26.*abs(t) ) ;
%
% **********************************************************************
% **********************************************************************
% STEP BY STEP INTEGRATION
%
while( true );
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.;
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. )
tolfac = max(tolfac,abs(y(k)).*(rer./et));
else;
tolfac = max(tolfac,rer./rtol(ktol));
end;
end; k = fix(neq+1);
if( tolfac>1. )
break;
end;
%
% 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..*abs(h) )
if( abs(dt)>abs(h)./0.9 )
%
h = 0.5.*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( ksteps<=mxstep );
%
% ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H
%
[f,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar]=defehl(f,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.;
estiff = 0.;
ktol = 1;
for k = 1 : neq;
yavg = 0.5.*(abs(y(k))+abs(ys(k)));
if( info(2)==1 )
ktol = fix(k);
end;
et = rtol(ktol).*yavg + atol(ktol);
if( et<=0. )
%
% PURE RELATIVE ERROR INAPPROPRIATE WHEN SOLUTION
% VANISHES
idid = -3;
gt(4)=1;
break;
else;
%
ee = abs((-2090..*yp(k)+(21970..*f3(k)-15048..*f4(k)))+(22528..*f2(k)-27360..*f5(k)));
if( ~(stiff || nonstf) )
es = abs(h.*(0.055455.*yp(k)-0.035493.*f1(k)-0.036571.*f2(k)+0.023107.*f3(k)-0.009515.*f4(k)+0.003017.*f5(k)));
estiff = max(estiff,es./et);
end;
eeoet = max(eeoet,ee./et);
end;
end;
if(gt(4)==1)
break;
end;
%
esttol = abs(h).*eeoet./752400.;
%
if( esttol<=1. )
%
%.......................................................................
%
% 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]=f(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.;
if( esttol>1.889568e-4 )
s = 0.9./esttol.^0.2;
end;
if( hfaild )
s = min(s,1.);
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.
%
if( ~(stiff) )
ntstep = fix(rem(ntstep+1,50));
if( ntstep==1 )
nonstf = false;
end;
if( ~(nonstf) )
if( estiff<=1. )
%
% 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
%
if( output )
idid = 2;
t = tout;
told = t;
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
%
if( info(3)==0 )
gt(2)=1;
break;
end;
%
% **********************************************************************
% **********************************************************************
%
% INTEGRATION SUCCESSFULLY COMPLETED
%
% ONE-STEP MODE
idid = 1;
told = t;
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;
elseif( abs(h)<=hmin ) ;
%
% REQUESTED ERROR UNATTAINABLE AT SMALLEST
% ALLOWABLE STEP SIZE
tolfac = 1.69.*esttol;
idid = -2;
gt(4)=1;
break;
else;
%
% REDUCE THE STEP SIZE , TRY AGAIN
% THE DECREASE IS LIMITED TO A FACTOR OF 1/10
%
hfaild = true;
output = false;
s = 0.1;
if( esttol<59049. )
s = 0.9./esttol.^0.2;
end;
h = (abs(max(s.*abs(h),hmin)).*sign(h));
end;
end;
if(gt(2)==1)
gt(2)=0;
continue;
end;
if(gt(4)==1)
break;
end;
%
% A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED
idid = -1;
ksteps = 0;
if( stiff )
%
% PROBLEM APPEARS TO BE STIFF
idid = -4;
stiff = false;
nonstf = false;
ntstep = 0;
nstifs = 0;
end;
gt(4)=1;
break;
end;
if(gt(4)==1)
break;
end;
%
% REQUESTED ERROR UNATTAINABLE DUE TO LIMITED
% PRECISION AVAILABLE
tolfac = 2..*tolfac;
idid = -2;
%to100
break;
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]=f(a,y,yp,rpar,ipar);
ksteps = fix(ksteps + 1);
end;
%
% INTERVAL MODE
idid = 2;
t = tout;
told = t;
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;
elseif( iquit==(-33) ) ;
xermsg('SLATEC','DERKFS',['IN DERKF, 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;
else;
iquit = -33;
end;
break;
end;
%
% INTEGRATION TASK INTERRUPTED
%
info(1) = -1;
told = t;
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);
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 derkfs
%DECK DES
|
|
