d21c4427cb6319029d71bf70b1289e1866a1a63c
[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / ode.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2007-2008 - INRIA
4 // Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
5 // Copyright (C) 2013 - Scilab Enterprises - Adeline CARNIS
6 //
7 //  This file is distributed under the same license as the Scilab package.
8 // =============================================================================
9
10 // <-- CLI SHELL MODE -->
11 warning("off")
12 ilib_verbose(0);
13
14 version = getversion("scilab");
15
16 // to check that ode works
17 // ---------- Simple one dimension ODE (Scilab function external)
18 // dy/dt=y^2-y sin(t)+cos(t), y(0)=0
19 function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
20 y0=0;t0=0;t=0:0.1:%pi;
21 y=ode(y0,t0,t,f);
22 assert_checkalmostequal(size(y), [1 32] , %eps, [], "matrix");
23 clear y;
24 clear t;
25 //*************************** function F and lsoda ********************************/
26 // create functions
27 cd TMPDIR;
28
29 CC=["void fex1(int* neq, double* t, double* y, double* ydot)"
30 "{"
31 "   ydot[0] = -0.04*y[0] + 1.0e+4*y[1]*y[2];"
32 "   ydot[2] = 3.0e+7*y[1]*y[1];"
33 "   ydot[1] = -ydot[0] - ydot[2];"
34 "}"];
35 mputl(CC,TMPDIR+"/fex1.c");
36 ilib_for_link("fex1","fex1.c",[],"c");
37 exec loader.sce;
38
39 C=[ "void fex2(int* neq, double* t, double* y, double* ydot)"
40 "{"
41 "   ydot[0] = y[4]*y[0] + y[5]*y[1]*y[2];"
42 "   ydot[2] = y[3]*y[1]*y[1];"
43 "   ydot[1] = -ydot[0] - ydot[2];"
44 "}"];
45
46 mputl(C,TMPDIR+"/fex2.c");
47 ilib_for_link("fex2","fex2.c",[],"c");
48 exec loader.sce;
49 clear f;
50 function ydot = f(t,yin)
51     ydot(1)=-0.040*yin(1) + 1.0D4*yin(2)*yin(3);
52     ydot(3)=3.0D7*yin(2)**2;
53     ydot(2)=-ydot(1) - ydot(3);
54 endfunction
55
56 function ydot = f1(t,yin,a,b,c)
57     ydot(1)=b*yin(1) + c*yin(2)*yin(3);
58     ydot(3)=a*yin(2)**2;
59     ydot(2)=-ydot(1) - ydot(3);
60 endfunction
61
62 function ydot = f2(t,yin,a)
63     ydot(1)=a(2)*yin(1) + a(3)*yin(2)*yin(3);
64     ydot(3)=a(1)*yin(2)**2;
65     ydot(2)=-ydot(1) - ydot(3);
66 endfunction
67
68 // init variables
69 y(1)    = 1;
70 y(2)    = 0;
71 y(3)    = 0;
72 t       = 0;
73 tout    = 0.4*exp((0:11)*log(10));
74 rtol    = 1.0d-4;
75 atol(1) = 1.0d-6;
76 atol(2) = 1.0d-10;
77 atol(3) = 1.0d-6;
78
79 // result provide by lsoda documentation.
80 // on a cdc-7600 in single precision.
81 //   at t =  4.0000e-01
82 resDoc(:,1) = [ 9.851712e-01 ; 3.386380e-05 ; 1.479493e-02 ];
83 //   at t =  4.0000e+00
84 resDoc(:,2) = [ 9.055333e-01 ; 2.240655e-05 ; 9.444430e-02 ];
85 //   at t =  4.0000e+01
86 resDoc(:,3) = [ 7.158403e-01 ; 9.186334e-06 ; 2.841505e-01 ];
87 //   at t =  4.0000e+02
88 resDoc(:,4) = [ 4.505250e-01 ; 3.222964e-06 ; 5.494717e-01 ];
89 //   at t =  4.0000e+03
90 resDoc(:,5) = [ 1.831975e-01 ; 8.941774e-07 ; 8.168016e-01 ];
91 //   at t =  4.0000e+04
92 resDoc(:,6) = [ 3.898730e-02 ; 1.621940e-07 ; 9.610125e-01 ];
93 //   at t =  4.0000e+05
94 resDoc(:,7) = [ 4.936363e-03 ; 1.984221e-08 ; 9.950636e-01 ];
95 //   at t =  4.0000e+06
96 resDoc(:,8) = [ 5.161831e-04 ; 2.065786e-09 ; 9.994838e-01 ];
97 //   at t =  4.0000e+07
98 resDoc(:,9) = [ 5.179817e-05 ; 2.072032e-10 ; 9.999482e-01 ];
99 //   at t =  4.0000e+08
100 resDoc(:,10) = [ 5.283401e-06 ; 2.113371e-11 ; 9.999947e-01 ];
101 //   at t =  4.0000e+09
102 resDoc(:,11) = [ 4.659031e-07 ; 1.863613e-12 ; 9.999995e-01 ];
103 //   at t =  4.0000e+10
104 resDoc(:,12) = [ 1.404280e-08 ; 5.617126e-14 ; 1.000000e+00 ];
105
106 // f as a string (dynamic link function)
107 res  = ode(y, t,tout, rtol, atol, "fex1");
108 // f as a list(string,...
109 if version(1) > 5 then
110     res2 = ode(y, t, tout, rtol, atol, list("fex2", 3.0d+7, -0.04, 1.0d+4));
111 end
112 // f as a macro
113 res3 = ode(y, t, tout, rtol, atol, f);
114 // f as a list(macro,...
115 res4 = ode(y, t, tout, rtol, atol, list(f1, 3.0d+7, -0.04, 1.0d+4));
116 args = [ 3.0d+7, -0.04, 1.0d+4 ];
117 res5 = ode(y, t, tout, rtol, atol, list(f2, args));
118 // f as a string (static link function)
119 res6 = ode(y, t,tout, rtol, atol, "fex");
120
121 // check results
122 assert_checkalmostequal(resDoc, res, 2d-7, [], "matrix"); // There are a little diff between resDoc and res
123 if version(1) > 5 then
124     assert_checkalmostequal(res, res2, 2d-7, [], "matrix"); // because results provides by lsoda
125 end
126 assert_checkalmostequal(res, res3, 2d-7, [], "matrix"); // documentation are in single precision.
127 assert_checkalmostequal(res, res4, 2d-7, [], "matrix");
128 assert_checkalmostequal(res, res5, 2d-7, [], "matrix");
129 assert_checkalmostequal(res, res6, 2d-7, [], "matrix");
130
131 //*************************** w iw ********************************/
132 tout2 = 0.4*exp(12*log(10));
133 tout3 = 0.4*exp((0:12)*log(10));
134 [yout w iw] = ode(y, t, tout, rtol, atol, f);
135 yout1 = ode(y, t, tout2, rtol, atol, f, w, iw);
136 yout2 = ode(y, t, tout3, rtol, atol, f);
137
138 assert_checkalmostequal(yout2(3*12+1:3*13), yout1, %eps, [], "matrix");
139
140 //*************************** Polynom ********************************/
141 //y(1) = 1;
142 //y(2) = 2;
143 //y(3) = 3;
144 //yy   = 1+2*%s+3*%s*%s;
145
146 //res  = ode(y, t,tout, rtol, atol, 'fex1');
147
148 //res6  = ode(yy, t, tout, rtol, atol, 'fex1');
149 //for i=1:12, assert_checkalmostequal(poly(res(:,i), "s", "coeff"),res6(i), %eps); end
150 //res7 = ode(yy, t, tout, rtol, atol, list('fex2', 3.0d+7, -0.04, 1.0d+4));
151 //for i=1:12, assert_checkequal(poly(res(:,i), "s", "coeff"), res7(i), %eps); end
152
153 //*************************** function Jac and lsode ********************************/
154 CC=["void jac(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)"
155 "{"
156 "  j[0] = y[6]; "
157 "  j[1] = y[7]; "
158 "  j[2] = y[8]; "
159 "  j[3] = y[9]; "
160 "}"];
161
162 mputl(CC,TMPDIR+"/jac.c");
163 ilib_for_link("jac","jac.c",[],"c");
164 exec loader.sce;
165
166 CC=["void jac2(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)"
167 "{"
168 "  j[0] = 10; "
169 "  j[1] = 0; "
170 "  j[2] = 0; "
171 "  j[3] = -1; "
172 "}"];
173 mputl(CC,TMPDIR+"/jac2.c");
174 ilib_for_link("jac2","jac2.c",[],"c");
175 exec loader.sce;
176
177 // ydot = A * y
178 C=[ "void fext(int* neq, double* t, double* y, double* ydot)"
179 "{"
180 "   ydot[0] = y[2]*y[0] + y[4]*y[1];"
181 "   ydot[1] = y[3]*y[0] + y[5]*y[1];"
182 "}"];
183
184 mputl(C,TMPDIR+"/fext.c");
185 ilib_for_link("fext","fext.c",[],"c");
186 exec loader.sce;
187
188 C=[ "void fext2(int* neq, double* t, double* y, double* ydot)"
189 "{"
190 "   ydot[0] = 10*y[0] + 0*y[1];"
191 "   ydot[1] = 0*y[0] + (-1)*y[1];"
192 "}"];
193
194 mputl(C,TMPDIR+"/fext2.c");
195 ilib_for_link("fext2","fext2.c",[],"c");
196 exec loader.sce;
197 clear f;
198
199 function ydot=f(t, y)
200     ydot=A*y
201 endfunction
202
203 function J=Jacobian(t, y)
204     J=A
205 endfunction
206
207 A=[10,0;0,-1];
208 y0=[0;1];
209 t0=0;
210 t=1;
211
212 res  = ode("stiff", y0, t0, t, f, Jacobian);
213 if version(1) > 5 then
214     res1 = ode("stiff", y0, t0, t, list("fext", 10, 0, 0,-1), list("jac", A));
215 end
216 res2 = ode("stiff", y0, t0, t, "fext2", "jac2");
217 res3 = ode("stiff", y0, t0, t, f, "jac2");
218 res4 = ode("stiff", y0, t0, t, "fext2", Jacobian);
219
220 assert_checkalmostequal(res, expm(A*t)*y0, 1.0D-7, [], "matrix");
221 if version(1) > 5 then
222     assert_checkalmostequal(res, res1, %eps, [], "matrix");
223 end
224 assert_checkalmostequal(res, res2, %eps, [], "matrix");
225 assert_checkalmostequal(res, res3, %eps, [], "matrix");
226 assert_checkalmostequal(res, res4, %eps, [], "matrix");
227
228 //*************************** discrete ********************************/
229 function yp=a_function(k,y)
230     yp=A*y+B*u(k);
231 endfunction
232
233 y1 = [1;2;3];
234 A  = diag([0.2,0.5,0.9]);
235 B  = [1;1;1];
236 u  = 1:10;
237 n  = 5;
238
239 y =ode("discrete", y1, 1, 1:n, a_function);
240 for i = 1:4, y1(:,i+1) = A * y1(:,i) + B * u(i); end
241 assert_checkalmostequal(y, y1, %eps, [], "matrix");
242
243 // Now y evaluates  at [y3,y5,y7,y9]
244 y1 = [1;2;3];
245 t  = 3:2:9;
246 y  = ode("discrete", y1, 1, t, a_function);
247 for i=1:9, y1(:,i+1) = A * y1(:,i) + B * u(i); end
248 y1 = y1(:,t);
249 assert_checkalmostequal(y, y1, %eps, [], "matrix");
250
251 //*************************** root ********************************/
252 y0=1;
253 ng=1;
254
255 C=[ "void fextern(int* neq, double* t, double* y, double* ydot)"
256 "{"
257 "int i = 0;"
258 "for( i = 0; i < *neq; i++)"
259 "   ydot[i] = y[i];"
260 "}"];
261
262 mputl(C,TMPDIR+"/fextern.c");
263 ilib_for_link("fextern","fextern.c",[],"c");
264 exec loader.sce;
265 clear f;
266 function ydot=f(t,y)
267     ydot=y;
268 endfunction
269
270 clear g;
271 // check rd result
272 function z=g(t,y)
273     z=y-2;
274 endfunction
275 [y,rd]=ode("root",y0,0,2,f,ng,g);
276 assert_checkequal( rd(2) <> 1 , %f);
277
278 clear g;
279 function z=g(t,y)
280     z=y-[2;2;33];
281 endfunction
282 [y,rd]=ode("root",y0,0,2,f,3,g);
283 assert_checkequal( rd(2) <> 1 , %f);
284 assert_checkequal( rd(3) <> 2 , %f);
285
286 clear g;
287 function z=g(t,y)
288     z=y-[2;2;2];
289 endfunction
290 [y,rd]=ode("root",y0,0,2,f,3,g);
291 assert_checkequal( rd(2) <> 1 , %f);
292 assert_checkequal( rd(3) <> 2 , %f);
293 assert_checkequal( rd(4) <> 3 , %f);
294
295 clear g;
296 function z=g(t,y)
297     z=y-[2;6;2;2];
298 endfunction
299 [y,rd]=ode("root",y0,0,2,f,4,g);
300 assert_checkequal( rd(2) <> 1 , %f);
301 assert_checkequal( rd(3) <> 3 , %f);
302 assert_checkequal( rd(4) <> 4 , %f);
303
304 // check y result
305
306 // result provide by lsodar documentation.
307 // on a cdc-7600 in single precision.
308 //   at t =  2.6400e-01
309 resDocRoot(:,1) = [ 9.899653e-01 ; 3.470563e-05 ; 1.000000e-02 ];
310 //        the above line is a root,  jroot =    0    1
311 //   at t =  4.0000e-01
312 resDoc(:,1) = [ 9.851712e-01 ; 3.386380e-05 ; 1.479493e-02 ];
313 //   at t =  4.0000e+00
314 resDoc(:,2) = [ 9.055333e-01 ; 2.240655e-05 ; 9.444430e-02 ];
315 //   at t =  4.0000e+01
316 resDoc(:,3) = [ 7.158403e-01 ; 9.186334e-06 ; 2.841505e-01 ];
317 //   at t =  4.0000e+02
318 resDoc(:,4) = [ 4.505250e-01 ; 3.222964e-06 ; 5.494717e-01 ];
319 //   at t =  4.0000e+03
320 resDoc(:,5) = [ 1.831975e-01 ; 8.941774e-07 ; 8.168016e-01 ];
321 //   at t =  4.0000e+04
322 resDoc(:,6) = [ 3.898730e-02 ; 1.621940e-07 ; 9.610125e-01 ];
323 //   at t =  4.0000e+05
324 resDoc(:,7) = [ 4.936363e-03 ; 1.984221e-08 ; 9.950636e-01 ];
325 //   at t =  4.0000e+06
326 resDoc(:,8) = [ 5.161831e-04 ; 2.065786e-09 ; 9.994838e-01 ];
327 //   at t =  2.0745e+07
328 resDocRoot(:,2) = [ 1.000000e-04 ; 4.000395e-10 ; 9.999000e-01 ];
329 //        the above line is a root,  jroot =    1    0
330 //   at t =  4.0000e+07
331 resDoc(:,9) = [ 5.179817e-05 ; 2.072032e-10 ; 9.999482e-01 ];
332 //   at t =  4.0000e+08
333 resDoc(:,10) = [ 5.283401e-06 ; 2.113371e-11 ; 9.999947e-01 ];
334 //   at t =  4.0000e+09
335 resDoc(:,11) = [ 4.659031e-07 ; 1.863613e-12 ; 9.999995e-01 ];
336 //   at t =  4.0000e+10
337 resDoc(:,12) = [ 1.404280e-08 ; 5.617126e-14 ; 1.000000e+00 ];
338
339 G=[ "void gex(int* neq, double* t, double* y, int* ng, double* gout)"
340 "{"
341 "gout[0] = y[0] - 1.0e-4;"
342 "gout[1] = y[2] - 1.0e-2;"
343 "}"];
344
345 mputl(G,TMPDIR+"/gex.c");
346 ilib_for_link("gex","gex.c",[],"c");
347 exec loader.sce;
348
349 y(1) = 1;
350 y(2) = 0;
351 y(3) = 0;
352 t0   = 0;
353
354 [yout1,rd1,w,iw] = ode("root", y, t0, tout, "fex1", 2, "gex");
355 assert_checkalmostequal(rd1(1), 2.64d-01, 1d-4);
356 [yout2,rd2,w,iw] = ode("root", y, t0, tout, "fex1", 2, "gex", w, iw);
357 assert_checkalmostequal(rd2(1), 2.0795776d+07, 4d-5);
358 err = execstr("[yout3,rd,w,iw] = ode(""root"", y, t0, tout, ""fex1"", 2, ""gex"", w, iw);","errcatch");
359 assert_checkequal( err == 0 , %f);
360
361 // check results
362 assert_checkalmostequal(resDocRoot(:,1), yout1, 2.0D-8, [], "matrix");
363 assert_checkalmostequal(resDocRoot(:,2), yout2(:,9), 2.0D-8, [], "matrix");
364 assert_checkalmostequal(resDoc(:,1:8), yout2(:,1:8), 2.0D-4, [], "matrix");
365
366 //*************************** rk/rkf/fix ********************************/
367 function ydot=functionF(t, y)
368     ydot=y^2-y*sin(t)+cos(t)
369 endfunction
370
371 y0 = 0;
372 t0 = 0;
373 t  = 0:0.1:%pi;
374
375 rk   = ode("rk",  y0, t0, t, functionF);
376 rkf  = ode("rkf", y0, t0, t, functionF);
377 fixx = ode("fix", y0, t0, t, functionF);
378
379 rkRes = [0.    0.09983341664683527    0.19866933079512300    0.29552020666153711 0.38941834230905709    0.47942553860493514    0.56464247339623352    0.64421768723947215 0.71735609090195429    0.78332690963060869    0.8414709848117728     0.89120736006615475 0.93203908597292051    0.96355818542403104    0.98544972999664104    0.99749498661376];
380 rkRes(17:32) = [0.99957360305287668    0.99166481046564392    0.97384763089330029    0.94630008770455387 0.90929742684492987    0.86320936667026738    0.80849640384312127    0.74570521220233155 0.67546318057873300    0.59847214413334937    0.51550137185246248    0.42737988026620155 0.33498815018939587    0.23924932924832210    0.14112000809477554    0.04158066246848785];
381
382 rkfRes = [0.    0.09983341667063099    0.19866933087782307    0.2955202068043328    0.38941834246046680     0.47942553864900261    0.56464247314940919    0.64421768645637856    0.71735608928892303      0.78332690686480611    0.84147098056308622    0.89120735401875306    0.93203907784361117      0.96355817497519225    0.98544971704249074    0.99749497101988072    0.99957358472991464];
383 rkfRes(18:32) = [    0.99166478935902302    0.97384760697150774    0.94630006094857066    0.90929739724116376     0.86320933420871493    0.80849636852157181    0.74570517403636893    0.67546313961625104      0.59847210047141619    0.51550132565379336    0.42737983177229999    0.33498809972769594      0.23924927723128114    0.14111995500988161    0.04158060885936282];
384
385 assert_checkalmostequal(rkRes, rk, %eps * 20, [], "matrix");
386 assert_checkalmostequal(rkfRes, rkf, %eps, [], "matrix");
387 assert_checkalmostequal(rkf, fixx, %eps, [], "matrix");
388 warning("on")