← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/integ/rksimp.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/integ/rksimp.f
 
1
      subroutine rksimp(fydot2,neqn,y,t,tout,itol,rerr,aerr,
 
2
     1    itask,iflag,iopt,work,lrw,iwork,liw,bjac,mf)
 
3
c
 
4
c     fehlberg fourth-fifth order runge-kutta method
 
5
c
 
6
c
 
7
      integer neqn,iflag,iwork(*)
 
8
      double precision y(*),t,tout,rerr,aerr,work(*),h
 
9
c
 
10
      double precision ae,scale,eeoet,et,esttol,ee
 
11
      common/ierode/iero
 
12
      external fydot2
 
13
c
 
14
      integer k1,k2,k3,k4,k5,k6
 
15
      iero=0
 
16
c
 
17
c     compute indices for the splitting of the work array
 
18
c
 
19
      scale=2.0d0/rerr
 
20
      ae=scale*aerr
 
21
      k1=1+neqn
 
22
      k2=k1+neqn
 
23
      k3=k2+neqn
 
24
      k4=k3+neqn
 
25
      k5=k4+neqn
 
26
      k6=k5+neqn
 
27
c
 
28
c     this interfacing routine merely relieves the user of a long
 
29
c     calling list via the splitting apart of two working storage
 
30
c     arrays. if this is not compatible with the users compiler,
 
31
c     he must use rkfs directly.
 
32
c
 
33
      h=tout-t
 
34
      do 33 k=1,neqn
 
35
      work(k6+k-1)=y(k)
 
36
 33   continue
 
37
      call fehl2(fydot2,neqn,y,t,h,work(1),
 
38
     1          work(k1),work(k2),work(k3),work(k4),work(k5),
 
39
     2          work(k6))
 
40
c
 
41
      eeoet=0.0d0
 
42
      do 250 k=1,neqn
 
43
        et=dabs(work(k6+k-1))+dabs(work(k1-1+k))+ae
 
44
        if (et .gt. 0.0d0) go to 240
 
45
c
 
46
c       inappropriate error tolerance
 
47
        iflag=4
 
48
        return
 
49
c
 
50
 240    continue
 
51
        ee=dabs((-2090.0d0*work(k)+(21970.0d0*
 
52
     1      work(k3-1+k)-15048.0d0*work(k4-1+k)))+
 
53
     2    (22528.0d0*work(k2-1+k)-27360.0d0*work(k5-1+k)))
 
54
  250   eeoet=dmax1(eeoet,ee/et)
 
55
c
 
56
      esttol=dabs(h)*eeoet*scale/752400.0d0
 
57
c
 
58
      if (esttol .le. 1.0d0) then
 
59
c        OK
 
60
      iflag=2
 
61
      t=tout
 
62
c      write(6,*) esttol,scale,eeoet,ee,et
 
63
      return
 
64
      endif
 
65
      iflag=3
 
66
c
 
67
      return
 
68
      end
 
69
 
 
70
      subroutine fehl2(fydot2,neqn,y,t,h,yp,f1,f2,f3,f4,f5,s)
 
71
c
 
72
c     fehlberg fourth-fifth order runge-kutta method
 
73
c
 
74
c
 
75
      integer  neqn
 
76
      double precision  y(*),t,h,yp(neqn),f1(neqn),f2(neqn),
 
77
     1  f3(neqn),f4(neqn),f5(neqn),s(neqn)
 
78
c
 
79
      double precision  ch
 
80
      integer  k
 
81
      external fydot2
 
82
      common/ierode/iero
 
83
c
 
84
c      write(6,*) 'inputfelh2:',y(1),y(2)
 
85
      call fydot2(neqn,t,y,yp)
 
86
      if(iero.gt.0) return
 
87
      ch=h/4.0d0
 
88
      do 221 k=1,neqn
 
89
  221   y(k)=y(k)+ch*yp(k)
 
90
      call fydot2(neqn,t+ch,y,f1)
 
91
      if(iero.gt.0) return
 
92
c
 
93
      ch=3.0d0*h/32.0d0
 
94
      do 222 k=1,neqn
 
95
  222   y(k)=s(k)+ch*(yp(k)+3.0d0*f1(k))
 
96
      call fydot2(neqn,t+3.0d0*h/8.0d0,y,f2)
 
97
      if(iero.gt.0) return
 
98
c
 
99
      ch=h/2197.0d0
 
100
      do 223 k=1,neqn
 
101
  223   y(k)=s(k)+ch*(1932.0d0*yp(k)+(7296.0d0*f2(k)-7200.0d0*f1(k)))
 
102
      call fydot2(neqn,t+12.0d0*h/13.0d0,y,f3)
 
103
      if(iero.gt.0) return
 
104
c
 
105
      ch=h/4104.0d0
 
106
      do 224 k=1,neqn
 
107
  224   y(k)=s(k)+ch*((8341.0d0*yp(k)-845.0d0*f3(k))+
 
108
     1                            (29440.0d0*f2(k)-32832.0d0*f1(k)))
 
109
      call fydot2(neqn,t+h,y,f4)
 
110
      if(iero.gt.0) return
 
111
c
 
112
      ch=h/20520.0d0
 
113
      do 225 k=1,neqn
 
114
  225   y(k)=s(k)+ch*((-6080.0d0*yp(k)+(9295.0d0*f3(k)-
 
115
     1         5643.0d0*f4(k)))+(41040.0d0*f1(k)-28352.0d0*f2(k)))
 
116
      call fydot2(neqn,t+h/2.0d0,y,f5)
 
117
      if(iero.gt.0) return
 
118
c
 
119
c     compute approximate solution at t+h
 
120
c
 
121
      ch=h/7618050.0d0
 
122
      do 230 k=1,neqn
 
123
       y(k)=s(k)+ch*((902880.0d0*yp(k)+(3855735.0d0*f3(k)-
 
124
     1        1371249.0d0*f4(k)))+(3953664.0d0*f2(k)+
 
125
     2        277020.0d0*f5(k)))
 
126
c
 
127
 230  continue
 
128
c      write(6,*) 'endfelh2:',y(1),y(2)
 
129
      return
 
130
      end