scilab/modules/differential_equations/tests/unit_tests/dassldasrt.dia.ref

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