Bug #10574 fixed - Runge-Kutta-Method failed for vector [x, 1] with x > 9.
[scilab.git] / scilab / modules / differential_equations / src / fortran / rkqc.f
1 c     ====================================
2       subroutine rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs,rwork)
3 C      implicit undefined (a-z)
4 c ................................................................
5 c .   Subroutine rkqc is a fifth-order Runge-Kutta step with 
6 c .   monitoring of local truncation error to ensure accuracy
7 c .   and adjust stepsize.   
8 c .
9 c .   Inputs:
10 c .   y(1:n) = the dependent variable vector
11 c .   dydx(1:n) = deriv of y wrt x at starting value of indep var x
12 c .   htry = attempted step size
13 c .   eps = required accuracy
14 c .   yscal(1:n) = scaling vector for the error
15 c .   
16 c .   Outputs:
17 c .   y = new value of y
18 c .   x = new value of x
19 c .   hdid = actual stepsize accomplished
20 c .   hnext = estimated next stepsize   
21 c .
22 c .   derivs is the user-supplied subroutine that computes the 
23 c .    right-hand side derivatives
24 c ................................................................
25 c
26 c     The original version has been modified to replace statically
27 c     allocated arrays dysav, ytemp and ysav by rwork arguments parts
28 c     array + blas use. Serge Steer INRIA- feb 2012
29
30       integer n,i
31       double precision fcor,one,safety,errcon
32       parameter (fcor=.0666666667,one=1.0,safety=0.9,errcon=6.e-4)
33       double precision x,htry,eps,hdid,hnext,pgrow,pshrnk,xsav,hh
34       double precision errmax,h
35       double precision y(n),dydx(n),yscal(n)
36       double precision rwork(*)
37
38       external derivs
39       integer iero
40       common/ierode/iero
41
42       lysav=1
43       ldysav=lysav+n
44       lytemp=ldysav+n
45       lwork=lytemp+n
46
47       iero=0
48       pgrow=-0.20d0
49       pshrnk=-0.25d0
50       xsav=x
51
52       call dcopy(n,y,1,rwork(lysav),1)
53       call dcopy(n,dydx,1,rwork(ldysav),1)
54       h=htry
55 1     hh=0.5*h
56       call rk4(rwork(lysav),rwork(ldysav),n,xsav,hh,rwork(lytemp),
57      $     derivs,rwork(lwork))
58       x=xsav+hh
59       call derivs(n,x,rwork(lytemp),dydx)
60       if (iero.gt.0) return 
61       call rk4(rwork(lytemp),dydx,n,x,hh,y,derivs,rwork(lwork))
62       x=xsav+h
63       if(x.eq.xsav) then 
64          iero=1
65          return
66       endif
67       call rk4(rwork(lysav),rwork(ldysav),n,xsav,h,rwork(lytemp),
68      $     derivs,rwork(lwork))
69       errmax=0.0d0
70        do 12 i=0,n-1
71         rwork(lytemp+i)=y(i+1)-rwork(lytemp+i)
72         errmax=max(errmax,abs(rwork(lytemp+i)/(yscal(i+1)*eps)))
73 12    continue
74        
75       if(errmax.gt.one) then
76         h=safety*h*(errmax**pshrnk)
77         goto 1
78       else
79         hdid=h
80         if(errmax.gt.errcon) then
81           hnext=safety*h*(errmax**pgrow)
82         else
83           hnext=4.0d0*h
84         endif
85       endif
86       call daxpy(n,fcor,rwork(lytemp),1,y,1)
87
88       return
89       end
90 c     ====================================