Line data Source code
1 0 : subroutine eiginh(x,y,rhs,ii,iy,ig,isn,nd1,nn)
2 : c
3 : c Initial value integrator.
4 : c =========================
5 : c Integrates a set of two first-order ordinary
6 : c linear homogeneous differential equations, by means of
7 : c constant-coefficient technique (Gabriel & Noels, A and A, vol. 53,
8 : c p. 149, 1976).
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))
13 : c j
14 : c
15 : c The coefficients are set up in the routine rhs, which must be
16 : c supplied by the user as external.
17 : c
18 : c The routine finds the solution
19 : c at every isn-th point in the input mesh. For ordinary use
20 : c isn is set to 1.
21 : c
22 : c nd1 specifies the initial value of the mesh index passed into the
23 : c routine rhs (see below).
24 : c
25 : c iy is the first dimension of y in the calling programme.
26 : c
27 : c The arguments ii and ig are included for consistency with
28 : c subroutine linint.
29 : c
30 : c On input x(1+isn*(n-1)), n = 1,nn, must contain the mesh in the
31 : c independent variable. The initial conditions for set k must be
32 : c set into y(i,1), i = 1,ii.
33 : c
34 : c The routine returns the solution for initial value set k in
35 : c y(i,1+isn*(n-1)), i=1,ii, n=1,nn.
36 : c
37 : c The right hand side subroutine rhs:
38 : c -----------------------------------
39 : c
40 : c This is called by linint and should be defined as
41 : c
42 : c subroutine rhs(x,aa,finh,iaa,nd)
43 : c dimension aa(iaa,1)
44 : c .
45 : c .
46 : c .
47 : c
48 : c A single call must set the coefficient matrix at a single point.
49 : c On input x (a scalar) gives the value of the independent variable
50 : c at the given point; iaa is the first dimension of aa;
51 : c nd is passed from linint and may be used to address a data
52 : c array in rhs.
53 : c
54 : c finh may later be used to set inhomogeneous term
55 : c
56 : c The routine should return
57 : c
58 : c aa(i,j) = a(i,j; x)
59 : c
60 : c where a(i,j; x) as defined above is the coefficient matrix in
61 : c the equations.
62 : c
63 : c The definition of nd is slightly convoluted. It is assumed that the
64 : c use might need certain variables to define the right hand side
65 : c of the equations. Assume, for example, that these are passed into
66 : c the routine in
67 : c
68 : c common/rhsdat/ data(10,200)
69 : c
70 : c When rhs is called at the meshpoint x(1+isn*(n-1)), nd is
71 : c set to nd1+isn*(n-1). The meshpoint should therefore correspond
72 : c to the data in data(i,nd). In most cases presumably x and the
73 : c data would be given on the same mesh, and nd1 would be 1.
74 : c
75 : c Note: this version is consistent with s/r lininh to integrate
76 : c inhomogeneous equations, in the call of rhs.
77 : c However, a non-zero inhomogeneous terms has not yet been
78 : c implemented.
79 : c
80 : c Original version: 9/7/95
81 : c
82 : implicit double precision(a-h,o-z)
83 : common/clscon/ icls
84 : c
85 : dimension x(nn),y(iy,nn),f(2,4),finh(2),cc(2,20),al(2),a(2,2)
86 : external rhs
87 0 : ifd=2
88 0 : icls=0
89 : c right hand sides at first point
90 0 : 5 n=1
91 0 : nc=1
92 0 : nd=nd1
93 0 : if1=1
94 0 : if2=3
95 0 : x2=x(1)
96 0 : call rhs(x2,f(1,if1),finh,ifd,nd)
97 : c right hand sides at next point
98 0 : 10 nd=nd+isn
99 0 : nc=nc+1
100 0 : n1=n
101 0 : n=n+isn
102 0 : x1=x2
103 0 : x2=x(n)
104 0 : dx=x2-x1
105 : c
106 0 : call rhs(x2,f(1,if2),finh,ifd,nd)
107 : c set mean coefficient matrix
108 0 : j1=if1-1
109 0 : j2=if2-1
110 0 : do 15 i=1,2
111 0 : do 15 j=1,2
112 0 : 15 a(i,j)=0.5d0*(f(i,j1+j)+f(i,j2+j))
113 : c
114 0 : do 17 i=1,2
115 0 : 17 y(i,n)=y(i,n1)
116 : c integrate constant coefficient equation
117 0 : call ieigst(a,y(1,n),dx,ig,al,cc,ie)
118 : c increment icls by contribution from this interval
119 0 : icls=icls+iclcon(x2,a,dx,al,cc,ie)
120 : c.. write(72,72090) n, x(n), y(1,n), y(2,n), icls
121 : c..72090 format(i5,f12.7,1p2e13.5,i5)
122 : c
123 0 : if(nc.eq.nn) return
124 0 : i=if1
125 0 : if1=if2
126 0 : if2=i
127 : c
128 0 : go to 10
129 : c
130 : end
131 0 : subroutine ieigst(a,yy,tt,ig,al,cc,ie)
132 : c integrates second order constant coefficient equation from 0
133 : c to tt
134 : implicit double precision(a-h,o-z)
135 : dimension a(2,2),yy(1),al(2),cc(2,1),x(2,2)
136 : c eigenvalues of a
137 0 : do 12 i=1,2
138 0 : al(i)=0
139 0 : do 12 j=1,2
140 0 : 12 x(i,j)=0
141 0 : call egenv2(a,x,al,ie)
142 : c set constants in analytical solution and solution at tt
143 0 : dd=x(1,1)*x(2,2)-x(1,2)*x(2,1)
144 : c
145 0 : i1=-1
146 0 : do 30 k=1,ig
147 0 : i1=i1+2
148 0 : i2=i1+1
149 0 : i0=i1-1
150 0 : y1=yy(i1)
151 0 : y2=yy(i2)
152 0 : ca=(x(2,2)*y1-x(2,1)*y2)/dd
153 0 : cb=(x(1,1)*y2-x(1,2)*y1)/dd
154 : c
155 0 : if(ie-1) 14,17,22
156 0 : 14 do 15 i=1,2
157 0 : cc(i,i1)=ca*x(1,i)+cb*x(2,i)
158 0 : 15 cc(i,i2)=cb*x(1,i)-ca*x(2,i)
159 0 : amt=al(2)*tt
160 0 : elt=exp(al(1)*tt)
161 0 : c=cos(amt)*elt
162 0 : s=sin(amt)*elt
163 0 : do 16 i=1,2
164 0 : 16 yy(i+i0)=cc(i,i1)*c+cc(i,i2)*s
165 0 : go to 30
166 : c
167 0 : 17 all=0.5d0*(a(1,1)-a(2,2))
168 0 : if(x(2,1).eq.0) go to 18
169 0 : cc(1,i2)=all
170 0 : cc(2,i2)=a(2,1)
171 0 : go to 19
172 0 : 18 cc(2,i2)=-all
173 0 : cc(1,i2)=a(1,2)
174 0 : 19 do 20 i=1,2
175 0 : cc(i,i1)=ca*x(1,i)+cb*x(2,i)
176 0 : 20 cc(i,i2)=cb*cc(i,i2)
177 0 : elt=exp(al(1)*tt)
178 0 : do 21 i=1,2
179 0 : 21 yy(i+i0)=(cc(i,i1)+tt*cc(i,i2))*elt
180 0 : go to 30
181 : c
182 0 : 22 do 23 i=1,2
183 0 : cc(i,i1)=ca*x(1,i)
184 0 : 23 cc(i,i2)=cb*x(2,i)
185 0 : el1=exp(al(1)*tt)
186 0 : el2=exp(al(2)*tt)
187 0 : do 24 i=1,2
188 0 : 24 yy(i+i0)=cc(i,i1)*el1+cc(i,i2)*el2
189 0 : 30 continue
190 0 : return
191 : end
192 0 : subroutine egenv2(a,x,al,i)
193 : c finds eigenvectors and eigenvalues of real 2*2 matrix a.
194 : c output:
195 : c ======
196 : c for i.ge.1: i (1 or 2) real eigenvalues corresponding to different
197 : c eigenvectors (the eigenvalues may be the same). the j-th eigenvalue
198 : c is in al(j) and the j-th eigenvector in x(j,1),x(j,2), for j=1,i.
199 : c
200 : c for i = -1: two complex conjugated eigenvalues. the eigenvalues are
201 : c al(1) +- i*al(2), and the eigenvectors have the k-th component
202 : c x(1,k) +- i*x(2,k), k=1,2.
203 : implicit double precision(a-h,o-z)
204 : dimension a(2,2),x(2,2),al(2)
205 0 : eps=1.d-10
206 0 : eps2=eps*eps
207 0 : alp=(a(1,1)+a(2,2))/2
208 0 : y=a(2,2)-a(1,1)
209 0 : x1=a(1,2)*a(2,1)
210 0 : if(abs(x1).lt.eps2)goto 2
211 0 : dis=y*y+4*x1
212 0 : if(abs(dis).gt.eps2) goto 6
213 : c dis = 0, one real eigenvalue
214 0 : i=1
215 0 : al(1)=alp
216 0 : x(1,1)=2*a(1,2)/y
217 0 : x(1,2)=1.d0
218 : c set x(2,.)
219 0 : if(abs(a(1,2)).gt.abs(a(2,1))) go to 9
220 0 : x(2,1)=1.d0
221 0 : x(2,2)=0.d0
222 0 : return
223 0 : 9 x(2,1)=0.d0
224 0 : x(2,2)=1.d0
225 0 : return
226 : c dis.ne.0
227 0 : 6 if(dis.lt.0.0d0) goto 1
228 : c two real eigenvalues
229 0 : dis=sqrt(dis)/2
230 0 : d=a(1,1)*a(2,2)-x1
231 0 : if(alp) 11,11,12
232 0 : 11 al(2)=alp-dis
233 0 : al(1)=d/al(2)
234 0 : go to 14
235 0 : 12 al(1)=alp+dis
236 0 : al(2)=d/al(1)
237 0 : 14 do 16 j=1,2
238 0 : x(j,2)=1
239 0 : 16 x(j,1)=a(1,2)/(al(j)-a(1,1))
240 0 : i=2
241 0 : return
242 : c two complex eigenvalues
243 0 : 1 i=-1
244 0 : al(1)=alp
245 0 : al(2)=sqrt(-dis)/2
246 0 : xx=alp-a(1,1)
247 0 : xx1=xx*xx-dis/4
248 0 : x(1,1)=a(1,2)*xx/xx1
249 0 : x(1,2)=1.0d0
250 0 : x(2,1)=-a(1,2)*al(2)/xx1
251 0 : x(2,2)=0.0d0
252 0 : return
253 : c zero off-diagonal element
254 0 : 2 if(abs(y).lt.eps) goto 3
255 : c two different real eigenvalues
256 0 : al(1)=a(1,1)
257 0 : al(2)=a(2,2)
258 0 : x(1,1)=1.d0
259 0 : x(2,2)=1.d0
260 0 : x(1,2)=-a(2,1)/y
261 0 : x(2,1)=a(1,2)/y
262 0 : i=2
263 0 : return
264 : c diagonal elements equal. one or two eigenvectors
265 0 : 3 i=0
266 0 : if(abs(a(2,1)).gt.eps) goto 4
267 0 : i=1
268 0 : al(1)=a(1,1)
269 0 : x(1,1)=1.d0
270 0 : x(1,2)=0.d0
271 0 : if(abs(a(1,2)).lt.eps) goto 5
272 : c a(1,2).ne.0, set x(2,.)
273 0 : x(2,1)=0.d0
274 0 : x(2,2)=1.d0
275 0 : return
276 : c a(2,1).ne.0., set x(2,.)
277 0 : 4 x(2,1)=1.d0
278 0 : x(2,2)=0.d0
279 : c
280 0 : 5 i=i+1
281 0 : al(i)=a(2,2)
282 0 : x(i,1)=0.d0
283 0 : x(i,2)=1.d0
284 0 : return
285 : end
286 0 : integer function iclcon(x,a,tt,al,cc,ie)
287 : c finds contribution to classification index from interval (0,tt),
288 : c from solution of constant coefficient equation.
289 : implicit double precision(a-h,o-z)
290 : dimension a(2,2),al(2),cc(2,*)
291 : data pi/3.14159265358979d0/
292 : c find number of zeros of y1
293 0 : c1=cc(1,1)
294 0 : c2=cc(1,2)
295 0 : icl=0
296 0 : if(c1.eq.0.and.c2.eq.0) go to 50
297 : c
298 0 : 10 if(ie-1) 20,30,40
299 : c oscillatory solution
300 0 : 20 if(al(2).eq.0) go to 50
301 0 : delta=atan2(c2,c1)
302 0 : zt0=-delta/pi-0.5d0
303 0 : ztt=al(2)*tt/pi+zt0
304 : c test for direction of integration
305 0 : if(tt) 22,50,24
306 0 : 22 zt1=ztt
307 0 : zt2=zt0
308 0 : go to 26
309 : c
310 0 : 24 zt1=zt0
311 0 : zt2=ztt
312 : c now zt2 corresponds to the larger x, zt1 to the smaller
313 0 : 26 izt1=intgpt(-zt1)
314 0 : izt2=intgpt(zt2)
315 0 : icl=1+izt1+izt2
316 : c test for zero at smallest x
317 0 : if(izt1.eq.-zt1) icl=icl-1
318 0 : go to 50
319 : c degenerate and exponential cases
320 : c
321 : c find zero
322 0 : 30 if(c2.eq.0) go to 50
323 0 : tz=-c1/c2
324 0 : go to 45
325 : c
326 0 : 40 if(c1*c2.ge.0.or.al(1).eq.al(2)) go to 50
327 0 : tz=log(-c2/c1)/(al(1)-al(2))
328 : c test for direction of integration
329 0 : 45 if(tt) 46,50,48
330 : c test for position of zero
331 0 : 46 if(tt.lt.tz.and.tz.le.0) icl=1
332 0 : go to 50
333 0 : 48 if(0.lt.tz.and.tz.le.tt) icl=1
334 : c find sign of a(1,2)
335 0 : 50 isg=0
336 0 : if(a(1,2).ne.0) isg=sign(1.d0,a(1,2))
337 : c contribution to index
338 0 : iclcon=-isg*icl
339 : c.. write(73,73090) x,zt1,zt2,izt1,izt2,icl
340 : c..73090 format(f12.7,1p2e13.5,3i10)
341 0 : return
342 : end
|
