Line data Source code
1 0 : subroutine lininh(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 linear differential equation solver.
80 : c differs from original version
81 : c of lininh by incorporating gaussian elimination without pivoting
82 : c in subroutine.
83 : c
84 : c original version: 27/8/1991
85 : c
86 : implicit double precision(a-h,o-z)
87 : dimension x(nn),y(iy,nn),fd(20,20),finh(20),w(20,20),y1(20)
88 : common/modfac/ r21,ntld
89 : c
90 : c common defining standard input and output
91 : c
92 : common/cstdio/ istdin, istdou, istdpr, istder
93 : external rhs
94 : c
95 : c.. write(6,*) 'Enter lininh with ii, iy, ig =',ii, ig, iy
96 : c.. write(6,*) 'y(1-2,1) =',(y(i,1),i=1,2)
97 : c
98 : c explicitly disable check for linear dependence
99 : c
100 0 : ifd=20
101 : c
102 0 : iw=20
103 : c right hand sides at first point
104 0 : 5 n=1
105 0 : nc=1
106 0 : nd=nd1
107 0 : ntst=nn-2
108 0 : iig=ii*ig
109 0 : x2=x(1)
110 0 : call rhs(x2,fd,finh,ifd,nd)
111 0 : go to 20
112 : c right hand sides at next point
113 0 : 10 nd=nd+isn
114 0 : call rhs(x2,fd,finh,ifd,nd)
115 : c set coefficient matrix for linear equations
116 0 : do 14 i=1,ii
117 0 : do 13 j=1,ii
118 0 : 13 w(i,j)=dx*fd(i,j)
119 0 : 14 w(i,i)=1+w(i,i)
120 : c
121 : c include inhomogeneous contribution from next point
122 : c
123 0 : do 15 i=1,iig
124 0 : 15 y(i,n)=y(i,n)-dx*finh(i)
125 : c
126 : c.. write(6,15091) n,x(n),((w(i,j),j=1,2),y(i,n),i=1,2)
127 : c..15091 format(' Equations at n, x =',i5,f10.5/(1p3e13.5))
128 : c
129 : c solve linear equations by gaussian elimination without pivoting
130 : c
131 0 : ii1=ii-1
132 0 : if(ii1.gt.0) then
133 : c
134 : c triangularize matrix
135 : c
136 0 : do 17 i=1,ii1
137 0 : r=w(i,i)
138 : c
139 : c test for non-zero diagonal element
140 : c
141 0 : if(r.eq.0) then
142 0 : if(istdpr.gt.0) write(istdpr,110) i,n
143 0 : return
144 : end if
145 : c
146 0 : r=-1.d0/r
147 : c
148 0 : i1=i+1
149 : c
150 0 : do 17 j=i1,ii
151 0 : rj=r*w(j,i)
152 0 : do 16 k=i1,ii
153 0 : 16 w(j,k)=w(j,k)+rj*w(i,k)
154 : c
155 0 : js=j
156 0 : is=i
157 : c
158 0 : do 17 l=1,ig
159 0 : y(js,n)=y(js,n)+rj*y(is,n)
160 0 : js=js+ii
161 0 : 17 is=is+ii
162 : c
163 : end if
164 : c
165 : c now matrix is on triangular form
166 : c
167 : c start solution
168 : c
169 0 : i=ii+1
170 0 : do 18 icnt=1,ii
171 0 : i=i-1
172 : c
173 0 : r=w(i,i)
174 : c
175 : c test for non-zero diagonal element
176 : c
177 0 : if(r.eq.0) then
178 0 : if(istdpr.gt.0) write(istdpr,110) i,n
179 0 : return
180 : end if
181 : c
182 0 : r=1.d0/r
183 0 : is=i
184 0 : i1=i+1
185 : c
186 0 : do 18 l=1,ig
187 0 : sum=y(is,n)
188 : c
189 0 : if(i.lt.ii) then
190 : c
191 0 : j1=is
192 0 : do 17100 j=i1,ii
193 0 : j1=j1+1
194 0 : 17100 sum=sum-w(i,j)*y(j1,n)
195 : c
196 : end if
197 : c
198 0 : y(is,n)=r*sum
199 : c
200 0 : 18 is=is+ii
201 0 : if(nc.eq.nn) return
202 : c
203 : c set up right hand side of linear equation at next point
204 : c
205 0 : 20 n1=n
206 0 : n=n+isn
207 0 : nc=nc+1
208 0 : x1=x2
209 0 : x2=x(n)
210 0 : dx=0.5d0*(x1-x2)
211 0 : kg=-ii
212 0 : do 25 l=1,ig
213 0 : kg=kg+ii
214 0 : do 25 i=1,ii
215 0 : k=kg+i
216 0 : sum=0
217 0 : do 22 j=1,ii
218 0 : sum=sum+fd(i,j)*y(kg+j,n1)
219 : c.. write(6,*) 'i,j,kg+j,n1,fd(i,j),y(kg+j,n1),sum:'
220 : c.. write(6,*) i,j,kg+j,n1,fd(i,j),y(kg+j,n1),sum
221 0 : 22 continue
222 0 : y(k,n)=y(k,n1)-dx*(sum+finh(k))
223 : c.. write(6,*) 'k, y(k,n1), dx, sum, finh(k):'
224 : c.. write(6,*) k, y(k,n1), dx, sum, finh(k)
225 0 : 25 continue
226 0 : if(ntld.ne.1) go to 10
227 0 : if(nc.ne.ntst.or.ig.ne.2) go to 10
228 : c test for linear dependence (only implemented for ig = 2)
229 0 : do 30 i=1,iig
230 0 : 30 y1(i)=y(i,n1)
231 0 : aym=0
232 0 : im=0
233 0 : aym2=0
234 0 : do 32 i=1,ii
235 0 : ay=abs(y1(i))
236 0 : aym2=max(aym2,abs(y1(i+ii)))
237 0 : if(ay.le.aym) go to 32
238 0 : aym=ay
239 0 : im=i
240 0 : 32 continue
241 : c
242 0 : r21=y1(ii+im)/y1(im)
243 0 : s1=0
244 0 : s2=0
245 0 : j=ii
246 0 : do 34 i=1,ii
247 0 : j=j+1
248 0 : yi=y1(i)
249 0 : if(yi.ne.0) go to 33
250 0 : if(abs(y1(j)).gt.1.d-5*aym2) go to 10
251 0 : go to 34
252 0 : 33 ri=y1(j)/y1(i)
253 0 : s1=s1+abs(r21-ri)
254 0 : s2=s2+abs(ri)
255 0 : 34 continue
256 : c test
257 0 : if(s1/s2.gt.1.d-7) go to 10
258 : c nearly linear dependence. modify initial values and start again
259 0 : j=0
260 0 : aym=0
261 0 : do 36 i=1,ii
262 0 : y1(i)=y(i+ii,1)-r21*y(i,1)
263 0 : 36 aym=max(aym,abs(y1(i)))
264 : c normalize initial y's
265 0 : aym=1.d0/aym
266 0 : do 38 i=1,ii
267 0 : 38 y(ii+i,1)=aym*y1(i)
268 : c
269 0 : if(istdpr.gt.0) write(istdpr,100) (y(i,1),i=1,iig)
270 0 : go to 5
271 : c
272 : 100 format(//' initial solution has been modified. new y(n=1):'/
273 : * (1p10e13.5))
274 : c
275 : 110 format(//' ***** in s/r lininh zero diagonal element at i =',
276 : * i5,' n =',i5)
277 : end
|
