| [neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,f,jac,rpar,ipar]=stod(neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,f,jac,rpar,ipar); |
function [neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,f,jac,rpar,ipar]=stod(neq,y,yh,nyh,yh1,ewt,savf,acor,wm,iwm,f,jac,rpar,ipar);
persistent dcon ddn del delp dsm dup exdn exsm exup gt10 gt11 gt12 gt2 gt3 gt4 gt5 gt6 gt7 gt8 gt9 i i1 iredo iret j jb m ncf newq r rh rhdn rhsm rhup told ;
global debdf1_19; if isempty(debdf1_19), debdf1_19=zeros(1,6); end;
global debdf1_18; if isempty(debdf1_18), debdf1_18=0; end;
%***BEGIN PROLOGUE STOD
%***SUBSIDIARY
%***PURPOSE Subsidiary to DEBDF
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (STOD-S, DSTOD-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% STOD integrates a system of first order odes over one step in the
% integrator package DEBDF.
% ----------------------------------------------------------------------
% STOD performs one step of the integration of an initial value
% problem for a system of ordinary differential equations.
% Note.. STOD is independent of the value of the iteration method
% indicator MITER, when this is .NE. 0, and hence is independent
% of the type of chord method used, or the Jacobian structure.
% Communication with STOD is done with the following variables..
%
% Y = An array of length .GE. debdf1_33 used as the Y argument in
% all calls to F and JAC.
% NEQ = Integer array containing problem size in NEQ(1), and
% passed as the NEQ argument in all calls to F and JAC.
% YH = An NYH by LMAX array containing the dependent variables
% and their approximate scaled derivatives, where
% LMAX = MAXORD + 1. YH(I,J+1) contains the approximate
% J-th derivative of Y(I), scaled by H**J/Factorial(j)
% (J = 0,1,...,NQ). On entry for the first step, the first
% two columns of YH must be set from the initial values.
% NYH = A constant integer .GE. N, the first dimension of YH.
% YH1 = A one-dimensional array occupying the same space as YH.
% EWT = An array of N elements with which the estimated local
% errors in YH are compared.
% SAVF = An array of working storage, of length N.
% ACOR = A work array of length N, used for the accumulated
% corrections. On a successful return, ACOR(I) contains
% the estimated one-step local error in Y(I).
% WM,IWM = Real and integer work arrays associated with matrix
% operations in chord iteration (MITER .NE. 0).
% PJAC = Name of routine to evaluate and preprocess Jacobian matrix
% if a chord method is being used.
% SLVS = Name of routine to solve linear system in chord iteration.
% H = The step size to be attempted on the next step.
% H is altered by the error control algorithm during the
% problem. H can be either positive or negative, but its
% sign must remain constant throughout the problem.
% HMIN = The minimum absolute value of the step size H to be used.
% HMXI = Inverse of the maximum absolute value of H to be used.
% HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
% HMIN and HMXI may be changed at any time, but will not
% take effect until the next change of H is considered.
% TN = The independent variable. TN is updated on each step taken.
% JSTART = An integer used for input only, with the following
% values and meanings..
% 0 Perform the first step.
% .GT.0 Take a new step continuing from the last.
% -1 Take the next step with a new value of H, MAXORD,
% N, METH, MITER, and/or matrix parameters.
% -2 Take the next step with a new value of H,
% but with other inputs unchanged.
% On return, JSTART is set to 1 to facilitate continuation.
% KFLAG = a completion code with the following meanings..
% 0 The step was successful.
% -1 The requested error could not be achieved.
% -2 Corrector convergence could not be achieved.
% A return with KFLAG = -1 or -2 means either
% ABS(H) = HMIN or 10 consecutive failures occurred.
% On a return with KFLAG negative, the values of TN and
% the YH array are as of the beginning of the last
% step, and H is the last step size attempted.
% MAXORD = The maximum order of integration method to be allowed.
% METH/MITER = The method flags. See description in driver.
% N = The number of first-order differential equations.
% ----------------------------------------------------------------------
%
%***SEE ALSO DEBDF
%***ROUTINES CALLED CFOD, PJAC, SLVS, VNWRMS
%***COMMON BLOCKS DEBDF1
%***REVISION HISTORY (YYMMDD)
% 800901 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 910722 Updated AUTHOR section. (ALS)
% 920422 Changed DIMENSION statement. (WRB)
%***end PROLOGUE STOD
%
%LLL. OPTIMIZE
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
global debdf1_20; if isempty(debdf1_20), debdf1_20=0; end;
global debdf1_26; if isempty(debdf1_26), debdf1_26=0; end;
global debdf1_17; if isempty(debdf1_17), debdf1_17=zeros(1,7); end;
if isempty(iredo), iredo=0; end;
if isempty(iret), iret=0; end;
global debdf1_21; if isempty(debdf1_21), debdf1_21=0; end;
if isempty(j), j=0; end;
if isempty(jb), jb=0; end;
global debdf1_27; if isempty(debdf1_27), debdf1_27=0; end;
global debdf1_28; if isempty(debdf1_28), debdf1_28=0; end;
global debdf1_29; if isempty(debdf1_29), debdf1_29=0; end;
global debdf1_22; if isempty(debdf1_22), debdf1_22=0; end;
if isempty(m), m=0; end;
global debdf1_32; if isempty(debdf1_32), debdf1_32=0; end;
global debdf1_23; if isempty(debdf1_23), debdf1_23=0; end;
global debdf1_30; if isempty(debdf1_30), debdf1_30=0; end;
global debdf1_31; if isempty(debdf1_31), debdf1_31=0; end;
global debdf1_33; if isempty(debdf1_33), debdf1_33=0; end;
if isempty(ncf), ncf=0; end;
if isempty(newq), newq=0; end;
global debdf1_36; if isempty(debdf1_36), debdf1_36=0; end;
global debdf1_37; if isempty(debdf1_37), debdf1_37=0; end;
global debdf1_34; if isempty(debdf1_34), debdf1_34=0; end;
global debdf1_24; if isempty(debdf1_24), debdf1_24=0; end;
global debdf1_38; if isempty(debdf1_38), debdf1_38=0; end;
global debdf1_35; if isempty(debdf1_35), debdf1_35=0; end;
global debdf1_25; if isempty(debdf1_25), debdf1_25=0; end;
if isempty(gt2), gt2=0; end;
if isempty(gt3), gt3=0; end;
if isempty(gt4), gt4=0; end;
if isempty(gt5), gt5=0; end;
if isempty(gt6), gt6=0; end;
if isempty(gt7), gt7=0; end;
if isempty(gt8), gt8=0; end;
if isempty(gt9), gt9=0; end;
if isempty(gt10), gt10=0; end;
if isempty(gt11), gt11=0; end;
if isempty(gt12), gt12=0; end;
global debdf1_1; if isempty(debdf1_1), debdf1_1=0; end;
global debdf1_2; if isempty(debdf1_2), debdf1_2=0; end;
global debdf1_3; if isempty(debdf1_3), debdf1_3=0; end;
global debdf1_4; if isempty(debdf1_4), debdf1_4=zeros(1,13); end;
global debdf1_5; if isempty(debdf1_5), debdf1_5=zeros(13,12); end;
global debdf1_6; if isempty(debdf1_6), debdf1_6=0; end;
global debdf1_7; if isempty(debdf1_7), debdf1_7=0; end;
global debdf1_8; if isempty(debdf1_8), debdf1_8=0; end;
global debdf1_9; if isempty(debdf1_9), debdf1_9=zeros(3,12); end;
global debdf1_10; if isempty(debdf1_10), debdf1_10=0; end;
global debdf1_11; if isempty(debdf1_11), debdf1_11=0; end;
global debdf1_12; if isempty(debdf1_12), debdf1_12=0; end;
global debdf1_13; if isempty(debdf1_13), debdf1_13=0; end;
global debdf1_14; if isempty(debdf1_14), debdf1_14=0; end;
global debdf1_15; if isempty(debdf1_15), debdf1_15=0; end;
global debdf1_16; if isempty(debdf1_16), debdf1_16=0; end;
if isempty(dcon), dcon=0; end;
if isempty(ddn), ddn=0; end;
if isempty(del), del=0; end;
if isempty(delp), delp=0; end;
if isempty(dsm), dsm=0; end;
if isempty(dup), dup=0; end;
if isempty(exdn), exdn=0; end;
if isempty(exsm), exsm=0; end;
if isempty(exup), exup=0; end;
if isempty(r), r=0; end;
if isempty(rh), rh=0; end;
if isempty(rhdn), rhdn=0; end;
if isempty(rhsm), rhsm=0; end;
if isempty(rhup), rhup=0; end;
if isempty(told), told=0; end;
y_shape=size(y);y=reshape(y,1,[]);
yh_shape=size(yh);yh=reshape([yh(:).',zeros(1,ceil(numel(yh)./prod([nyh])).*prod([nyh])-numel(yh))],nyh,[]);
yh1_shape=size(yh1);yh1=reshape(yh1,1,[]);
ewt_shape=size(ewt);ewt=reshape(ewt,1,[]);
savf_shape=size(savf);savf=reshape(savf,1,[]);
acor_shape=size(acor);acor=reshape(acor,1,[]);
wm_shape=size(wm);wm=reshape(wm,1,[]);
iwm_shape=size(iwm);iwm=reshape(iwm,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
% common :: ;
%% common /debdf1/ rownd,conit,crate,el(13),elco(13,12),hold,rc,rmax,tesco(3,12),el0,h,hmin,hmxi,hu,tn,uround,iownd(7),ksteps,iod(6),ialth,ipup,lmax,meo,nqnyh,nstepj,ier,jstart,kflag,l,meth,miter ,maxord,n,nq,nst,nfe,nje,nqu;
%% common /debdf1/ debdf1_1,debdf1_2,debdf1_3,debdf1_4(13),debdf1_5(13,12),debdf1_6,debdf1_7,debdf1_8,debdf1_9(3,12),debdf1_10,debdf1_11,debdf1_12,debdf1_13,debdf1_14,debdf1_15,debdf1_16,debdf1_17(7),debdf1_18,debdf1_19(6),debdf1_20,debdf1_21,debdf1_22,debdf1_23,debdf1_24,debdf1_25,debdf1_26,debdf1_27,debdf1_28,debdf1_29,debdf1_30,debdf1_31 ,debdf1_32,debdf1_33,debdf1_34,debdf1_35,debdf1_36,debdf1_37,debdf1_38;
%
%
%***FIRST EXECUTABLE STATEMENT STOD
debdf1_28 = 0;
told = debdf1_15;
ncf = 0;
gt2=0;
gt3=0;
gt4=0;
gt5=0;
gt6=0;
gt7=0;
gt8=0;
gt9=0;
gt10=0;
gt11=0;
gt12=0;
while (1);
if( debdf1_27>0 )
gt4=1;
break;
end;
if( debdf1_27==-1 )
%-----------------------------------------------------------------------
% THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
% IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
% IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
% IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
% IF THE CALLER HAS CHANGED METH, CFOD IS CALLED TO RESET
% THE COEFFICIENTS OF THE METHOD.
% IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
% ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
% IF H IS TO BE CHANGED, YH MUST BE RESCALED.
% IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
% TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
%-----------------------------------------------------------------------
debdf1_21 = fix(debdf1_31);
debdf1_22 = fix(debdf1_32 + 1);
if( debdf1_20==1 )
debdf1_20 = 2;
end;
if( debdf1_30~=debdf1_23 )
[debdf1_30,debdf1_5,debdf1_9]=cfod(debdf1_30,debdf1_5,debdf1_9);
debdf1_23 = fix(debdf1_30);
if( debdf1_34<=debdf1_32 )
debdf1_20 = fix(debdf1_29);
iret = 1;
break;
end;
elseif( debdf1_34<=debdf1_32 ) ;
gt2=1;
break;
end;
debdf1_34 = fix(debdf1_32);
debdf1_29 = fix(debdf1_22);
for i = 1 : debdf1_29;
debdf1_4(i) = debdf1_5(i,debdf1_34);
end; i = fix(debdf1_29+1);
debdf1_24 = fix(debdf1_34.*nyh);
debdf1_7 = debdf1_7.*debdf1_4(1)./debdf1_10;
debdf1_10 = debdf1_4(1);
debdf1_2 = 0.5e0./(debdf1_34+2);
ddn = vnwrms(debdf1_33,savf,ewt)./debdf1_9(1,debdf1_29);
exdn = 1.0e0./debdf1_29;
rhdn = 1.0e0./(1.3e0.*ddn.^exdn+0.0000013e0);
rh = min(rhdn,1.0e0);
iredo = 3;
if( debdf1_11==debdf1_6 )
rh = max(rh,debdf1_12./abs(debdf1_11));
else;
rh = min(rh,abs(debdf1_11./debdf1_6));
debdf1_11 = debdf1_6;
end;
gt3=1;
break;
else;
if( debdf1_27==-2 )
gt2=1;
break;
end;
%-----------------------------------------------------------------------
% ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
% INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
% IN A SINGLE STEP. IT IS INITIALLY 1.0E4 TO COMPENSATE FOR THE SMALL
% INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
% OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
% FOR THE NEXT INCREASE.
%-----------------------------------------------------------------------
debdf1_22 = fix(debdf1_32 + 1);
debdf1_34 = 1;
debdf1_29 = 2;
debdf1_20 = 2;
debdf1_8 = 10000.0e0;
debdf1_7 = 0.0e0;
debdf1_10 = 1.0e0;
debdf1_3 = 0.7e0;
delp = 0.0e0;
debdf1_6 = debdf1_11;
debdf1_23 = fix(debdf1_30);
debdf1_25 = 0;
iret = 3;
%-----------------------------------------------------------------------
% CFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
% CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
% WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
%-----------------------------------------------------------------------
[debdf1_30,debdf1_5,debdf1_9]=cfod(debdf1_30,debdf1_5,debdf1_9);
end;
break;
end;
while (1);
if(gt12==0)
if(gt11==0)
if(gt10==0)
if(gt9==0)
if(gt8==0)
if(gt7==0)
if(gt6==0)
if(gt5==0)
if(gt4==0)
if(gt3==0)
if(gt2==0)
for i = 1 : debdf1_29;
debdf1_4(i) = debdf1_5(i,debdf1_34);
end; i = fix(debdf1_29+1);
debdf1_24 = fix(debdf1_34.*nyh);
debdf1_7 = debdf1_7.*debdf1_4(1)./debdf1_10;
debdf1_10 = debdf1_4(1);
debdf1_2 = 0.5e0./(debdf1_34+2);
if( iret==2 )
rh = max(rh,debdf1_12./abs(debdf1_11));
gt3 =1;
continue;
elseif( iret==3 ) ;
gt4 =1;
continue;
end;
%-----------------------------------------------------------------------
% IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
% RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO
% L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
% FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
%-----------------------------------------------------------------------
%gt2
end;
gt2=0;
if( debdf1_11==debdf1_6 )
gt4=1;
continue;
end;
rh = debdf1_11./debdf1_6;
debdf1_11 = debdf1_6;
iredo = 3;
%gt3
end;
gt3=0;
rh = min(rh,debdf1_8);
rh = rh./max(1.0e0,abs(debdf1_11).*debdf1_13.*rh);
r = 1.0e0;
for j = 2 : debdf1_29;
r = r.*rh;
for i = 1 : debdf1_33;
yh(i,j) = yh(i,j).*r;
end; i = fix(debdf1_33+1);
end; j = fix(debdf1_29+1);
debdf1_11 = debdf1_11.*rh;
debdf1_7 = debdf1_7.*rh;
debdf1_20 = fix(debdf1_29);
if( iredo==0 )
debdf1_8 = 10.0e0;
gt12 =1;
continue;
end;
%-----------------------------------------------------------------------
% THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
% MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
% RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
% WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER
% TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
% IN ANY CASE, PJAC IS CALLED AT LEAST EVERY 20-TH STEP.
%-----------------------------------------------------------------------
%gt4
end;
gt4=0;
if( abs(debdf1_7-1.0e0)>0.3e0 )
debdf1_21 = fix(debdf1_31);
end;
if( debdf1_35>=debdf1_25+20 )
debdf1_21 = fix(debdf1_31);
end;
debdf1_15 = debdf1_15 + debdf1_11;
i1 = fix(debdf1_24 + 1);
for jb = 1 : debdf1_34;
i1 = fix(i1 - nyh);
for i = i1 : debdf1_24;
yh1(i) = yh1(i) + yh1(i+nyh);
end; i = fix(debdf1_24+1);
end; jb = fix(debdf1_34+1);
debdf1_18 = fix(debdf1_18 + 1);
%-----------------------------------------------------------------------
% UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS
% MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR
% WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE
% VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
%-----------------------------------------------------------------------
%gt5
end;
gt5=0;
m = 0;
for i = 1 : debdf1_33;
y(i) = yh(i,1);
end; i = fix(debdf1_33+1);
[debdf1_15,y,savf,rpar,ipar]=f(debdf1_15,y,savf,rpar,ipar);
debdf1_36 = fix(debdf1_36 + 1);
if( debdf1_21>0 )
%-----------------------------------------------------------------------
% IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND
% PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET
% TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
%-----------------------------------------------------------------------
debdf1_21 = 0;
debdf1_7 = 1.0e0;
debdf1_25 = fix(debdf1_35);
debdf1_3 = 0.7e0;
[neq,y,yh,nyh,ewt,acor,savf,wm,iwm,f,jac,rpar,ipar]=pjac(neq,y,yh,nyh,ewt,acor,savf,wm,iwm,f,jac,rpar,ipar);
if( debdf1_26~=0 )
gt8=1;
continue;
end;
end;
for i = 1 : debdf1_33;
acor(i) = 0.0e0;
end; i = fix(debdf1_33+1);
%gt6
end;
gt6=0;
if( debdf1_31~=0 )
%-----------------------------------------------------------------------
% IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
% AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
% P AS COEFFICIENT MATRIX.
%-----------------------------------------------------------------------
for i = 1 : debdf1_33;
y(i) = debdf1_11.*savf(i) -(yh(i,2)+acor(i));
end; i = fix(debdf1_33+1);
[wm,iwm,y,savf]=slvs(wm,iwm,y,savf);
if( debdf1_26~=0 )
gt7=1;
continue;
end;
[del ,debdf1_33,y,ewt]=vnwrms(debdf1_33,y,ewt);
for i = 1 : debdf1_33;
acor(i) = acor(i) + y(i);
y(i) = yh(i,1) + debdf1_4(1).*acor(i);
end; i = fix(debdf1_33+1);
else;
%-----------------------------------------------------------------------
% IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
% THE RESULT OF THE LAST FUNCTION EVALUATION.
%-----------------------------------------------------------------------
for i = 1 : debdf1_33;
savf(i) = debdf1_11.*savf(i) - yh(i,2);
y(i) = savf(i) - acor(i);
end; i = fix(debdf1_33+1);
[del ,debdf1_33,y,ewt]=vnwrms(debdf1_33,y,ewt);
for i = 1 : debdf1_33;
y(i) = yh(i,1) + debdf1_4(1).*savf(i);
acor(i) = savf(i);
end; i = fix(debdf1_33+1);
end;
%-----------------------------------------------------------------------
% TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
% RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
%-----------------------------------------------------------------------
if( m~=0 )
debdf1_3 = max(0.2e0.*debdf1_3,del./delp);
end;
dcon = del.*min(1.0e0,1.5e0.*debdf1_3)./(debdf1_9(2,debdf1_34).*debdf1_2);
if( dcon<=1.0e0 )
%-----------------------------------------------------------------------
% THE CORRECTOR HAS CONVERGED. IPUP IS SET TO -1 IF MITER .NE. 0,
% TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER.
% THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
% IF IT FAILS.
%-----------------------------------------------------------------------
if( debdf1_31~=0 )
debdf1_21 = -1;
end;
if( m==0 )
dsm = del./debdf1_9(2,debdf1_34);
end;
if( m>0 )
dsm = vnwrms(debdf1_33,acor,ewt)./debdf1_9(2,debdf1_34);
end;
if( dsm>1.0e0 )
%-----------------------------------------------------------------------
% THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
% RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
% TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
% ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
% BY A FACTOR OF 0.2 OR LESS.
%-----------------------------------------------------------------------
debdf1_28 = fix(debdf1_28 - 1);
debdf1_15 = told;
i1 = fix(debdf1_24 + 1);
for jb = 1 : debdf1_34;
i1 = fix(i1 - nyh);
for i = i1 : debdf1_24;
yh1(i) = yh1(i) - yh1(i+nyh);
end; i = fix(debdf1_24+1);
end; jb = fix(debdf1_34+1);
debdf1_8 = 2.0e0;
if( abs(debdf1_11)<=debdf1_12.*1.00001e0 )
%-----------------------------------------------------------------------
% ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
% TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
%-----------------------------------------------------------------------
debdf1_28 = -1;
break;
elseif( debdf1_28<=-3 ) ;
%-----------------------------------------------------------------------
% CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURRED.
% IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
% IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
% YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
% DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
% H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
% UNTIL IT SUCCEEDS OR H REACHES HMIN.
%-----------------------------------------------------------------------
if( debdf1_28==-10 )
debdf1_28 = -1;
break;
else;
rh = 0.1e0;
rh = max(debdf1_12./abs(debdf1_11),rh);
debdf1_11 = debdf1_11.*rh;
for i = 1 : debdf1_33;
y(i) = yh(i,1);
end; i = fix(debdf1_33+1);
[debdf1_15,y,savf,rpar,ipar]=f(debdf1_15,y,savf,rpar,ipar);
debdf1_36 = fix(debdf1_36 + 1);
for i = 1 : debdf1_33;
yh(i,2) = debdf1_11.*savf(i);
end; i = fix(debdf1_33+1);
debdf1_21 = fix(debdf1_31);
debdf1_20 = 5;
if( debdf1_34==1 )
gt4=1;
continue;
end;
debdf1_34 = 1;
debdf1_29 = 2;
iret = 3;
continue;
end;
else;
iredo = 2;
rhup = 0.0e0;
gt9 =1;
continue;
end;
else;
%-----------------------------------------------------------------------
% AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
% CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1.
% IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
% use IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
% IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
% BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
% FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT
% TESTING FOR THAT MANY STEPS.
%-----------------------------------------------------------------------
debdf1_28 = 0;
iredo = 0;
debdf1_35 = fix(debdf1_35 + 1);
debdf1_14 = debdf1_11;
debdf1_38 = fix(debdf1_34);
for j = 1 : debdf1_29;
for i = 1 : debdf1_33;
yh(i,j) = yh(i,j) + debdf1_4(j).*acor(i);
end; i = fix(debdf1_33+1);
end; j = fix(debdf1_29+1);
debdf1_20 = fix(debdf1_20 - 1);
if( debdf1_20==0 )
%-----------------------------------------------------------------------
% REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
% RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
% AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
% IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
% THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
% ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
% ADDITIONAL SCALED DERIVATIVE.
%-----------------------------------------------------------------------
rhup = 0.0e0;
if( debdf1_29~=debdf1_22 )
for i = 1 : debdf1_33;
savf(i) = acor(i) - yh(i,debdf1_22);
end; i = fix(debdf1_33+1);
dup = vnwrms(debdf1_33,savf,ewt)./debdf1_9(3,debdf1_34);
exup = 1.0e0./(debdf1_29+1);
rhup = 1.0e0./(1.4e0.*dup.^exup+0.0000014e0);
end;
gt9 =1;
continue;
else;
if( debdf1_20<=1 )
if( debdf1_29~=debdf1_22 )
for i = 1 : debdf1_33;
yh(i,debdf1_22) = acor(i);
end; i = fix(debdf1_33+1);
end;
end;
gt12 =1;
continue;
end;
end;
else;
m = fix(m + 1);
if( m~=3 )
if( m<2 || del<=2.0e0.*delp )
delp = del;
[debdf1_15,y,savf,rpar,ipar]=f(debdf1_15,y,savf,rpar,ipar);
debdf1_36 = fix(debdf1_36 + 1);
gt6 =1;
continue;
end;
end;
end;
%-----------------------------------------------------------------------
% THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES.
% IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
% THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
% BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE
% REDUCED OR 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
%-----------------------------------------------------------------------
%gt7
end;
gt7=0;
if( debdf1_21~=0 )
debdf1_21 = fix(debdf1_31);
gt5 =1;
continue;
end;
%gt8
end;
gt8=0;
debdf1_15 = told;
ncf = fix(ncf + 1);
debdf1_8 = 2.0e0;
i1 = fix(debdf1_24 + 1);
for jb = 1 : debdf1_34;
i1 = fix(i1 - nyh);
for i = i1 : debdf1_24;
yh1(i) = yh1(i) - yh1(i+nyh);
end; i = fix(debdf1_24+1);
end; jb = fix(debdf1_34+1);
if( abs(debdf1_11)<=debdf1_12.*1.00001e0 )
debdf1_28 = -2;
break;
elseif( ncf==10 ) ;
debdf1_28 = -2;
break;
else;
rh = 0.25e0;
debdf1_21 = fix(debdf1_31);
iredo = 1;
rh = max(rh,debdf1_12./abs(debdf1_11));
gt3 =1;
continue;
end;
%gt9
end;
gt9=0;
exsm = 1.0e0./debdf1_29;
rhsm = 1.0e0./(1.2e0.*dsm.^exsm+0.0000012e0);
rhdn = 0.0e0;
if( debdf1_34~=1 )
ddn = vnwrms(debdf1_33,yh(sub2ind(size(yh),1,debdf1_29):end),ewt)./debdf1_9(1,debdf1_34);
exdn = 1.0e0./debdf1_34;
rhdn = 1.0e0./(1.3e0.*ddn.^exdn+0.0000013e0);
end;
if( rhsm>=rhup )
if( rhsm>=rhdn )
newq = fix(debdf1_34);
rh = rhsm;
gt10 =1;
continue;
end;
elseif( rhup>rhdn ) ;
newq = fix(debdf1_29);
rh = rhup;
if( rh<1.1e0 )
debdf1_20 = 3;
gt12 =1;
continue;
else;
r = debdf1_4(debdf1_29)./debdf1_29;
for i = 1 : debdf1_33;
yh(i,newq+1) = acor(i).*r;
end; i = fix(debdf1_33+1);
gt11 =1;
continue;
end;
end;
newq = fix(debdf1_34 - 1);
rh = rhdn;
if( debdf1_28<0 && rh>1.0e0 )
rh = 1.0e0;
end;
%gt10
end;
gt10=0;
if((debdf1_28==0) &&(rh<1.1e0) )
debdf1_20 = 3;
gt12 =1;
continue;
else;
if( debdf1_28<=-2 )
rh = min(rh,0.2e0);
end;
%-----------------------------------------------------------------------
% IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
% IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
% THEN EXIT FROM 680 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
%-----------------------------------------------------------------------
if( newq==debdf1_34 )
rh = max(rh,debdf1_12./abs(debdf1_11));
gt3 =1;
continue;
end;
end;
%gt11
end;
gt11=0;
debdf1_34 = fix(newq);
debdf1_29 = fix(debdf1_34 + 1);
iret = 2;
continue;
%gt12
end;
gt12=0;
r = 1.0e0./debdf1_9(2,debdf1_38);
for i = 1 : debdf1_33;
acor(i) = acor(i).*r;
end; i = fix(debdf1_33+1);
break;
end;
debdf1_6 = debdf1_11;
debdf1_27 = 1;
%----------------------- end OF SUBROUTINE STOD -----------------------
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yh_shape=zeros(yh_shape);yh_shape(:)=yh(1:numel(yh_shape));yh=yh_shape;
yh1_shape=zeros(yh1_shape);yh1_shape(:)=yh1(1:numel(yh1_shape));yh1=yh1_shape;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
savf_shape=zeros(savf_shape);savf_shape(:)=savf(1:numel(savf_shape));savf=savf_shape;
acor_shape=zeros(acor_shape);acor_shape(:)=acor(1:numel(acor_shape));acor=acor_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_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 stod
%DECK STOR1
|
|
