edee5f2b333ad5cfe93a390335db1103a49fdb01
[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / daskr.dia.ref
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING,
7 // which you should have received as part of this distribution.
8 // The terms are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 // =============================================================================
11 // Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
12 //C-----------------------------------------------------------------------
13 //C First problem.
14 //C The initial value problem is..
15 //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
16 //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
17 //C The two root functions are..
18 //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
19 //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
20 //C-----------------------------------------------------------------------
21 y0=1;t=2:6;t0=1;y0d=3;
22 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
23 atol=1.d-6;rtol=0;ng=2;
24 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
25 if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
26 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
27 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
28 if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
29 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
30 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
31 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
32 if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
33 // Same problem, but using macro for the derivative evaluation function 'res1'
34 deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2.*log(y)+8)./t-5).*y')
35 deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
36 y0=1;t=2:6;t0=1;y0d=3;
37 atol=1.d-6;rtol=0;ng=2;
38 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
39 if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
40 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
41 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
42 if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
43 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
44 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
45 if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
46 // Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
47 // pjac uses the macro res1 defined above.
48 function [r, ier] = psol(wp, iwp, b)
49     ier = 0;
50     //Compute the LU factorization of R.
51     sp = sparse(iwp, wp);
52     [h, rk] = lufact(sp);
53     //Solve the system LU*X = b
54     r = lusolve(h, b);
55     ludel(h);
56 endfunction
57 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
58     ires = 0;
59     SQuround = 1.490D-08;
60     tx = t;
61     nrow = 0;
62     e = zeros(1, neq);
63     wp = zeros(neq*neq, 1);
64     iwp = zeros(neq*neq, 2);
65     for i=1:neq
66         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
67         if h*ydot(i) < 0 then del = -del; end
68         ysave = y(i);
69         ypsave = ydot(i);
70         y(i) = y(i) + del;
71         ydot(i) = ydot(i) + cj*del;
72         [e ires]=res1(tx, y, ydot);
73         if ires < 0 then return; end
74         delinv = 1/del;
75         for j=1:neq
76             wp(nrow+j) = delinv*(e(j)-savr(j));
77             iwp(nrow+j,1) = i;
78             iwp(nrow+j,2) = j;
79         end
80         nrow = nrow + neq;
81         y(i) = ysave;
82         ydot(i) = ypsave;
83     end
84 endfunction
85 y0=1;t=2:6;t0=1;y0d=3;
86 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
87 atol=1.d-6;rtol=0;ng=2;
88 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
89 if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
90 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
91 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
92 if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
93 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
94 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
95 if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
96 //C
97 //C-----------------------------------------------------------------------
98 //C Second problem (Van Der Pol oscillator).
99 //C The initial value problem is..
100 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
101 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
102 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
103 //C The root function is  G = Y1.
104 //C An analytic solution is not known, but the zeros of Y1 are known
105 //C to 15 figures for purposes of checking the accuracy.
106 //C-----------------------------------------------------------------------
107 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
108 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
109 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
110 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
111 if abs(nn(1)-81.163512)>0.001 then bugmes();quit;end
112 deff('[delta,ires]=res2(t,y,ydot)',...
113 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
114 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
115 deff('J=jac2(t,y,ydot,c)','y1=y(1);y2=y(2);J=[c,-1;200*y1*y2+1,c-100*(1-y1*y1)]')
116 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
117 deff('s=gr2(t,y,yd)','s=y(1)')
118 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
119 // Same problem, with psol and pjac example routines
120 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
121 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,'psol1','pjac1');
122 if abs(nn(1)-81.163512)>0.009 then bugmes();quit;end
123 deff('s=gr2(t,y,yd)','s=y(1)')
124 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,'psol1','pjac1');
125 if abs(nn(1)-81.163512)>0.009 then bugmes();quit;end
126 // Same problem, with psol and pjac macros
127 // Redefine pjac to use res2
128 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
129     ires = 0;
130     SQuround = 1.490D-08;
131     tx = t;
132     nrow = 0;
133     e = zeros(1, neq);
134     wp = zeros(neq*neq, 1);
135     iwp = zeros(neq*neq, 2);
136     for i=1:neq
137         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
138         if h*ydot(i) < 0 then del = -del; end
139         ysave = y(i);
140         ypsave = ydot(i);
141         y(i) = y(i) + del;
142         ydot(i) = ydot(i) + cj*del;
143         [e ires]=res2(tx, y, ydot);
144         if ires < 0 then return; end
145         delinv = 1/del;
146         for j=1:neq
147             wp(nrow+j) = delinv*(e(j)-savr(j));
148             iwp(nrow+j,1) = i;
149             iwp(nrow+j,2) = j;
150         end
151         nrow = nrow + neq;
152         y(i) = ysave;
153         ydot(i) = ypsave;
154     end
155 endfunction
156 Warning : redefining function: pjac                    . Use funcprot(0) to avoid this message
157
158 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,psol,pjac);
159 if abs(nn(1)-81.163512)>0.003 then bugmes();quit;end
160 deff('s=gr2(t,y,yd)','s=y(1)')
161 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,psol,pjac);
162 if abs(nn(1)-81.163512)>0.003 then bugmes();quit;end
163 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
164 //           Hot Restart
165 [yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
166 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
167 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hotd);
168 if abs(nn(1)-162.57763)>0.004 then bugmes();quit;end
169 //Test of Dynamic link (Require f77!)
170 //         1 making the routines
171 res22=[...
172 '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
173 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
174 '      INTEGER NEQ';
175 '      DIMENSION Y(*), YDOT(*), DELTA(*)';
176 '      NEQ=2';
177 '      CALL F2(NEQ,T,Y,DELTA)';
178 '      DO 10 I = 1,NEQ';
179 '         DELTA(I) = YDOT(I) - DELTA(I)';
180 ' 10   CONTINUE';
181 '      RETURN';
182 '      END';
183 '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
184 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
185 '      INTEGER NEQ';
186 '      DOUBLE PRECISION T, Y, YDOT';
187 '      DIMENSION Y(*), YDOT(*)';
188 '      YDOT(1) = Y(2)';
189 '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
190 '      RETURN';
191 '      END';]
192  res22  =
193  
194 !      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)    !
195 !                                                         !
196 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
197 !                                                         !
198 !      INTEGER NEQ                                        !
199 !                                                         !
200 !      DIMENSION Y(*), YDOT(*), DELTA(*)                  !
201 !                                                         !
202 !      NEQ=2                                              !
203 !                                                         !
204 !      CALL F2(NEQ,T,Y,DELTA)                             !
205 !                                                         !
206 !      DO 10 I = 1,NEQ                                    !
207 !                                                         !
208 !         DELTA(I) = YDOT(I) - DELTA(I)                   !
209 !                                                         !
210 ! 10   CONTINUE                                           !
211 !                                                         !
212 !      RETURN                                             !
213 !                                                         !
214 !      END                                                !
215 !                                                         !
216 !      SUBROUTINE F2 (NEQ, T, Y, YDOT)                    !
217 !                                                         !
218 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
219 !                                                         !
220 !      INTEGER NEQ                                        !
221 !                                                         !
222 !      DOUBLE PRECISION T, Y, YDOT                        !
223 !                                                         !
224 !      DIMENSION Y(*), YDOT(*)                            !
225 !                                                         !
226 !      YDOT(1) = Y(2)                                     !
227 !                                                         !
228 !      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)  !
229 !                                                         !
230 !      RETURN                                             !
231 !                                                         !
232 !      END                                                !
233 jac22=[...
234 '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
235 ' ';
236 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
237 '      INTEGER  NROWPD';
238 '      DOUBLE PRECISION T, Y, PD';
239 '      PARAMETER (NROWPD=2)';
240 '      DIMENSION Y(2), PD(NROWPD,2)';
241 '      PD(1,1) = 0.0D0';
242 '      PD(1,2) = 1.0D0';
243 '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
244 '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
245 '      PD(1,1) = CJ - PD(1,1)';
246 '      PD(1,2) =    - PD(1,2)';
247 '      PD(2,1) =    - PD(2,1)';
248 '      PD(2,2) = CJ - PD(2,2)';
249 '      RETURN';
250 '      END';]
251  jac22  =
252  
253 !      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)  !
254 !                                                         !
255 !                                                         !
256 !                                                         !
257 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
258 !                                                         !
259 !      INTEGER  NROWPD                                    !
260 !                                                         !
261 !      DOUBLE PRECISION T, Y, PD                          !
262 !                                                         !
263 !      PARAMETER (NROWPD=2)                               !
264 !                                                         !
265 !      DIMENSION Y(2), PD(NROWPD,2)                       !
266 !                                                         !
267 !      PD(1,1) = 0.0D0                                    !
268 !                                                         !
269 !      PD(1,2) = 1.0D0                                    !
270 !                                                         !
271 !      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0               !
272 !                                                         !
273 !      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))              !
274 !                                                         !
275 !      PD(1,1) = CJ - PD(1,1)                             !
276 !                                                         !
277 !      PD(1,2) =    - PD(1,2)                             !
278 !                                                         !
279 !      PD(2,1) =    - PD(2,1)                             !
280 !                                                         !
281 !      PD(2,2) = CJ - PD(2,2)                             !
282 !                                                         !
283 !      RETURN                                             !
284 !                                                         !
285 !      END                                                !
286 gr22=[...
287 '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
288 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
289 '      INTEGER NEQ, NG';
290 '      DOUBLE PRECISION T, Y, GROOT';
291 '      DIMENSION Y(*), GROOT(*)';
292 '      GROOT(1) = Y(1)';
293 '      RETURN';
294 '      END';]
295  gr22  =
296  
297 !      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)  !
298 !                                                          !
299 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                 !
300 !                                                          !
301 !      INTEGER NEQ, NG                                     !
302 !                                                          !
303 !      DOUBLE PRECISION T, Y, GROOT                        !
304 !                                                          !
305 !      DIMENSION Y(*), GROOT(*)                            !
306 !                                                          !
307 !      GROOT(1) = Y(1)                                     !
308 !                                                          !
309 !      RETURN                                              !
310 !                                                          !
311 !      END                                                 !
312 //Uncomment lines below: link may be machine dependent if some f77 libs are
313 //needed for linking
314 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
315 //unix_g('cd /tmp;rm -f /tmp/jac22.f');
316 //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
317 //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
318 //unix_g('cd /tmp;make /tmp/jac22.o');
319 //          2  Linking the routines
320 //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
321 //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
322 //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
323 //          3 Calling the routines by dasrt
324 //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
325 // hot restart
326 //[yy,nn,hotd]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
327 //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
328 //[yy,nn,hotd]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hotd);