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