| [f,neq,t,y,tout,rtol,atol,idid,ypout,yh,yh1,ewt,savf,acor,wm,iwm,jac,intout,tstop,tolfac,delsgn,rpar,ipar]=lsod(f,neq,t,y,tout,rtol,atol,idid,ypout,yh,yh1,ewt,savf,acor,wm,iwm,jac,intout,tstop,tolfac,delsgn,rpar,ipar); |
function [f,neq,t,y,tout,rtol,atol,idid,ypout,yh,yh1,ewt,savf,acor,wm,iwm,jac,intout,tstop,tolfac,delsgn,rpar,ipar]=lsod(f,neq,t,y,tout,rtol,atol,idid,ypout,yh,yh1,ewt,savf,acor,wm,iwm,jac,intout,tstop,tolfac,delsgn,rpar,ipar);
persistent absdel big del dt firstCall gt ha intflg k l ltol maxnum natolp nrtolp tol xern1 xern3 xern4 ; if isempty(firstCall),firstCall=1;end;
if isempty(absdel), absdel=0; end;
if isempty(big), big=0; end;
if isempty(del), del=0; end;
if isempty(dt), dt=0; end;
global debdf1_3; if isempty(debdf1_3), debdf1_3=0; end;
global debdf1_4; if isempty(debdf1_4), debdf1_4=0; end;
if isempty(ha), ha=0; end;
global debdf1_5; if isempty(debdf1_5), debdf1_5=0; 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_2; if isempty(debdf1_2), debdf1_2=zeros(1,210); end;
if isempty(tol), tol=0; end;
global debdf1_1; if isempty(debdf1_1), debdf1_1=0; end;
global debdf1_9; if isempty(debdf1_9), debdf1_9=0; end;
global debdf1_8; if isempty(debdf1_8), debdf1_8=0; end;
global debdf1_23; if isempty(debdf1_23), debdf1_23=0; end;
global debdf1_18; if isempty(debdf1_18), debdf1_18=0; end;
global debdf1_25; if isempty(debdf1_25), debdf1_25=0; end;
global debdf1_20; if isempty(debdf1_20), debdf1_20=0; end;
global debdf1_22; if isempty(debdf1_22), debdf1_22=0; end;
global debdf1_11; if isempty(debdf1_11), debdf1_11=0; end;
if isempty(intflg), intflg=0; end;
global debdf1_24; if isempty(debdf1_24), debdf1_24=zeros(1,6); end;
global debdf1_10; if isempty(debdf1_10), debdf1_10=0; end;
global debdf1_19; if isempty(debdf1_19), debdf1_19=0; end;
global debdf1_21; if isempty(debdf1_21), debdf1_21=0; end;
global debdf1_26; if isempty(debdf1_26), debdf1_26=0; end;
if isempty(k), k=0; end;
global debdf1_27; if isempty(debdf1_27), debdf1_27=0; end;
global debdf1_17; if isempty(debdf1_17), debdf1_17=0; end;
if isempty(l), l=0; end;
global debdf1_14; if isempty(debdf1_14), debdf1_14=0; end;
global debdf1_28; if isempty(debdf1_28), debdf1_28=0; end;
global debdf1_13; if isempty(debdf1_13), debdf1_13=0; end;
global debdf1_15; if isempty(debdf1_15), debdf1_15=0; end;
if isempty(ltol), ltol=0; end;
global debdf1_16; if isempty(debdf1_16), debdf1_16=0; end;
global debdf1_12; if isempty(debdf1_12), debdf1_12=0; end;
if isempty(maxnum), maxnum=0; end;
global debdf1_31; if isempty(debdf1_31), debdf1_31=0; end;
global debdf1_29; if isempty(debdf1_29), debdf1_29=0; end;
global debdf1_30; if isempty(debdf1_30), debdf1_30=0; end;
global debdf1_32; if isempty(debdf1_32), debdf1_32=0; end;
if isempty(natolp), natolp=0; end;
global debdf1_35; if isempty(debdf1_35), debdf1_35=0; end;
global debdf1_36; if isempty(debdf1_36), debdf1_36=0; end;
global debdf1_33; if isempty(debdf1_33), debdf1_33=0; end;
global debdf1_37; if isempty(debdf1_37), debdf1_37=0; end;
if isempty(nrtolp), nrtolp=0; end;
global debdf1_34; if isempty(debdf1_34), debdf1_34=0; end;
if isempty(gt), gt=zeros(1,3); end;
%***BEGIN PROLOGUE LSOD
%***SUBSIDIARY
%***PURPOSE Subsidiary to DEBDF
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (LSOD-S, DLSOD-D)
%***AUTHOR (UNKNOWN)
%***DESCRIPTION
%
% DEBDF merely allocates storage for LSOD to relieve the user of
% the inconvenience of a long call list. Consequently LSOD is used
% as described in the comments for DEBDF .
%
%***SEE ALSO DEBDF
%***ROUTINES CALLED HSTART, INTYD, R1MACH, STOD, VNWRMS, XERMSG
%***COMMON BLOCKS DEBDF1
%***REVISION HISTORY (YYMMDD)
% 800901 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 900510 Convert XERRWV calls to XERMSG calls. (RWC)
%***end PROLOGUE LSOD
%
%
%
y_shape=size(y);y=reshape(y,1,[]);
ypout_shape=size(ypout);ypout=reshape(ypout,1,[]);
yh_orig=yh;yh_shape=[neq,6];yh=reshape([yh_orig(1:min(prod(yh_shape),numel(yh_orig))),zeros(1,max(0,prod(yh_shape)-numel(yh_orig)))],yh_shape);
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,[]);
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,[]);
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,16); end;
if isempty(xern4), xern4=repmat(' ',1,16); end;
%
% common :: ;
%% common /debdf1/ told , rowns(210) , el0 , h , hmin , hmxi , hu ,x , u , iquit , init , lyh , lewt , lacor ,lsavf , lwm , ksteps , ibegin , itol , iinteg ,itstop , ijac , iband , iowns(6) , ier , jstart ,kflag , ldum , meth , miter , maxord , n , nq ,nst , nfe , nje , nqu;
%% common /debdf1/ debdf1_1 , debdf1_2(210) , debdf1_3 , debdf1_4 , debdf1_5 , debdf1_6 , debdf1_7 ,debdf1_8 , debdf1_9 , debdf1_10 , debdf1_11 , debdf1_12 , debdf1_13 , debdf1_14 ,debdf1_15 , debdf1_16 , debdf1_17 , debdf1_18 , debdf1_19 , debdf1_20 ,debdf1_21 , debdf1_22 , debdf1_23 , debdf1_24(6) , 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;
%
%
%.......................................................................
%
% THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE
% NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MAXNUM, THE COUNTER
% IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE EXCESSIVE
% WORK.
%
if firstCall, maxnum=[500]; end;
firstCall=0;
%
%.......................................................................
%
%***FIRST EXECUTABLE STATEMENT LSOD
if( debdf1_18==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.
%
[debdf1_9 ]=r1mach(4);
% -- SET ASSOCIATED MACHINE DEPENDENT PARAMETER
wm(1) = sqrt(debdf1_9);
% -- SET TERMINATION FLAG
debdf1_10 = 0;
% -- SET INITIALIZATION INDICATOR
debdf1_11 = 0;
% -- SET COUNTER FOR ATTEMPTED STEPS
debdf1_17 = 0;
% -- SET INDICATOR FOR INTERMEDIATE-OUTPUT
intout = false;
% -- SET START INDICATOR FOR STOD CODE
debdf1_26 = 0;
% -- SET BDF METHOD INDICATOR
debdf1_29 = 2;
% -- SET MAXIMUM ORDER FOR BDF METHOD
debdf1_31 = 5;
% -- SET ITERATION MATRIX INDICATOR
%
if( debdf1_22==0 && debdf1_23==0 )
debdf1_30 = 2;
end;
if( debdf1_22==1 && debdf1_23==0 )
debdf1_30 = 1;
end;
if( debdf1_22==0 && debdf1_23==1 )
debdf1_30 = 5;
end;
if( debdf1_22==1 && debdf1_23==1 )
debdf1_30 = 4;
end;
%
% -- SET OTHER NECESSARY ITEMS IN COMMON BLOCK
debdf1_32 = fix(neq);
debdf1_34 = 0;
debdf1_36 = 0;
debdf1_6 = 0.;
debdf1_33 = 1;
debdf1_4 = 1.;
% -- RESET IBEGIN FOR SUBSEQUENT CALLS
debdf1_18 = 1;
end;
%
%.......................................................................
%
% CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY
%
if( neq<1 )
xern1=sprintf(['%8i'], neq);
xermsg('SLATEC','LSOD',['IN DEBDF, THE NUMBER OF EQUATIONS 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 )
if( rtol(k)<0. )
xern1=sprintf(['%8i'], k);
xern3=sprintf([repmat('%15.6f',1,1)], rtol(k));
xermsg('SLATEC','LSOD',['IN DEBDF, THE RELATIVE ERROR TOLERANCES 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;
if( natolp>0 )
break;
end;
nrtolp = 1;
elseif( natolp>0 ) ;
if( debdf1_19==0 )
break;
end;
continue;
end;
end;
%
if( atol(k)<0. )
xern1=sprintf(['%8i'], k);
xern3=sprintf([repmat('%15.6f',1,1)], atol(k));
xermsg('SLATEC','LSOD',['IN DEBDF, THE ABSOLUTE ERROR ',['TOLERANCES 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;
if( nrtolp>0 )
break;
end;
natolp = 1;
end;
if( debdf1_19==0 )
break;
end;
end;
%
if( debdf1_21==1 )
if( (abs(1.).*sign(tout-t))~=(abs(1.).*sign(tstop-t)) || abs(tout-t)>abs(tstop-t) )
xern3=sprintf([repmat('%15.6f',1,1)], tout);
xern4=sprintf([repmat('%15.6f',1,1)], tstop);
xermsg('SLATEC','LSOD',['IN DEBDF, YOU HAVE CALLED THE ',['CODE WITH TOUT = ',[xern3,['$$BUT YOU HAVE ',['ALSO TOLD THE CODE NOT TO INTEGRATE PAST THE POINT ',['TSTOP = ',[xern4,[' BY SETTING INFO(4) = 1. ','THESE INSTRUCTIONS CONFLICT.']]]]]]]],14,1);
idid = -33;
end;
end;
%
% CHECK SOME CONTINUATION POSSIBILITIES
%
if( debdf1_11~=0 )
if( t==tout )
xern3=sprintf([repmat('%15.6f',1,1)], t);
xermsg('SLATEC','LSOD',['IN DEBDF, YOU HAVE CALLED THE CODE WITH T = TOUT = ',[xern3,' THIS IS NOT ALLOWED ON CONTINUATION CALLS.']],9,1);
idid = -33;
end;
%
if( t~=debdf1_1 )
xern3=sprintf([repmat('%15.6f',1,1)], debdf1_1);
xern4=sprintf([repmat('%15.6f',1,1)], t);
xermsg('SLATEC','LSOD',['IN DEBDF, YOU HAVE CHANGED THE VALUE OF T FROM ',[xern3,[' TO ',[xern4,' THIS IS NOT ALLOWED ON CONTINUATION CALLS.']]]],10,1);
idid = -33;
end;
%
if( debdf1_11~=1 )
if( delsgn.*(tout-t)<0. )
xern3=sprintf([repmat('%15.6f',1,1)], tout);
xermsg('SLATEC','LSOD',['IN DEBDF, 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;
%
if( idid==(-33) )
if( debdf1_10~=(-33) )
% INVALID INPUT DETECTED
debdf1_10 = -33;
debdf1_18 = -1;
else;
xermsg('SLATEC','LSOD',['IN DEBDF, 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);
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
yh_orig(1:prod(yh_shape))=yh;yh=yh_orig;
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;
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;
%
%.......................................................................
%
% 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
% 100*U WHICH IS LIKELY TO BE REASONABLE FOR THIS METHOD AND MACHINE
%
gt(:)=0;
for k = 1 : neq;
if( rtol(k)+atol(k)<=0. )
rtol(k) = 100..*debdf1_9;
idid = -2;
end;
if( debdf1_19==0 )
break;
end;
end;
%
while (1);
if( idid~=(-2) )
%
% BRANCH ON STATUS OF INITIALIZATION INDICATOR
% INIT=0 MEANS INITIAL DERIVATIVES AND NOMINAL STEP SIZE
% AND DIRECTION NOT YET SET
% INIT=1 MEANS NOMINAL STEP SIZE AND DIRECTION NOT YET SET
% INIT=2 MEANS NO FURTHER INITIALIZATION REQUIRED
%
if( debdf1_11==0 )
%
%.......................................................................
%
% MORE INITIALIZATION --
% -- EVALUATE INITIAL DERIVATIVES
%
debdf1_11 = 1;
[t,y,yh(1:end,2:end),rpar,ipar]=f(t,y,yh(1:end,2:end),rpar,ipar);
debdf1_35 = 1;
if( t==tout )
idid = 2;
for l = 1 : neq;
ypout(l) = yh(l,2);
end; l = fix(neq+1);
debdf1_1 = t;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
yh_orig(1:prod(yh_shape))=yh;yh=yh_orig;
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;
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( debdf1_11~=1 ) ;
break;
end;
%
% -- COMPUTE INITIAL STEP SIZE
% -- SAVE SIGN OF INTEGRATION DIRECTION
% -- SET INDEPENDENT AND DEPENDENT VARIABLES
% X AND YH(*) FOR STOD
%
ltol = 1;
for l = 1 : neq;
if( debdf1_19==1 )
ltol = fix(l);
end;
tol = rtol(ltol).*abs(y(l)) + atol(ltol);
if( tol==0. )
gt(3)=1;
break;
end;
ewt(l) = tol;
end;
if(gt(3)==1)
break;
end;
%
big = sqrt(r1mach(2));
[f,neq,t,tout,y,dumvar6,ewt,dumvar8,debdf1_9,big,dumvar11,dumvar12,dumvar13,dumvar14,rpar,ipar,debdf1_4]=hstart(f,neq,t,tout,y,yh(sub2ind(size(yh),1,2):end),ewt,1,debdf1_9,big,yh(sub2ind(size(yh),1,3):end),yh(sub2ind(size(yh),1,4):end),yh(sub2ind(size(yh),1,5):end),yh(sub2ind(size(yh),1,6):end),rpar,ipar,debdf1_4); yh(sub2ind(size(yh),1,2):end)=dumvar6; yh(sub2ind(size(yh),1,3):end)=dumvar11; yh(sub2ind(size(yh),1,4):end)=dumvar12; yh(sub2ind(size(yh),1,5):end)=dumvar13; yh(sub2ind(size(yh),1,6):end)=dumvar14;
%
delsgn = (abs(1.0).*sign(tout-t));
debdf1_8 = t;
for l = 1 : neq;
yh(l,1) = y(l);
yh(l,2) = debdf1_4.*yh(l,2);
end; l = fix(neq+1);
debdf1_11 = 2;
else;
% RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
% SMALL POSITIVE VALUE
debdf1_18 = -1;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
yh_orig(1:prod(yh_shape))=yh;yh=yh_orig;
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;
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;
break;
end;
%
%.......................................................................
%
% ON EACH CALL SET INFORMATION WHICH DETERMINES THE ALLOWED INTERVAL
% OF INTEGRATION BEFORE RETURNING WITH AN ANSWER AT TOUT
%
while (1);
if(gt(3)==0)
if(gt(2)==0)
del = tout - t;
absdel = abs(del);
%
%.......................................................................
%
% IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN
%
end;
gt(2)=0;
if( abs(debdf1_8-t)<absdel )
%
% IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE,
% EXTRAPOLATE AND RETURN
%
if( debdf1_21==1 )
if( abs(tstop-debdf1_8)<100..*debdf1_9.*abs(debdf1_8) )
dt = tout - debdf1_8;
for l = 1 : neq;
y(l) = yh(l,1) +(dt./debdf1_4).*yh(l,2);
end; l = fix(neq+1);
[tout,y,ypout,rpar,ipar]=f(tout,y,ypout,rpar,ipar);
debdf1_35 = fix(debdf1_35 + 1);
idid = 3;
t = tout;
debdf1_1 = t;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
yh_orig(1:prod(yh_shape))=yh;yh=yh_orig;
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;
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;
%
if( ~(debdf1_20==0 || ~intout) )
%
% INTERMEDIATE-OUTPUT MODE
%
idid = 1;
break;
%
%.......................................................................
%
% MONITOR NUMBER OF STEPS ATTEMPTED
%
elseif( debdf1_17<=maxnum ) ;
%
%.......................................................................
%
% LIMIT STEP SIZE AND SET WEIGHT VECTOR
%
debdf1_5 = 100..*debdf1_9.*abs(debdf1_8);
ha = max(abs(debdf1_4),debdf1_5);
if( debdf1_21==1 )
ha = min(ha,abs(tstop-debdf1_8));
end;
debdf1_4 = (abs(ha).*sign(debdf1_4));
ltol = 1;
for l = 1 : neq;
if( debdf1_19==1 )
ltol = fix(l);
end;
ewt(l) = rtol(ltol).*abs(yh(l,1)) + atol(ltol);
if( ewt(l)<=0.0 )
gt(3)=1;
continue;
end;
end; l = fix(neq+1);
tolfac = debdf1_9.*vnwrms(neq,yh,ewt);
if( tolfac<=1. )
%
%.......................................................................
%
% TAKE A STEP
%
neq_orig=neq; [neq,y,yh,dumvar4,yh1,ewt,savf,acor,wm,iwm,f,jac,rpar,ipar]=stod(neq,y,yh,neq,yh1,ewt,savf,acor,wm,iwm,f,jac,rpar,ipar); neq(dumvar4~=neq_orig)=dumvar4(dumvar4~=neq_orig);
%
debdf1_26 = -2;
intout = true;
if( debdf1_27==0 )
gt(2)=1;
continue;
end;
%
%.......................................................................
%
if( debdf1_27==-1 )
%
% REPEATED ERROR TEST FAILURES
idid = -7;
debdf1_18 = -1;
else;
%
% REPEATED CORRECTOR CONVERGENCE FAILURES
idid = -6;
debdf1_18 = -1;
end;
break;
else;
%
% TOLERANCES TOO SMALL
idid = -2;
tolfac = 2..*tolfac;
rtol(1) = tolfac.*rtol(1);
atol(1) = tolfac.*atol(1);
if( debdf1_19~=0 )
for l = 2 : neq;
rtol(l) = tolfac.*rtol(l);
atol(l) = tolfac.*atol(l);
end; l = fix(neq+1);
end;
debdf1_18 = -1;
break;
end;
else;
%
% A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED
idid = -1;
debdf1_17 = 0;
debdf1_18 = -1;
break;
end;
else;
[tout,dumvar2,yh,neq,y,intflg]=intyd(tout,0,yh,neq,y,intflg);
[tout,dumvar2,yh,neq,ypout,intflg]=intyd(tout,1,yh,neq,ypout,intflg);
idid = 3;
if( debdf1_8==tout )
idid = 2;
intout = false;
end;
t = tout;
debdf1_1 = t;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
yh_orig(1:prod(yh_shape))=yh;yh=yh_orig;
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;
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;
%
% RELATIVE ERROR CRITERION INAPPROPRIATE
end;
gt(3)=0;
idid = -3;
debdf1_18 = -1;
break;
end;
%
%.......................................................................
%
% STORE VALUES BEFORE RETURNING TO DEBDF
for l = 1 : neq;
y(l) = yh(l,1);
ypout(l) = yh(l,2)./debdf1_4;
end; l = fix(neq+1);
t = debdf1_8;
debdf1_1 = t;
intout = false;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
yh_orig(1:prod(yh_shape))=yh;yh=yh_orig;
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;
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 lsod
%DECK LSSODS
|
|
