Line data Source code
1 1784 : subroutine lininh4(x,y,rhs,ii,iy,ig,isn,nd1,nn)
2 : c
3 : c
4 : c Initial value integrator.
5 : c =========================
6 : c Integrates a set of first-order ordinary
7 : c linear inhomogeneous differential equations, by means of second-order
8 : c centred difference approximation, from given initial values.
9 : c
10 : c The equations are on the form
11 : c
12 : c d y(i; x)/dx = sum ( a(i,j; x) * y(j; x)) + f(i; x)
13 : c j
14 : c
15 : c The order of the equations is ii.
16 : c
17 : c The coefficients and the inhomogeneous terms
18 : c are set up in the routine rhs, which must be
19 : c supplied by the user as external.
20 : c
21 : c The routine has the option for determining in a single call the
22 : c solutions to a given set of equations for several initial conditions;
23 : c ig specifies the number of such sets.
24 : c
25 : c For use with Richardson extrapolation the routine finds the solution
26 : c at every isn-th point in the input mesh. For ordinary use
27 : c isn is set to 1.
28 : c
29 : c nd1 specifies the initial value of the mesh index passed into the
30 : c routine rhs (see below).
31 : c
32 : c iy is the first dimension of y in the calling programme.
33 : c
34 : c On input x(1+isn*(n-1)), n = 1,nn, must contain the mesh in the
35 : c independent variable. The initial conditions for set k must be
36 : c set into y(i + ii*(k-1),1), i = 1,ii, for k = 1,ig.
37 : c
38 : c The routine returns the solution for initial value set k in
39 : c y(i+ii*(k-1),1+isn*(n-1)), i=1,ii, n=1,nn, k=1,ig.
40 : c
41 : c The right hand side subroutine rhs:
42 : c -----------------------------------
43 : c
44 : c This is called by lininh and should be defined as
45 : c
46 : c subroutine rhs(x,aa,finh,iaa,nd)
47 : c dimension aa(iaa,1),finh(1)
48 : c .
49 : c .
50 : c .
51 : c
52 : c A single call must set the coefficient matrix and inhomogeneous
53 : c term at a single point.
54 : c On input x (a scalar) gives the value of the independent variable
55 : c at the given point; iaa is the first dimension of aa;
56 : c nd is passed from lininh and may be used to address a data
57 : c array in rhs.
58 : c
59 : c The routine should return
60 : c
61 : c aa(i,j) = a(i,j; x)
62 : c finh(i + ii*(k-1)) = f(i; x) for the k-th set, k = 1,ig.
63 : c
64 : c where a(i,j; x) and f(i; x) as defined above are the
65 : c coefficient matrix and inhomogeneous term in the equations.
66 : c
67 : c The definition of nd is slightly convoluted. It is assumed that the
68 : c use might need certain variables to define the right hand side
69 : c of the equations. Assume, for example, that these are passed into
70 : c the routine in
71 : c
72 : c common/rhsdat/ data(10,200)
73 : c
74 : c When rhs is called at the meshpoint x(1+isn*(n-1)), nd is
75 : c set to nd1+isn*(n-1). The meshpoint should therefore correspond
76 : c to the data in data(i,nd). In most cases presumably x and the
77 : c data would be given on the same mesh, and nd1 would be 1.
78 : c
79 : c 4th order version, using algorithm of Cash & Moore (1980; BIT 20:1, 44 - 52)
80 : c
81 : c original version: 11/3/05
82 : c
83 : implicit double precision(a-h,o-z)
84 : dimension x(nn),y(iy,nn),fd(20,20),fdp(20,20),fdh(20,20),finh(20),
85 : * fp(20),w(20,20),y1(20)
86 : common/modfac/ r21,ntld
87 : c
88 : c common defining standard input and output
89 : c
90 : common/cstdio/ istdin, istdou, istdpr, istder
91 : external rhs
92 : c
93 : c.. write(6,*) 'Enter lininh with ii, iy, ig =',ii, ig, iy
94 : c.. write(6,*) 'y(1-2,1) =',(y(i,1),i=1,2)
95 : c
96 : c explicitly disable check for linear dependence
97 : c
98 1784 : ifd=20
99 : c
100 1784 : iw=20
101 : c right hand sides at first point
102 1784 : 5 n=1
103 1784 : nc=1
104 1784 : nd=nd1
105 1784 : ntst=nn-2
106 1784 : iig=ii*ig
107 1784 : x2=x(1)
108 1784 : call rhs(x2,fd,finh,ifd,nd)
109 : c
110 : c right hand sides at next point
111 : c
112 4280720 : 10 do i=1,ii
113 16665488 : do j=1,ii
114 15715520 : fdp(i,j)=fd(i,j)
115 : end do
116 : end do
117 : c
118 : c set up right hand side of linear equation at next point
119 : c
120 949968 : n1=n
121 949968 : n=n+isn
122 949968 : nc=nc+1
123 949968 : x1=x2
124 949968 : x2=x(n)
125 : c
126 949968 : nd=nd+isn
127 : c
128 949968 : call rhs(x2,fd,finh,ifd,nd)
129 : c
130 4280720 : do i=1,ii
131 16665488 : do j=1,ii
132 15715520 : fdh(i,j)=0.5d0*(fd(i,j)+fdp(i,j))
133 : end do
134 : end do
135 949968 : dx=x2-x1
136 949968 : kg=-ii
137 2615344 : do 25 l=1,ig
138 1665376 : kg=kg+ii
139 7857760 : do 23 i=1,ii
140 6192384 : k=kg+i
141 6192384 : sum=0
142 30023680 : do 22 j=1,ii
143 23831296 : sum=sum+fdp(i,j)*y(kg+j,n1)
144 6192384 : 22 continue
145 6192384 : fp(i)=sum
146 1665376 : 23 continue
147 : c
148 7857760 : do i=1,ii
149 6192384 : sum1=0
150 6192384 : sum2=0
151 30023680 : do j=1,ii
152 23831296 : sum1=sum1+fdh(i,j)*y(kg+j,n1)
153 30023680 : sum2=sum2+fdh(i,j)*fp(j)
154 : end do
155 : y(kg+i,n)=y(kg+i,n1)+dx*fp(i)/6.d0+dx*sum1/3.d0
156 7857760 : * +dx*dx*sum2/12.d0
157 : end do
158 : c
159 949968 : 25 continue
160 : c set coefficient matrix for linear equations
161 4280720 : do i=1,ii
162 15715520 : do j=1,ii
163 12384768 : sum=0
164 60047360 : do k=1,ii
165 60047360 : sum=sum+fdh(i,k)*fd(k,j)
166 : end do
167 15715520 : w(i,j)=-(fd(i,j)+2.d0*fdh(i,j))*dx/6.d0+sum*dx*dx/12.d0
168 : end do
169 4280720 : w(i,i)=1.d0+w(i,i)
170 : end do
171 : c
172 : c include inhomogeneous contribution from next point (needs to be added)
173 : c
174 : c.. do 15 i=1,iig
175 : c.. 15 y(i,n)=y(i,n)-dx*finh(i)
176 : c
177 : c.. write(6,15091) n,x(n),((w(i,j),j=1,2),y(i,n),i=1,2)
178 : c..15091 format(' Equations at n, x =',i5,f10.5/(1p3e13.5))
179 : c
180 : c solve linear equations by gaussian elimination without pivoting
181 : c
182 949968 : ii1=ii-1
183 949968 : if(ii1.gt.0) then
184 : c
185 : c triangularize matrix
186 : c
187 3330752 : do 27 i=1,ii1
188 2380784 : r=w(i,i)
189 : c
190 : c test for non-zero diagonal element
191 : c
192 2380784 : if(r.eq.0) then
193 0 : if(istdpr.gt.0) write(istdpr,110) i,n
194 0 : return
195 : end if
196 : c
197 2380784 : r=-1.d0/r
198 : c
199 2380784 : i1=i+1
200 : c
201 7857760 : do 27 j=i1,ii
202 4527008 : rj=r*w(j,i)
203 14777280 : do 26 k=i1,ii
204 14777280 : 26 w(j,k)=w(j,k)+rj*w(i,k)
205 : c
206 4527008 : js=j
207 4527008 : is=i
208 : c
209 15727248 : do 27 l=1,ig
210 8819456 : y(js,n)=y(js,n)+rj*y(is,n)
211 8819456 : js=js+ii
212 13346464 : 27 is=is+ii
213 : c
214 : end if
215 : c
216 : c now matrix is on triangular form
217 : c
218 : c start solution
219 : c
220 949968 : i=ii+1
221 4280720 : do 28 icnt=1,ii
222 3330752 : i=i-1
223 : c
224 3330752 : r=w(i,i)
225 : c
226 : c test for non-zero diagonal element
227 : c
228 3330752 : if(r.eq.0) then
229 0 : if(istdpr.gt.0) write(istdpr,110) i,n
230 0 : return
231 : end if
232 : c
233 3330752 : r=1.d0/r
234 3330752 : is=i
235 3330752 : i1=i+1
236 : c
237 10473104 : do 28 l=1,ig
238 6192384 : sum=y(is,n)
239 : c
240 6192384 : if(i.lt.ii) then
241 : c
242 4527008 : j1=is
243 13346464 : do j=i1,ii
244 8819456 : j1=j1+1
245 13346464 : sum=sum-w(i,j)*y(j1,n)
246 : end do
247 : c
248 : end if
249 : c
250 6192384 : y(is,n)=r*sum
251 : c
252 9523136 : 28 is=is+ii
253 949968 : if(nc.eq.nn) return
254 : c
255 948184 : go to 10
256 : c
257 : 110 format(//' ***** in s/r lininh zero diagonal element at i =',
258 : * i5,' n =',i5)
259 : end
|
