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