scilab/modules/differential_equations/tests/unit_tests/dae.tst

   1 // ===========================================================================
   2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 // Copyright (C) 2013 - Scilab Enterprises - Paul Bignier : added "root2" (daskr)
   4 // Copyright (C) 2008 - INRIA - Sabinere Gauzere
   5 //
   6 //  This file is distributed under the same license as the Scilab package.
   7 // ===========================================================================
   8
   9 // <-- CLI SHELL MODE -->
  10 // <-- ENGLISH IMPOSED -->
  11
  12 ilib_verbose(0);
  13
  14 //DASSL
  15 // PROBLEM 1..   LINEAR DIFFERENTIAL/ALGEBRAIC SYSTEM
  16 //
  17 //X1DOT + 10.0*X1 = 0
  18 //X1 + X2 = 1
  19 //X1(0) = 1.0, X2(0) = 0.0
  20 //
  21 t=1:10;t0=0;y0=[1;0];y0d=[-10;0];
  22 info=list([],0,[],[],[],0,0);
  23 //    Calling Scilab functions
  24 deff("[r,ires]=dres1(t,y,ydot)","r=[ydot(1)+10*y(1);y(2)+y(1)-1];ires=0")
  25 deff("pd=djac1(t,y,ydot,cj)","pd=[cj+10,0;1,1]")
  26 //   scilab function, without jacobian
  27 yy0=dae([y0,y0d],t0,t,dres1);
  28 //   scilab functions, with jacobian
  29 yy1=dae([y0,y0d],t0,t,dres1,djac1);
  30 // fortran routine, without jacobian
  31 yy2=dae([y0,y0d],t0,t,"dres1");   //=yy0
  32 if norm(yy2-yy0,1)>1E-5 then pause,end
  33 // fortran routines, with jacobian
  34 yy3=dae([y0,y0d],t0,t,"dres1","djac1");  //=yy1
  35 if norm(yy3-yy1,1)>1E-5 then pause,end
  36 yy3bis=dae([y0,y0d],t0,t,"dres1",djac1);
  37 // call fortran dres1 and scilab's djac1
  38 yy3ter=dae([y0,y0d],t0,t,dres1,"djac1");
  39 //
  40 // with specific atol and rtol parameters
  41 atol=1.d-6;rtol=0;
  42 yy4=dae([y0,y0d],t0,t,rtol,atol,dres1);
  43 yy5=dae([y0,y0d],t0,t,rtol,atol,"dres1"); //=yy4
  44 if norm(yy5-yy4,1)>1E-9 then pause,end
  45 yy6=dae([y0,y0d],t0,t,rtol,atol,dres1,djac1);
  46 yy7=dae([y0,y0d],t0,t,rtol,atol,"dres1","djac1"); //==yy6
  47 if norm(yy7-yy6,1)>1E-12 then pause,end
  48 //
  49 //   Testing E xdot - A x=0
  50 //   x(0)=x0;   xdot(0)=xd0
  51 rand("seed",0);
  52 nx=5;
  53 E=rand(nx,1)*rand(1,nx);
  54 A=rand(nx,nx);
  55 //         Index 1
  56 [Si,Pi,Di,o]=penlaur(E,A);
  57 pp=Si*E;
  58 [q,m]=fullrf(pp);
  59 x0=q(:,1);
  60 x0d=pinv(E)*A*x0;
  61 deff("[r,ires]=g(t,x,xdot)","r=E*xdot-A*x;ires=0");
  62 t=[1,2,3];t0=0;
  63 %DAEOPTIONS=list([],0,[],[],[],0,0);
  64 x=dae([x0,x0d],t0,t,g);
  65 x=[t;x];
  66 x1=x(2:nx+1,:);
  67 if norm(pp*x1-x1,1)>1.d-5 then pause,end // Bug because we don't have the first line anymore
  68 //x(4) goes through 1 at  t=1.5409711;
  69 %DAEOPTIONS=list([],0,[],[],[],0,0);
  70 t=1.5409711;
  71 ww=dae([x0,x0d],t0,t,g);
  72 ww=[t;ww];
  73 if abs(ww(5)-1)>0.001 then pause,end
  74 deff("[rt]=surface(t,y,yd)","rt=y(4)-1");
  75 nbsurf=1;
  76 [yyy,nnn]=dae("root",[x0,x0d],t0,t,g,nbsurf,surface);
  77 deff("pd=j(t,y,ydot,cj)","pd=cj*E-A");
  78 x=dae([x0,x0d],t0,t,g,j);
  79 x=[t;x];
  80 x2=x(2:nx+1,1);
  81 if norm(x2-ww(2:nx+1,1),1)>0.0001 then pause,end
  82 [yyy1,nnn]=dae("root",[x0,x0d],t0,t,g,j,nbsurf,surface);
  83 //x0d is not known:
  84 x0d=ones(x0);
  85 %DAEOPTIONS=list([],0,[],[],[],0,1);
  86 x=dae([x0,x0d],t0,t,g);
  87 xn=dae([x0,x0d],t0,t,g,j);
  88 if norm(x-xn,1)>0.00001 then pause,end
  89 //PROBLEM 2..
  90 DAEOPTIONS=list([],0,[],[],[],0,0);
  91 y0=zeros(25,1);y0(1)=1;
  92 delta=0*y0;
  93 //link('dres2.o','dres2');
  94 //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');
  95 y0d=zeros(y0);y0d(1)=-2;y0d(2)=1;y0d(6)=1;
  96 t0=0;t=[0.01,0.1,1,10,100];
  97 rtol=0;atol=1.d-6;
  98 y=dae([y0,y0d],t0,t,rtol,atol,"dres2");
  99
 100 //                 Root finder (dasrt)
 101 //
 102 //-----------------------------------------------------------------------
 103 // First problem.
 104 // The initial value problem is..
 105 //   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
 106 // The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
 107 // The two root functions are..
 108 //   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
 109 //   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
 110 //-----------------------------------------------------------------------
 111 y0=1;t=2:6;t0=1;y0d=3;
 112 %DAEOPTIONS=list([],0,[],[],[],0,0);
 113 atol=1.d-6;rtol=0;ng=2;
 114 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 115 if abs(nn(1)-2.47)>0.001 then pause,end
 116 y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
 117 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 118 if abs(nn(1)-2.5)>0.001 then pause,end
 119 y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
 120 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 121 if abs(nn(1)-2.53)>0.001 then pause,end
 122 deff("[delta,ires]=res1(t,y,ydot)","ires=0;delta=ydot-((2*log(y)+8)/t-5)*y")
 123 deff("[rts]=gr1(t,y,yd)","rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]")
 124 y0=1;t=2:6;t0=1;y0d=3;
 125 %DAEOPTIONS=list([],0,[],[],[],0,0);
 126 atol=1.d-6;rtol=0;ng=2;
 127 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
 128 if abs(nn(1)-2.47)>0.001 then pause,end
 129 y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
 130 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
 131 if abs(nn(1)-2.5)>0.001 then pause,end
 132 y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
 133 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
 134 if abs(nn(1)-2.53)>0.001 then pause,end
 135 //
 136 //-----------------------------------------------------------------------
 137 // Second problem (Van Der Pol oscillator).
 138 // The initial value problem is..
 139 //   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
 140 //   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
 141 //   Y1PRIME(0) = 0, Y2PRIME(0) = -2
 142 // The root function is  G = Y1.
 143 // An analytic solution is not known, but the zeros of Y1 are known
 144 // to 15 figures for purposes of checking the accuracy.
 145 //-----------------------------------------------------------------------
 146 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
 147 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 148 %DAEOPTIONS=list([],0,[],[],[],0,0);
 149 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
 150 if abs(nn(1)-81.163512)>0.001 then pause,end
 151 deff("[delta,ires]=res2(t,y,ydot)",...
 152 "ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]")
 153 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,"jac2",ng,"gr2");
 154 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)]")
 155 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2");
 156 deff("s=gr2(t,y,yd)","s=y(1)")
 157 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2);
 158 //           Hot Restart
 159 [yy,nn,hotd]=dae("root",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
 160 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(1:2,qq);y0d1=yy(2:3,qq);
 161 %DAEOPTIONS=list([],0,[],[],[],0,0);
 162 [yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,rtol,atol,"res2","jac2",ng,"gr2",hotd);
 163 if abs(nn(1)-162.57763)>0.004 then pause,end
 164
 165 //same with C code
 166 ilib_verbose(0);
 167
 168 cd TMPDIR;
 169 mkdir("dae_test1");
 170 cd("dae_test1");
 171
 172 code=["#include <math.h>"
 173 "void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)"
 174 "{res[0] = yd[0] - y[1];"
 175 " res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}"
 176 " "
 177 "void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)"
 178 "{pd[0]=*cj - 0.0;"
 179 " pd[1]=    - (-200.0*y[0]*y[1] - 1.0);"
 180 " pd[2]=    - 1.0;"
 181 " pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}"
 182 " "
 183 "void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)"
 184 "{ groot[0] = y[0];}"];
 185 mputl(code,"t22.c") ;
 186 ilib_for_link(["res22" "jac22" "gr22"],"t22.c","","c");
 187 exec("loader.sce");
 188
 189 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
 190 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 191 %DAEOPTIONS =list([],0,[],[],[],0,0);
 192 //hot restart
 193 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
 194 [yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,atol,rtol,"res22","jac22",ng,"gr22",hotd);
 195
 196 rtol=[1.d-6;1.d-6];
 197 atol=[1.d-6;1.d-4];
 198 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 199 %DAEOPTIONS =list([],0,[],[],[],0,0);
 200 [yy,nn]=dae("root",[y0,y0d],t0,t,atol,rtol,"res22","jac22",ng,"gr22");
 201 //hot restart
 202 [yy,nn,hotd]=dae("root",[y0,y0d],t0,t,atol,rtol,"res22","jac22",ng,"gr22");
 203 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
 204 [yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,atol,rtol,"res22","jac22",ng,"gr22",hotd);
 205
 206 //banded systems
 207 A=[-17,6,3,0,0,0,0,0,0,0;
 208 8,-12,9,4,0,0,0,0,0,0;
 209 0,7,-17,3,8,0,0,0,0,0;
 210 0,0,3,-13,2,2,0,0,0,0;
 211 0,0,0,4,-18,6,4,0,0,0;
 212 0,0,0,0,7,-13,7,0,0,0;
 213 0,0,0,0,0,7,-16,7,9,0;
 214 0,0,0,0,0,0,5,-17,8,1;
 215 0,0,0,0,0,0,0,4,-14,8;
 216 0,0,0,0,0,0,0,0,10,-10];
 217 n=size(A,1);
 218 //Full jacobian case for reference
 219 function [r,ires]=res(t,y,yd)
 220     r=yd-A*y; ires=0
 221 endfunction
 222 function pd=jac(x,y,yd,cj)
 223     pd=A+cj*eye()
 224 endfunction
 225
 226 y0=ones(n,1);yd0=0*y0;
 227 y=dae([y0,yd0],0,0:0.1:10,res,jac);
 228
 229 //banded estimated jacobian
 230 y1=dae([y0,yd0],0,0:0.1:10,res);
 231
 232 ml=1;mu=2;
 233 %DAEOPTIONS=list([],0,[ml,mu],[],[],0,0);
 234 yb1=dae([y0,yd0],0,0:0.1:10,res);
 235 norm(y1-yb1);
 236
 237 //banded  jacobian, C code
 238 //Residuals computation code
 239 cd TMPDIR;
 240 mkdir("dae_test2");
 241 cd("dae_test2");
 242 ccode=["#include <math.h>"
 243 "void myres(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)"
 244 "{"
 245 "  *ires =0;"
 246 "  res[0]=yd[0]+17.0*y[0]- 6.0*y[1]- 3.0*y[2];"
 247 "  res[1]=yd[1]-8.0*y[0]+12.0*y[1]- 9.0*y[2]- 4.0*y[3];"
 248 "  res[2]=yd[2]          -7.0*y[1]+17.0*y[2]- 3.0*y[3]- 8.0*y[4];"
 249 "  res[3]=yd[3]                    -3.0*y[2]+13.0*y[3]- 2.0*y[4]- 2.0*y[5];"
 250 "  res[4]=yd[4]                              -4.0*y[3]+18.0*y[4]- 6.0*y[5]- 4.0*y[6];"
 251 "  res[5]=yd[5]                                        -7.0*y[4]+13.0*y[5]- 7.0*y[6]- 0.0*y[7];"
 252 "  res[6]=yd[6]                                                  -7.0*y[5]+16.0*y[6]- 7.0*y[7]- 9.0*y[8];"
 253 "  res[7]=yd[7]                                                            -5.0*y[6]+17.0*y[7]- 8.0*y[8]- 1.0*y[9];"
 254 "  res[8]=yd[8]                                                                      -4.0*y[7]+14.0*y[8]- 8.0*y[9];"
 255 "  res[9]=yd[9]                                                                               -10.0*y[8]+10.0*y[9];"
 256 "}"
 257 "void myjac(double *t,double *y,double *yd,double *res,double *cj,double *rpar,int *ipar)"
 258 "{"
 259 "  res[0]=0.0;"
 260 "  res[1]=0.0;"
 261 "  res[2]=0.0;"
 262 "  res[3]=-17.0+*cj;"
 263 "  res[4]=8.0;"
 264 "  res[5]=0.0;"
 265 "  res[6]=0.0;"
 266 "  res[7]=6.0;"
 267 "  res[8]=-12.0+*cj;"
 268 "  res[9]=7.0;"
 269 "  res[10]=0.0;"
 270 "  res[11]=3.0;"
 271 "  res[12]=9.0;"
 272 "  res[13]=-17.0+*cj;"
 273 "  res[14]=3.0;"
 274 "  res[15]=0.0;"
 275 "  res[16]=4.0;"
 276 "  res[17]=3.0;"
 277 "  res[18]=-13.0+*cj;"
 278 "  res[19]=4.0;"
 279 "  res[20]=0.0;"
 280 "  res[21]=8.0;"
 281 "  res[22]=2.0;"
 282 "  res[23]=-18.0+*cj;"
 283 "  res[24]=7.0;"
 284 "  res[25]=0.0;"
 285 "  res[26]=2.0;"
 286 "  res[27]=6.0;"
 287 "  res[28]=-13.0+*cj;"
 288 "  res[29]=7.0;"
 289 "  res[30]=0.0;"
 290 "  res[31]=4.0;"
 291 "  res[32]=7.0;"
 292 "  res[33]=-16.0+*cj;"
 293 "  res[34]=5.0;"
 294 "  res[35]=0.0;"
 295 "  res[36]=0.0;"
 296 "  res[37]=7.0;"
 297 "  res[38]=-17.0+*cj;"
 298 "  res[39]=4.0;"
 299 "  res[40]=0.0;"
 300 "  res[41]=9.0;"
 301 "  res[42]=8.0;"
 302 "  res[43]=-14.0+*cj;"
 303 "  res[44]=10.0;"
 304 "  res[45]=0.0;"
 305 "  res[46]=1.0;"
 306 "  res[47]=8.0;"
 307 "  res[48]=-10.0+*cj;"
 308 "  res[49]=0.0;"
 309 "}"];
 310 mputl(ccode,"band.c"); //create the C file of myjac
 311 ilib_for_link(["myres","myjac"],"band.c","","c");//compile
 312 exec("loader.sce");
 313 y0=ones(n,1);yd0=0*y0;
 314 yb=dae([y0,yd0],0,0:0.1:10,"myres","myjac");
 315 a = norm(y-yb);
 316 if (a > %eps * 1e5) then pause,end
 317
 318 //                 Root finder (daskr)
 319 //
 320 y0=1;t=2:6;t0=1;y0d=3;
 321 %DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
 322 atol=1.d-6;rtol=0;ng=2;
 323 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1","psol1","pjac1");
 324 assert_checkalmostequal(nn(1),2.47,0.001);
 325 y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
 326 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1","psol1","pjac1");
 327 assert_checkalmostequal(nn(1),2.5,0.001);
 328 y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
 329 %DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 330 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 331 assert_checkalmostequal(nn(1),2.53,0.001);
 332
 333 // Same problem, but using macro for the derivative evaluation function 'res1'
 334 deff("[delta,ires]=res1(t,y,ydot)","ires=0;delta=ydot-((2.*log(y)+8)./t-5).*y")
 335 deff("[rts]=gr1(t,y,yd)","rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]")
 336 y0=1;t=2:6;t0=1;y0d=3;
 337 atol=1.d-6;rtol=0;ng=2;
 338 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
 339 assert_checkalmostequal(nn(1),2.47,0.001);
 340 y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
 341 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
 342 assert_checkalmostequal(nn(1),2.5,0.001);
 343 y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
 344 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
 345 assert_checkalmostequal(nn(1),2.53,0.001);
 346
 347 // Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
 348 // pjac uses the macro res1 defined above.
 349 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
 350     ires = 0;
 351     SQuround = 1.490D-08;
 352     tx = t;
 353     nrow = 0;
 354     e = zeros(1, neq);
 355     wp = zeros(neq*neq, 1);
 356     iwp = zeros(neq*neq, 2);
 357     for i=1:neq
 358         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
 359         if h*ydot(i) < 0 then del = -del; end
 360         ysave = y(i);
 361         ypsave = ydot(i);
 362         y(i) = y(i) + del;
 363         ydot(i) = ydot(i) + cj*del;
 364         [e ires] = res1(tx, y, ydot);
 365         if ires < 0 then
 366             ires = -1;
 367             return;
 368         end
 369         delinv = 1/del;
 370         for j=1:neq
 371             wp(nrow+j) = delinv*(e(j)-savr(j));
 372             if isnan(wp(nrow+j)) then
 373                 ires = -1;
 374                 return;
 375             end
 376             iwp(nrow+j, 1) = i;
 377             iwp(nrow+j, 2) = j;
 378         end
 379         nrow = nrow + neq;
 380         y(i) = ysave;
 381         ydot(i) = ypsave;
 382     end
 383 endfunction
 384 function [r, ier] = psol(wp, iwp, b)
 385     ier = 0;
 386     //Compute the LU factorization of R.
 387     sp = sparse(iwp, wp);
 388     [h, rk] = lufact(sp);
 389     //Solve the system LU*X = b
 390     r = lusolve(h, b);
 391     ludel(h);
 392 endfunction
 393 y0=1;t=2:6;t0=1;y0d=3;
 394 %DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
 395 atol=1.d-6;rtol=0;ng=2;
 396 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,"gr1",psol,pjac);
 397 assert_checkalmostequal(nn(1),2.47,0.001);
 398 y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
 399 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,"gr1",psol,pjac);
 400 assert_checkalmostequal(nn(1),2.5,0.001);
 401 y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
 402 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,"gr1",psol,pjac);
 403 assert_checkalmostequal(nn(1),2.53,0.001);
 404 //C
 405 //C-----------------------------------------------------------------------
 406 //C Second problem (Van Der Pol oscillator).
 407 //C The initial value problem is..
 408 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
 409 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
 410 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
 411 //C The root function is  G = Y1.
 412 //C An analytic solution is not known, but the zeros of Y1 are known
 413 //C to 15 figures for purposes of checking the accuracy.
 414 //C-----------------------------------------------------------------------
 415 %DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 416 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
 417 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
 418 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
 419 assert_checkalmostequal(nn(1),81.163512,0.001);
 420
 421 deff("[delta,ires]=res2(t,y,ydot)",...
 422 "ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]")
 423 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,"jac2",ng,"gr2");
 424 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)]")
 425 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2");
 426 deff("s=gr2(t,y,yd)","s=y(1)")
 427 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2);
 428
 429 // Same problem, with psol and pjac example routines
 430 %DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
 431 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2","psol1","pjac1");
 432 assert_checkalmostequal(nn(1),81.163512,0.009);
 433 deff("s=gr2(t,y,yd)","s=y(1)")
 434 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2,"psol1","pjac1");
 435 assert_checkalmostequal(nn(1),81.163512,0.009);
 436
 437 // Same problem, with psol and pjac macros
 438 // Redefine pjac to use res2
 439 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
 440     ires = 0;
 441     SQuround = 1.490D-08;
 442     tx = t;
 443     nrow = 0;
 444     e = zeros(1, neq);
 445     wp = zeros(neq*neq, 1);
 446     iwp = zeros(neq*neq, 2);
 447     for i=1:neq
 448         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
 449         if h*ydot(i) < 0 then del = -del; end
 450         ysave = y(i);
 451         ypsave = ydot(i);
 452         y(i) = y(i) + del;
 453         ydot(i) = ydot(i) + cj*del;
 454         [e ires]=res2(tx, y, ydot);
 455         if ires < 0 then return; end
 456         delinv = 1/del;
 457         for j=1:neq
 458             wp(nrow+j) = delinv*(e(j)-savr(j));
 459             iwp(nrow+j,1) = i;
 460             iwp(nrow+j,2) = j;
 461         end
 462         nrow = nrow + neq;
 463         y(i) = ysave;
 464         ydot(i) = ypsave;
 465     end
 466 endfunction
 467 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2",psol,pjac);
 468 assert_checkalmostequal(nn(1),81.163512,0.003);
 469 deff("s=gr2(t,y,yd)","s=y(1)")
 470 [yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2,psol,pjac);
 471 assert_checkalmostequal(nn(1),81.163512,0.003);
 472 %DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 473
 474 //           Hot Restart
 475 [yy,nn,hotd]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
 476 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
 477 [yy,nn,hotd]=dae("root2",[y01,y0d1],t01,t,rtol,atol,"res2","jac2",ng,"gr2",hotd);
 478 assert_checkalmostequal(nn(1),162.57763,0.004);