a89d494baa0ede761ce7e933dc02b3969ab5e181
[scilab.git] / 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);