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 // Execute with exec("SCI/modules/differential_equations/tests/unit_tests/daskr.tst");
  12 //  or test_run('differential_equations', 'daskr', ['no_check_error_output']);
  13 //C-----------------------------------------------------------------------
  14 //C First problem.
  15 //C The initial value problem is..
  16 //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
  17 //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
  18 //C The two root functions are..
  19 //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
  20 //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
  21 //C-----------------------------------------------------------------------
  22 y0=1;t=2:6;t0=1;y0d=3;
  23 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
  24 atol=1.d-6;rtol=0;ng=2;
  25 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
  26 if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
  27 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
  28 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
  29 if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
  30 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
  31 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
  32 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  33 if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
  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 //C
  47 //C-----------------------------------------------------------------------
  48 //C Second problem (Van Der Pol oscillator).
  49 //C The initial value problem is..
  50 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
  51 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
  52 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
  53 //C The root function is  G = Y1.
  54 //C An analytic solution is not known, but the zeros of Y1 are known
  55 //C to 15 figures for purposes of checking the accuracy.
  56 //C-----------------------------------------------------------------------
  57 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
  58 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
  59 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
  60 if abs(nn(1)-81.163512)>0.001 then bugmes();quit;end
  61 deff('[delta,ires]=res2(t,y,ydot)',...
  62 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
  63 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
  64 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)]')
  65 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
  66 deff('s=gr2(t,y,yd)','s=y(1)')
  67 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
  68 //           Hot Restart
  69 [yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
  70 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
  71 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hotd);
  72 if abs(nn(1)-162.57763)>0.004 then bugmes();quit;end
  73 //Test of Dynamic link (Require f77!)
  74 //         1 making the routines
  75 res22=[...
  76 '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
  77 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  78 '      INTEGER NEQ';
  79 '      DIMENSION Y(*), YDOT(*), DELTA(*)';
  80 '      NEQ=2';
  81 '      CALL F2(NEQ,T,Y,DELTA)';
  82 '      DO 10 I = 1,NEQ';
  83 '         DELTA(I) = YDOT(I) - DELTA(I)';
  84 ' 10   CONTINUE';
  85 '      RETURN';
  86 '      END';
  87 '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
  88 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  89 '      INTEGER NEQ';
  90 '      DOUBLE PRECISION T, Y, YDOT';
  91 '      DIMENSION Y(*), YDOT(*)';
  92 '      YDOT(1) = Y(2)';
  93 '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
  94 '      RETURN';
  95 '      END';]
  96  res22  =
  97
  98 !      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)    !
  99 !                                                         !
 100 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
 101 !                                                         !
 102 !      INTEGER NEQ                                        !
 103 !                                                         !
 104 !      DIMENSION Y(*), YDOT(*), DELTA(*)                  !
 105 !                                                         !
 106 !      NEQ=2                                              !
 107 !                                                         !
 108 !      CALL F2(NEQ,T,Y,DELTA)                             !
 109 !                                                         !
 110 !      DO 10 I = 1,NEQ                                    !
 111 !                                                         !
 112 !         DELTA(I) = YDOT(I) - DELTA(I)                   !
 113 !                                                         !
 114 ! 10   CONTINUE                                           !
 115 !                                                         !
 116 !      RETURN                                             !
 117 !                                                         !
 118 !      END                                                !
 119 !                                                         !
 120 !      SUBROUTINE F2 (NEQ, T, Y, YDOT)                    !
 121 !                                                         !
 122 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
 123 !                                                         !
 124 !      INTEGER NEQ                                        !
 125 !                                                         !
 126 !      DOUBLE PRECISION T, Y, YDOT                        !
 127 !                                                         !
 128 !      DIMENSION Y(*), YDOT(*)                            !
 129 !                                                         !
 130 !      YDOT(1) = Y(2)                                     !
 131 !                                                         !
 132 !      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)  !
 133 !                                                         !
 134 !      RETURN                                             !
 135 !                                                         !
 136 !      END                                                !
 137 jac22=[...
 138 '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
 139 ' ';
 140 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
 141 '      INTEGER  NROWPD';
 142 '      DOUBLE PRECISION T, Y, PD';
 143 '      PARAMETER (NROWPD=2)';
 144 '      DIMENSION Y(2), PD(NROWPD,2)';
 145 '      PD(1,1) = 0.0D0';
 146 '      PD(1,2) = 1.0D0';
 147 '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
 148 '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
 149 '      PD(1,1) = CJ - PD(1,1)';
 150 '      PD(1,2) =    - PD(1,2)';
 151 '      PD(2,1) =    - PD(2,1)';
 152 '      PD(2,2) = CJ - PD(2,2)';
 153 '      RETURN';
 154 '      END';]
 155  jac22  =
 156
 157 !      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)  !
 158 !                                                         !
 159 !                                                         !
 160 !                                                         !
 161 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
 162 !                                                         !
 163 !      INTEGER  NROWPD                                    !
 164 !                                                         !
 165 !      DOUBLE PRECISION T, Y, PD                          !
 166 !                                                         !
 167 !      PARAMETER (NROWPD=2)                               !
 168 !                                                         !
 169 !      DIMENSION Y(2), PD(NROWPD,2)                       !
 170 !                                                         !
 171 !      PD(1,1) = 0.0D0                                    !
 172 !                                                         !
 173 !      PD(1,2) = 1.0D0                                    !
 174 !                                                         !
 175 !      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0               !
 176 !                                                         !
 177 !      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))              !
 178 !                                                         !
 179 !      PD(1,1) = CJ - PD(1,1)                             !
 180 !                                                         !
 181 !      PD(1,2) =    - PD(1,2)                             !
 182 !                                                         !
 183 !      PD(2,1) =    - PD(2,1)                             !
 184 !                                                         !
 185 !      PD(2,2) = CJ - PD(2,2)                             !
 186 !                                                         !
 187 !      RETURN                                             !
 188 !                                                         !
 189 !      END                                                !
 190 gr22=[...
 191 '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
 192 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
 193 '      INTEGER NEQ, NG';
 194 '      DOUBLE PRECISION T, Y, GROOT';
 195 '      DIMENSION Y(*), GROOT(*)';
 196 '      GROOT(1) = Y(1)';
 197 '      RETURN';
 198 '      END';]
 199  gr22  =
 200
 201 !      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)  !
 202 !                                                          !
 203 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                 !
 204 !                                                          !
 205 !      INTEGER NEQ, NG                                     !
 206 !                                                          !
 207 !      DOUBLE PRECISION T, Y, GROOT                        !
 208 !                                                          !
 209 !      DIMENSION Y(*), GROOT(*)                            !
 210 !                                                          !
 211 !      GROOT(1) = Y(1)                                     !
 212 !                                                          !
 213 !      RETURN                                              !
 214 !                                                          !
 215 !      END                                                 !
 216 //Uncomment lines below: link may be machine dependent if some f77 libs are
 217 //needed for linking
 218 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
 219 //unix_g('cd /tmp;rm -f /tmp/jac22.f');
 220 //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
 221 //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
 222 //unix_g('cd /tmp;make /tmp/jac22.o');
 223 //          2  Linking the routines
 224 //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
 225 //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
 226 //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 227 //          3 Calling the routines by dasrt
 228 //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
 229 // hot restart
 230 //[yy,nn,hotd]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
 231 //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
 232 //[yy,nn,hotd]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hotd);