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