Line data Source code
1 : ! ----------------------------------------------------------------------
2 : !
3 : ! This file is part of the Test Set for IVP solvers
4 : ! http://www.dm.uniba.it/~testset/
5 : !
6 : ! Beam (ODE case)
7 : ! ODE of dimension 80
8 : !
9 : ! DISCLAIMER: see
10 : ! http://www.dm.uniba.it/~testset/disclaimer.php
11 : !
12 : ! The most recent version of this source file can be found at
13 : ! http://www.dm.uniba.it/~testset/src/problems/beam.f
14 : !
15 : ! This is revision
16 : ! $Id: beam.F,v 1.2 2006/10/02 10:29:13 testset Exp $
17 : !
18 : ! ----------------------------------------------------------------------
19 : module bari_beam
20 :
21 : use const_def, only: dp
22 :
23 : implicit none
24 :
25 : contains
26 :
27 0 : subroutine beam_init(neqn, y, yprime, consis)
28 : integer :: neqn
29 : real(dp) :: y(neqn), yprime(neqn)
30 : logical :: consis
31 :
32 : integer :: i
33 :
34 0 : do i = 1, neqn
35 0 : y(i) = 0d0
36 : end do
37 :
38 0 : end subroutine beam_init
39 : ! ----------------------------------------------------------------------
40 0 : subroutine beam_feval(nvar, t, th, df, ierr, rpar, ipar)
41 : use math_lib
42 : implicit real(dp) (A - H, O - Z)
43 : integer :: ierr, nvar, i, ipar(*)
44 : integer, parameter :: N = 40, NN = 2*N, NSQ = N*N, NQUATR = NSQ*NSQ
45 : real(dp) :: rpar(*), an, deltas
46 0 : dimension DF(NN), TH(150), U(150), V(150), W(150)
47 0 : dimension ALPHA(150), BETA(150), STH(150), CTH(150)
48 : ! --- SET DEFAULT VALUES
49 0 : if (nvar /= nn) stop 'bad nvar for beam_feval'
50 : AN = N
51 : DELTAS = 1.0D+0/AN
52 : ! ----- CALCUL DES TH(I) ET DES SIN ET COS -------------
53 0 : do I = 2, N
54 0 : THDIFF = TH(I) - TH(I - 1)
55 0 : STH(I) = sin(THDIFF)
56 0 : CTH(I) = cos(THDIFF)
57 : end do
58 : ! -------- CALCUL DU COTE DROIT DU SYSTEME LINEAIRE -----
59 0 : if (T > 3.14159265358979324D0) THEN
60 : ! --------- CASE T GREATER PI ------------
61 : ! ---------- I=1 ------------
62 0 : TERM1 = (-3.D0*TH(1) + TH(2))*NQUATR
63 0 : V(1) = TERM1
64 : ! -------- I=2,..,N-1 -----------
65 0 : do I = 2, N - 1
66 0 : TERM1 = (TH(I - 1) - 2.D0*TH(I) + TH(I + 1))*NQUATR
67 0 : V(I) = TERM1
68 : end do
69 : ! ----------- I=N -------------
70 0 : TERM1 = (TH(N - 1) - TH(N))*NQUATR
71 0 : V(N) = TERM1
72 : else
73 : ! --------- CASE T LESS EQUAL PI ------------
74 0 : FABS = 1.5D0*sin(T)*sin(T)
75 0 : FX = -FABS
76 0 : FY = FABS
77 : ! ---------- I=1 ------------
78 0 : TERM1 = (-3.D0*TH(1) + TH(2))*NQUATR
79 0 : TERM2 = NSQ*(FY*cos(TH(1)) - FX*sin(TH(1)))
80 0 : V(1) = TERM1 + TERM2
81 : ! -------- I=2,..,N-1 -----------
82 0 : do I = 2, N - 1
83 0 : TERM1 = (TH(I - 1) - 2.D0*TH(I) + TH(I + 1))*NQUATR
84 0 : TERM2 = NSQ*(FY*cos(TH(I)) - FX*sin(TH(I)))
85 0 : V(I) = TERM1 + TERM2
86 : end do
87 : ! ----------- I=N -------------
88 0 : TERM1 = (TH(N - 1) - TH(N))*NQUATR
89 0 : TERM2 = NSQ*(FY*cos(TH(N)) - FX*sin(TH(N)))
90 0 : V(N) = TERM1 + TERM2
91 : end if
92 : ! -------- COMPUTE PRODUCT DV=W ------------
93 0 : W(1) = STH(2)*V(2)
94 0 : do I = 2, N - 1
95 0 : W(I) = -STH(I)*V(I - 1) + STH(I + 1)*V(I + 1)
96 : end do
97 0 : W(N) = -STH(N)*V(N - 1)
98 : ! -------- TERME 3 -----------------
99 0 : do I = 1, N
100 0 : W(I) = W(I) + TH(N + I)*TH(N + I)
101 : end do
102 : ! ------- SOLVE SYSTEM CW=W ---------
103 0 : ALPHA(1) = 1.D0
104 0 : do I = 2, N
105 0 : ALPHA(I) = 2.D0
106 0 : BETA(I - 1) = -CTH(I)
107 : end do
108 0 : ALPHA(N) = 3.D0
109 0 : do I = N - 1, 1, -1
110 0 : Q = BETA(I)/ALPHA(I + 1)
111 0 : W(I) = W(I) - W(I + 1)*Q
112 0 : ALPHA(I) = ALPHA(I) - BETA(I)*Q
113 : end do
114 0 : W(1) = W(1)/ALPHA(1)
115 0 : do I = 2, N
116 0 : W(I) = (W(I) - BETA(I - 1)*W(I - 1))/ALPHA(I)
117 : end do
118 : ! -------- COMPUTE U=CV+DW ---------
119 0 : U(1) = V(1) - CTH(2)*V(2) + STH(2)*W(2)
120 0 : do I = 2, N - 1
121 0 : U(I) = 2.D0*V(I) - CTH(I)*V(I - 1) - CTH(I + 1)*V(I + 1) - STH(I)*W(I - 1) + STH(I + 1)*W(I + 1)
122 : end do
123 0 : U(N) = 3.D0*V(N) - CTH(N)*V(N - 1) - STH(N)*W(N - 1)
124 : ! -------- PUT DERIVATIVES IN RIGHT PLACE -------------
125 0 : do I = 1, N
126 0 : DF(I) = TH(N + I)
127 0 : DF(N + I) = U(I)
128 : end do
129 0 : end subroutine beam_feval
130 : ! ----------------------------------------------------------------------
131 0 : subroutine beam_jeval(ldim, neqn, t, y, yprime, dfdy, ierr, rpar, ipar)
132 : integer :: ldim, neqn, ierr, ipar(*)
133 : real(dp) :: t, y(neqn), yprime(neqn), dfdy(ldim, neqn), rpar(*)
134 : !
135 : ! dummy subroutine
136 : !
137 0 : end subroutine beam_jeval
138 : ! ----------------------------------------------------------------------
139 0 : subroutine beam_solut(neqn, t, y)
140 : integer :: neqn
141 : real(dp), intent(in) :: t
142 : real(dp), intent(out) :: y(neqn)
143 :
144 : ! computed using real(dp) RADAU on an
145 : ! Alphaserver DS20E, with a 667 MHz EV67 processor.
146 : !
147 : ! uround = 1.01d-19
148 : ! rtol = atol = h0 = 1.1d-18
149 :
150 0 : y(1) = -0.5792366591285007D-002
151 0 : y(2) = -0.1695298550721735D-001
152 0 : y(3) = -0.2769103312973094D-001
153 0 : y(4) = -0.3800815655879158D-001
154 0 : y(5) = -0.4790616859743763D-001
155 0 : y(6) = -0.5738710435274594D-001
156 0 : y(7) = -0.6645327313454617D-001
157 0 : y(8) = -0.7510730581979037D-001
158 0 : y(9) = -0.8335219765414992D-001
159 0 : y(10) = -0.9119134654647947D-001
160 0 : y(11) = -0.9862858700132091D-001
161 0 : y(12) = -0.1056682200378002D+000
162 0 : y(13) = -0.1123150395409595D+000
163 0 : y(14) = -0.1185743552727245D+000
164 0 : y(15) = -0.1244520128755561D+000
165 0 : y(16) = -0.1299544113264161D+000
166 0 : y(17) = -0.1350885180610398D+000
167 0 : y(18) = -0.1398618819194422D+000
168 0 : y(19) = -0.1442826441015242D+000
169 0 : y(20) = -0.1483595472463012D+000
170 0 : y(21) = -0.1521019429001447D+000
171 0 : y(22) = -0.1555197978061129D+000
172 0 : y(23) = -0.1586236993420229D+000
173 0 : y(24) = -0.1614248603702127D+000
174 0 : y(25) = -0.1639351238193223D+000
175 0 : y(26) = -0.1661669673440852D+000
176 0 : y(27) = -0.1681335081778558D+000
177 0 : y(28) = -0.1698485080602243D+000
178 0 : y(29) = -0.1713263782440888D+000
179 0 : y(30) = -0.1725821847462537D+000
180 0 : y(31) = -0.1736316537975654D+000
181 0 : y(32) = -0.1744911773840049D+000
182 0 : y(33) = -0.1751778187863392D+000
183 0 : y(34) = -0.1757093178712902D+000
184 0 : y(35) = -0.1761040960228576D+000
185 0 : y(36) = -0.1763812607175549D+000
186 0 : y(37) = -0.1765606097564671D+000
187 0 : y(38) = -0.1766626352260565D+000
188 0 : y(39) = -0.1767085270807460D+000
189 0 : y(40) = -0.1767201761075488D+000
190 0 : y(41) = 0.3747362681329794D-001
191 0 : y(42) = 0.1099117880217593D+000
192 0 : y(43) = 0.1798360474312799D+000
193 0 : y(44) = 0.2472427305442391D+000
194 0 : y(45) = 0.3121293820596567D+000
195 0 : y(46) = 0.3744947377019500D+000
196 0 : y(47) = 0.4343386073492798D+000
197 0 : y(48) = 0.4916620354601748D+000
198 0 : y(49) = 0.5464677854586807D+000
199 0 : y(50) = 0.5987609702624270D+000
200 0 : y(51) = 0.6485493611110740D+000
201 0 : y(52) = 0.6958435169132503D+000
202 0 : y(53) = 0.7406572668520808D+000
203 0 : y(54) = 0.7830081747813241D+000
204 0 : y(55) = 0.8229176659201515D+000
205 0 : y(56) = 0.8604110305190560D+000
206 0 : y(57) = 0.8955175502377805D+000
207 0 : y(58) = 0.9282708263127953D+000
208 0 : y(59) = 0.9587089334522034D+000
209 0 : y(60) = 0.9868747821728363D+000
210 0 : y(61) = 0.1012816579961883D+001
211 0 : y(62) = 0.1036587736679858D+001
212 0 : y(63) = 0.1058246826481355D+001
213 0 : y(64) = 0.1077857811433353D+001
214 0 : y(65) = 0.1095490222005369D+001
215 0 : y(66) = 0.1111219164294898D+001
216 0 : y(67) = 0.1125125269286501D+001
217 0 : y(68) = 0.1137294526609229D+001
218 0 : y(69) = 0.1147818025153607D+001
219 0 : y(70) = 0.1156792132004482D+001
220 0 : y(71) = 0.1164318845130183D+001
221 0 : y(72) = 0.1170505992596124D+001
222 0 : y(73) = 0.1175467424299550D+001
223 0 : y(74) = 0.1179323003228859D+001
224 0 : y(75) = 0.1182198586299667D+001
225 0 : y(76) = 0.1184226111223146D+001
226 0 : y(77) = 0.1185543909805575D+001
227 0 : y(78) = 0.1186297084203716D+001
228 0 : y(79) = 0.1186637618908127D+001
229 0 : y(80) = 0.1186724615113034D+001
230 :
231 0 : end subroutine beam_solut
232 :
233 : end module bari_beam
|
