0ccf6c01b1f648b5251f2e0b0d2e007dc4f44ec8
[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / dassldasrt.dia.ref
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7 // <-- CLI SHELL MODE -->
8 //DASSL
9 // PROBLEM 1..   LINEAR DIFFERENTIAL/ALGEBRAIC SYSTEM
10 //
11 //X1DOT + 10.0*X1 = 0  
12 //X1 + X2 = 1
13 //X1(0) = 1.0, X2(0) = 0.0
14 //
15 t=1:10;t0=0;y0=[1;0];y0d=[-10;0];
16 info=list([],0,[],[],[],0,0);
17 //    Calling Scilab functions
18 deff('[r,ires]=dres1(t,y,ydot)','r=[ydot(1)+10*y(1);y(2)+y(1)-1];ires=0')
19 deff('pd=djac1(t,y,ydot,cj)','pd=[cj+10,0;1,1]')
20 //   scilab function, without jacobian
21 yy0=dassl([y0,y0d],t0,t,dres1,info);
22 //   scilab functions, with jacobian
23 yy1=dassl([y0,y0d],t0,t,dres1,djac1,info);
24 // fortran routine, without jacobian
25 yy2=dassl([y0,y0d],t0,t,'dres1',info);   //=yy0
26 assert_checkalmostequal(norm(yy2,1),norm(yy0,1),1E-5);
27 // fortran routines, with jacobian
28 yy3=dassl([y0,y0d],t0,t,'dres1','djac1',info);  //=yy1
29 assert_checkalmostequal(norm(yy3,1),norm(yy1,1),1E-5);
30 yy3bis=dassl([y0,y0d],t0,t,'dres1',djac1,info);
31 // call fortran dres1 and scilab's djac1
32 yy3ter=dassl([y0,y0d],t0,t,dres1,'djac1',info);
33 //
34 // with specific atol and rtol parameters
35 atol=1.d-6;rtol=0;
36 yy4=dassl([y0,y0d],t0,t,atol,rtol,dres1,info);
37 yy5=dassl([y0,y0d],t0,t,atol,rtol,'dres1',info); //=yy4
38 assert_checkalmostequal(norm(yy5,1),norm(yy4,1),1E-9);
39 yy6=dassl([y0,y0d],t0,t,atol,rtol,dres1,djac1,info);
40 yy7=dassl([y0,y0d],t0,t,atol,rtol,'dres1','djac1',info); //==yy6
41 assert_checkalmostequal(norm(yy7,1),norm(yy6,1),1E-12);
42 //    
43 //   Testing E xdot - A x=0
44 //   x(0)=x0;   xdot(0)=xd0
45 rand('seed',0);
46 nx=5;
47 E=rand(nx,1)*rand(1,nx);A=rand(nx,nx);
48 //         Index 1
49 [Si,Pi,Di,o]=penlaur(E,A);pp=Si*E;[q,m]=fullrf(pp);x0=q(:,1);x0d=pinv(E)*A*x0;
50 deff('[r,ires]=g(t,x,xdot)','r=E*xdot-A*x;ires=0');
51 t=[1,2,3];t0=0;info=list([],0,[],[],[],0,0);
52 x=dassl([x0,x0d],t0,t,g,info);x1=x(2:nx+1,:);
53 assert_checkalmostequal(norm(pp*x1,1),norm(x1,1),1.d-5);
54 //x(4) goes through 1 at  t=1.5409711;
55 t=1.5409711;ww=dassl([x0,x0d],t0,t,g,info);
56 assert_checkalmostequal(ww(5),1,0.001);
57 deff('[rt]=surface(t,y,yd)','rt=y(4)-1');nbsurf=1;
58 [yyy,nnn]=dasrt([x0,x0d],t0,t,g,nbsurf,surface,info);
59 deff('pd=j(t,y,ydot,cj)','pd=cj*E-A');
60 x=dassl([x0,x0d],t0,t,g,j,info);x2=x(2:nx+1,1);
61 assert_checkalmostequal(norm(x2,1),norm(ww(2:nx+1,1),1),0.0001);
62 [yyy1,nnn]=dasrt([x0,x0d],t0,t,g,j,nbsurf,surface,info);
63 //x0d is not known:
64 x0d=ones(x0);info(7)=1;
65 x=dassl([x0,x0d],t0,t,g,info);
66 xn=dassl([x0,x0d],t0,t,g,j,info);
67 assert_checkalmostequal(norm(x,1),norm(xn,1),0.00001);
68 //PROBLEM 2..
69 info=list([],0,[],[],[],0,0);
70 y0=zeros(25,1);y0(1)=1;
71 delta=0*y0;
72 //link('dres2.o','dres2');
73 //y0d=call('dres2',0,1,'d',y0,2,'d',delta,3,'d',0,5,'i',0,6,'d',0,7,'d','out',[25,1],4,'d');
74 y0d=zeros(y0);y0d(1)=-2;y0d(2)=1;y0d(6)=1;
75 t0=0;t=[0.01,0.1,1,10,100];
76 rtol=0;atol=1.d-6;
77 y=dassl([y0,y0d],t0,t,atol,rtol,'dres2',info);
78 //                 DASRT
79 // 
80 //C-----------------------------------------------------------------------
81 //C First problem.
82 //C The initial value problem is..
83 //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
84 //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
85 //C The two root functions are..
86 //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
87 //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
88 //C-----------------------------------------------------------------------
89 y0=1;t=2:6;t0=1;y0d=3;
90 info=list([],0,[],[],[],0,0);
91 atol=1.d-6;rtol=0;ng=2;
92 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
93 assert_checkalmostequal(nn(1),2.47,0.001);
94 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
95 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
96 assert_checkalmostequal(nn(1),2.5,0.001);
97 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
98 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
99 assert_checkalmostequal(nn(1),2.53,0.001);
100 deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
101 deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
102 y0=1;t=2:6;t0=1;y0d=3;
103 info=list([],0,[],[],[],0,0);
104 atol=1.d-6;rtol=0;ng=2;
105 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
106 assert_checkalmostequal(nn(1),2.47,0.001);
107 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
108 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
109 assert_checkalmostequal(nn(1),2.5,0.001);
110 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
111 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
112 assert_checkalmostequal(nn(1),2.53,0.001);
113 //C
114 //C-----------------------------------------------------------------------
115 //C Second problem (Van Der Pol oscillator).
116 //C The initial value problem is..
117 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
118 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
119 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
120 //C The root function is  G = Y1.
121 //C An analytic solution is not known, but the zeros of Y1 are known
122 //C to 15 figures for purposes of checking the accuracy.
123 //C-----------------------------------------------------------------------
124 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
125 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
126 info=list([],0,[],[],[],0,0);
127 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
128 assert_checkalmostequal(nn(1),81.163512,0.001);
129 deff('[delta,ires]=res2(t,y,ydot)',...
130 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
131 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
132 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)]')
133 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
134 deff('s=gr2(t,y,yd)','s=y(1)')
135 [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
136 //           Hot Restart
137 [yy,nn,hotd]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
138 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
139 [yy,nn,hotd]=dasrt([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hotd);
140 assert_checkalmostequal(nn(1),162.57763,0.004);
141 //Test of Dynamic link (Require f77!)
142 //         1 making the routines
143 res22=[...
144 '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
145 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
146 '      INTEGER NEQ';
147 '      DIMENSION Y(*), YDOT(*), DELTA(*)';
148 '      NEQ=2';
149 '      CALL F2(NEQ,T,Y,DELTA)';
150 '      DO 10 I = 1,NEQ';
151 '         DELTA(I) = YDOT(I) - DELTA(I)';
152 ' 10   CONTINUE';
153 '      RETURN';
154 '      END';
155 '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
156 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
157 '      INTEGER NEQ';
158 '      DOUBLE PRECISION T, Y, YDOT';
159 '      DIMENSION Y(*), YDOT(*)';
160 '      YDOT(1) = Y(2)';
161 '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
162 '      RETURN';
163 '      END';]
164  res22  = 
165 !      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)    !
166 !                                                         !
167 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
168 !                                                         !
169 !      INTEGER NEQ                                        !
170 !                                                         !
171 !      DIMENSION Y(*), YDOT(*), DELTA(*)                  !
172 !                                                         !
173 !      NEQ=2                                              !
174 !                                                         !
175 !      CALL F2(NEQ,T,Y,DELTA)                             !
176 !                                                         !
177 !      DO 10 I = 1,NEQ                                    !
178 !                                                         !
179 !         DELTA(I) = YDOT(I) - DELTA(I)                   !
180 !                                                         !
181 ! 10   CONTINUE                                           !
182 !                                                         !
183 !      RETURN                                             !
184 !                                                         !
185 !      END                                                !
186 !                                                         !
187 !      SUBROUTINE F2 (NEQ, T, Y, YDOT)                    !
188 !                                                         !
189 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
190 !                                                         !
191 !      INTEGER NEQ                                        !
192 !                                                         !
193 !      DOUBLE PRECISION T, Y, YDOT                        !
194 !                                                         !
195 !      DIMENSION Y(*), YDOT(*)                            !
196 !                                                         !
197 !      YDOT(1) = Y(2)                                     !
198 !                                                         !
199 !      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)  !
200 !                                                         !
201 !      RETURN                                             !
202 !                                                         !
203 !      END                                                !
204 jac22=[...
205 '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
206 ' ';
207 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
208 '      INTEGER  NROWPD';
209 '      DOUBLE PRECISION T, Y, PD';
210 '      PARAMETER (NROWPD=2)';
211 '      DIMENSION Y(2), PD(NROWPD,2)';
212 '      PD(1,1) = 0.0D0';
213 '      PD(1,2) = 1.0D0';
214 '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
215 '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
216 '      PD(1,1) = CJ - PD(1,1)';
217 '      PD(1,2) =    - PD(1,2)';
218 '      PD(2,1) =    - PD(2,1)';
219 '      PD(2,2) = CJ - PD(2,2)';
220 '      RETURN';
221 '      END';]
222  jac22  = 
223 !      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)  !
224 !                                                         !
225 !                                                         !
226 !                                                         !
227 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
228 !                                                         !
229 !      INTEGER  NROWPD                                    !
230 !                                                         !
231 !      DOUBLE PRECISION T, Y, PD                          !
232 !                                                         !
233 !      PARAMETER (NROWPD=2)                               !
234 !                                                         !
235 !      DIMENSION Y(2), PD(NROWPD,2)                       !
236 !                                                         !
237 !      PD(1,1) = 0.0D0                                    !
238 !                                                         !
239 !      PD(1,2) = 1.0D0                                    !
240 !                                                         !
241 !      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0               !
242 !                                                         !
243 !      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))              !
244 !                                                         !
245 !      PD(1,1) = CJ - PD(1,1)                             !
246 !                                                         !
247 !      PD(1,2) =    - PD(1,2)                             !
248 !                                                         !
249 !      PD(2,1) =    - PD(2,1)                             !
250 !                                                         !
251 !      PD(2,2) = CJ - PD(2,2)                             !
252 !                                                         !
253 !      RETURN                                             !
254 !                                                         !
255 !      END                                                !
256 gr22=[...
257 '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
258 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
259 '      INTEGER NEQ, NG';
260 '      DOUBLE PRECISION T, Y, GROOT';
261 '      DIMENSION Y(*), GROOT(*)';
262 '      GROOT(1) = Y(1)';
263 '      RETURN';
264 '      END';]
265  gr22  = 
266 !      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)  !
267 !                                                          !
268 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                 !
269 !                                                          !
270 !      INTEGER NEQ, NG                                     !
271 !                                                          !
272 !      DOUBLE PRECISION T, Y, GROOT                        !
273 !                                                          !
274 !      DIMENSION Y(*), GROOT(*)                            !
275 !                                                          !
276 !      GROOT(1) = Y(1)                                     !
277 !                                                          !
278 !      RETURN                                              !
279 !                                                          !
280 !      END                                                 !
281 //Uncomment lines below: link may be machine dependent if some f77 libs are 
282 //needed for linking
283 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
284 //unix_g('cd /tmp;rm -f /tmp/jac22.f');
285 //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
286 //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
287 //unix_g('cd /tmp;make /tmp/jac22.o');
288 //          2  Linking the routines
289 //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
290 //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
291 //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
292 //info=list([],0,[],[],[],0,0);
293 //          3 Calling the routines by dasrt
294 //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
295 // hot restart
296 //[yy,nn,hotd]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
297 //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
298 //[yy,nn,hotd]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hotd);