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