| [df,neq,t,y,tout,info,rtol,atol,idid,ypout,yp,yy,wt,p,phi,alpha,beta,psi,v,w,sig,g,gi,h,eps,x,xold,hold,told,delsgn,tstop,twou,fouru,start,phase1,nornd,stiff,intout,ns,kord,kold,init,ksteps,kle4,iquit,kprev,ivc,iv,kgi,rpar,ipar]=ddes(df,neq,t,y,tout,info, |
function [df,neq,t,y,tout,info,rtol,atol,idid,ypout,yp,yy,wt,p,phi,alpha,beta,psi,v,w,sig,g,gi,h,eps,x,xold,hold,told,delsgn,tstop,twou,fouru,start,phase1,nornd,stiff,intout,ns,kord,kold,init,ksteps,kle4,iquit,kprev,ivc,iv,kgi,rpar,ipar]=ddes(df,neq,t,y,tout,info,rtol,atol,idid,ypout,yp,yy,wt,p,phi,alpha,beta,psi,v,w,sig,g,gi,h,eps,x,xold,hold,told,delsgn,tstop,twou,fouru,start,phase1,nornd,stiff,intout,ns,kord,kold,init,ksteps,kle4,iquit,kprev,ivc,iv,kgi,rpar,ipar);
%***BEGIN PROLOGUE DDES
%***SUBSIDIARY
%***PURPOSE Subsidiary to DDEABM
%***LIBRARY SLATEC
%***TYPE doubleprecision (DES-S, DDES-D)
%***AUTHOR Watts, H. A., (SNLA)
%***DESCRIPTION
%
% DDEABM merely allocates storage for DDES to relieve the user of the
% inconvenience of a long call list. Consequently DDES is used as
% described in the comments for DDEABM .
%
%***SEE ALSO DDEABM
%***ROUTINES CALLED D1MACH, DINTP, DSTEPS, XERMSG
%***REVISION HISTORY (YYMMDD)
% 820301 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, cvt GO TOs to
% IF-THEN-ELSE. (RWC)
% 910722 Updated AUTHOR section. (ALS)
%***end PROLOGUE DDES
%
persistent a absdel crash del dt firstCall gt ha k l ll ltol maxnum natolp nrtolp u xern1 xern3 xern4 ; if isempty(firstCall),firstCall=1;end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(ltol), ltol=0; end;
if isempty(maxnum), maxnum=0; end;
if isempty(natolp), natolp=0; end;
if isempty(nrtolp), nrtolp=0; end;
if isempty(ll), ll=0; end;
if isempty(gt), gt=0; end;
if isempty(a), a=0; end;
if isempty(absdel), absdel=0; end;
if isempty(del), del=0; end;
if isempty(dt), dt=0; end;
if isempty(ha), ha=0; end;
if isempty(u), u=0; end;
if isempty(crash), crash=false; end;
%
y_shape=size(y);y=reshape(y,1,[]);
yy_shape=size(yy);yy=reshape(yy,1,[]);
wt_shape=size(wt);wt=reshape(wt,1,[]);
phi_orig=phi;phi_shape=[neq,16];phi=reshape([phi_orig(1:min(prod(phi_shape),numel(phi_orig))),zeros(1,max(0,prod(phi_shape)-numel(phi_orig)))],phi_shape);
p_shape=size(p);p=reshape(p,1,[]);
yp_shape=size(yp);yp=reshape(yp,1,[]);
ypout_shape=size(ypout);ypout=reshape(ypout,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;
%
%
%.......................................................................
%
% 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 DDES
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
twou = 2.0d0.*u;
fouru = 4.0d0.*u;
% -- SET TERMINATION FLAG
iquit = 0;
% -- SET INITIALIZATION INDICATOR
init = 0;
% -- SET COUNTER FOR ATTEMPTED STEPS
ksteps = 0;
% -- SET INDICATOR FOR INTERMEDIATE-OUTPUT
intout = false;
% -- SET INDICATOR FOR STIFFNESS DETECTION
stiff = false;
% -- SET STEP COUNTER FOR STIFFNESS DETECTION
kle4 = 0;
% -- SET INDICATORS FOR STEPS CODE
start = true;
phase1 = true;
nornd = true;
% -- 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','DDES',['IN DDEABM, 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','DDES',['IN DDEABM, 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','DDES',['IN DDEABM, 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( info(4)~=0 && info(4)~=1 )
xern1=sprintf(['%8i'], info(4));
xermsg('SLATEC','DDES',['IN DDEABM, INFO(4) MUST BE ',['0 OR 1 INDICATING WHETHER OR NOT THE INTEGRATION ',['INTERVAL IS TO BE RESTRICTED BY A POINT TSTOP. YOU ',['HAVE CALLED THE CODE WITH INFO(4) = ',xern1]]]],14,1);
idid = -33;
end;
%
if( neq<1 )
xern1=sprintf(['%8i'], neq);
xermsg('SLATEC','DDES',['IN DDEABM, 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','DDES',['IN DDEABM, 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','DDES',['IN DDEABM, 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;
%
if( info(4)==1 )
if( (abs(1.0d0).*sign(tout-t))~=(abs(1.0d0).*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','DDES',['IN DDEABM, YOU HAVE ',['CALLED THE CODE WITH TOUT = ',[xern3,[' BUT ',['YOU HAVE ALSO TOLD THE CODE (INFO(4) = 1) NOT TO ',['INTEGRATE PAST THE POINT TSTOP = ',[xern4,' THESE INSTRUCTIONS CONFLICT.']]]]]]],14,1);
idid = -33;
end;
end;
%
% CHECK SOME CONTINUATION POSSIBILITIES
%
if( init~=0 )
if( t==tout )
xern3=sprintf([repmat('%15.6f',1,1)], t);
xermsg('SLATEC','DDES',['IN DDEABM, 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','DDES',['IN DDEABM, 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( delsgn.*(tout-t)<0.0d0 )
xern3=sprintf([repmat('%15.6f',1,1)], tout);
xermsg('SLATEC','DDES',['IN DDEABM, 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
%
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
% FOURU WHICH IS LIKELY TO BE REASONABLE FOR THIS METHOD AND MACHINE
%
for k = 1 : neq;
if( rtol(k)+atol(k)<=0.0d0 )
rtol(k) = fouru;
idid = -2;
end;
if( info(2)==0 )
break;
end;
end;
%
if( idid==(-2) )
% RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
% SMALL POSITIVE VALUE
info(1) = -1;
else;
%
% 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
%
gt=0;
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 )
idid = 2;
for l = 1 : neq;
ypout(l) = yp(l);
end; l = fix(neq+1);
told = t;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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;
end;
%
% -- SET INDEPENDENT AND DEPENDENT VARIABLES
% X AND YY(*) FOR STEPS
% -- SET SIGN OF INTEGRATION DIRECTION
% -- INITIALIZE THE STEP SIZE
%
if(gt==0)
init = 2;
x = t;
for l = 1 : neq;
yy(l) = y(l);
end; l = fix(neq+1);
delsgn = (abs(1.0d0).*sign(tout-t));
h = (abs(max(fouru.*abs(x),abs(tout-x))).*sign(tout-x));
end;
%
%.......................................................................
%
% ON EACH CALL SET INFORMATION WHICH DETERMINES THE ALLOWED INTERVAL
% OF INTEGRATION BEFORE RETURNING WITH AN ANSWER AT TOUT
%
del = tout - t;
absdel = abs(del);
%
%.......................................................................
%
% IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN
%
while( abs(x-t)<absdel );
%
% IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE,
% EXTRAPOLATE AND RETURN
%
if( info(4)==1 )
if( abs(tstop-x)<fouru.*abs(x) )
dt = tout - x;
for l = 1 : neq;
y(l) = yy(l) + dt.*yp(l);
end; l = fix(neq+1);
[tout,y,ypout,rpar,ipar]=df(tout,y,ypout,rpar,ipar);
idid = 3;
t = tout;
told = t;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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( ~(info(3)==0 || ~intout) )
%
% INTERMEDIATE-OUTPUT MODE
%
idid = 1;
for l = 1 : neq;
y(l) = yy(l);
ypout(l) = yp(l);
end; l = fix(neq+1);
t = x;
told = t;
intout = false;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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( ksteps>maxnum ) ;
%
% A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED
idid = -1;
ksteps = 0;
if( stiff )
%
% PROBLEM APPEARS TO BE STIFF
idid = -4;
stiff = false;
kle4 = 0;
end;
%
for l = 1 : neq;
y(l) = yy(l);
ypout(l) = yp(l);
end; l = fix(neq+1);
t = x;
told = t;
info(1) = -1;
intout = false;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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;
%
%.......................................................................
%
% LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP
%
ha = abs(h);
if( info(4)==1 )
ha = min(ha,abs(tstop-x));
end;
h = (abs(ha).*sign(h));
eps = 1.0d0;
ltol = 1;
for l = 1 : neq;
if( info(2)==1 )
ltol = fix(l);
end;
wt(l) = rtol(ltol).*abs(yy(l)) + atol(ltol);
if( wt(l)<=0.0d0 )
%
% RELATIVE ERROR CRITERION INAPPROPRIATE
idid = -3;
for ll = 1 : neq;
y(ll) = yy(ll);
ypout(ll) = yp(ll);
end; ll = fix(neq+1);
t = x;
told = t;
info(1) = -1;
intout = false;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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; l = fix(neq+1);
%
[df,neq,yy,x,h,eps,wt,start,hold,kord,kold,crash,phi,p,yp,psi,alpha,beta,sig,v,w,g,phase1,ns,nornd,ksteps,twou,fouru,xold,kprev,ivc,iv,kgi,gi,rpar,ipar]=dsteps(df,neq,yy,x,h,eps,wt,start,hold,kord,kold,crash,phi,p,yp,psi,alpha,beta,sig,v,w,g,phase1,ns,nornd,ksteps,twou,fouru,xold,kprev,ivc,iv,kgi,gi,rpar,ipar);
%
%.......................................................................
%
if( crash )
%
% TOLERANCES TOO SMALL
idid = -2;
rtol(1) = eps.*rtol(1);
atol(1) = eps.*atol(1);
if( info(2)~=0 )
for l = 2 : neq;
rtol(l) = eps.*rtol(l);
atol(l) = eps.*atol(l);
end; l = fix(neq+1);
end;
for l = 1 : neq;
y(l) = yy(l);
ypout(l) = yp(l);
end; l = fix(neq+1);
t = x;
told = t;
info(1) = -1;
intout = false;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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;
%
% (STIFFNESS TEST) COUNT NUMBER OF CONSECUTIVE STEPS TAKEN WITH THE
% ORDER OF THE METHOD BEING LESS OR EQUAL TO FOUR
%
kle4 = fix(kle4 + 1);
if( kold>4 )
kle4 = 0;
end;
if( kle4>=50 )
stiff = true;
end;
intout = true;
end;
end;
end;
[x,yy,tout,y,ypout,neq,kold,phi,ivc,iv,kgi,gi,alpha,g,w,xold,p]=dintp(x,yy,tout,y,ypout,neq,kold,phi,ivc,iv,kgi,gi,alpha,g,w,xold,p);
idid = 3;
if( x==tout )
idid = 2;
intout = false;
end;
t = tout;
told = t;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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','DDES',['IN DDEABM, 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);
else;
iquit = -33;
info(1) = -1;
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yy_shape=zeros(yy_shape);yy_shape(:)=yy(1:numel(yy_shape));yy=yy_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
phi_orig(1:prod(phi_shape))=phi;phi=phi_orig;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_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 ddes
%DECK DDNTL
|
|
