update compilation path with visual studio 2012 with Win8 sdk path and update some...
[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / daskr.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING,
7 // which you should have received as part of this distribution.
8 // The terms are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
10 // =============================================================================
11
12 //<-- ENGLISH IMPOSED -->
13
14 //C-----------------------------------------------------------------------
15 //C First problem.
16 //C The initial value problem is..
17 //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
18 //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
19 //C The two root functions are..
20 //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
21 //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
22 //C-----------------------------------------------------------------------
23 y0=1;t=2:6;t0=1;y0d=3;
24 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
25 atol=1.d-6;rtol=0;ng=2;
26 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,"res1",ng,"gr1",info,"psol1","pjac1");
27 assert_checkalmostequal(nn(1),2.47,0.001);
28 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
29 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,"res1",ng,"gr1",info,"psol1","pjac1");
30 assert_checkalmostequal(nn(1),2.5,0.001);
31 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
32 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
33 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,"res1",ng,"gr1",info);
34 assert_checkalmostequal(nn(1),2.500009,0.001);
35
36 // Same problem, but using macro for the derivative evaluation function 'res1'
37 deff("[delta,ires]=res1(t,y,ydot)","ires=0;delta=ydot-((2.*log(y)+8)./t-5).*y")
38 deff("[rts]=gr1(t,y,yd)","rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]")
39
40 y0=1;t=2:6;t0=1;y0d=3;
41 atol=1.d-6;rtol=0;ng=2;
42 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
43 assert_checkalmostequal(nn(1),2.47,0.001);
44 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
45 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
46 assert_checkalmostequal(nn(1),2.5,0.001);
47 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
48 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
49 assert_checkalmostequal(nn(1),2.53,0.001);
50
51 // Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
52 // pjac uses the macro res1 defined above.
53 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
54     ires = 0;
55     SQuround = 1.490D-08;
56     tx = t;
57     nrow = 0;
58     e = zeros(1, neq);
59     wp = zeros(neq*neq, 1);
60     iwp = zeros(neq*neq, 2);
61     for i=1:neq
62         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
63         if h*ydot(i) < 0 then del = -del; end
64         ysave = y(i);
65         ypsave = ydot(i);
66         y(i) = y(i) + del;
67         ydot(i) = ydot(i) + cj*del;
68         [e ires] = res1(tx, y, ydot);
69         if ires < 0 then
70             ires = -1;
71             return;
72         end
73         delinv = 1/del;
74         for j=1:neq
75             wp(nrow+j) = delinv*(e(j)-savr(j));
76             if isnan(wp(nrow+j)) then
77                 ires = -1;
78                 return;
79             end
80             iwp(nrow+j, 1) = i;
81             iwp(nrow+j, 2) = j;
82         end
83         nrow = nrow + neq;
84         y(i) = ysave;
85         ydot(i) = ypsave;
86     end
87 endfunction
88 function [r, ier] = psol(wp, iwp, b)
89     ier = 0;
90     //Compute the LU factorization of R.
91     sp = sparse(iwp, wp);
92     [h, rk] = lufact(sp);
93     //Solve the system LU*X = b
94     r = lusolve(h, b);
95     ludel(h);
96 endfunction
97
98 y0=1;t=2:6;t0=1;y0d=3;
99 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
100 atol=1.d-6;rtol=0;ng=2;
101 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,"gr1",info,psol,pjac);
102 assert_checkalmostequal(nn(1),2.47,0.001);
103 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
104 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,"gr1",info,psol,pjac);
105 assert_checkalmostequal(nn(1),2.5,0.001);
106 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
107 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,"gr1",info,psol,pjac);
108 assert_checkalmostequal(nn(1),2.53,0.001);
109
110 //C
111 //C-----------------------------------------------------------------------
112 //C Second problem (Van Der Pol oscillator).
113 //C The initial value problem is..
114 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
115 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
116 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
117 //C The root function is  G = Y1.
118 //C An analytic solution is not known, but the zeros of Y1 are known
119 //C to 15 figures for purposes of checking the accuracy.
120 //C-----------------------------------------------------------------------
121 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
122 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
123 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
124 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,"res2","jac2",ng,"gr2",info);
125 assert_checkalmostequal(nn(1),81.163512,0.001);
126
127 deff("[delta,ires]=res2(t,y,ydot)",...
128 "ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]")
129 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,"jac2",ng,"gr2",info);
130 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)]")
131 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,"gr2",info);
132 deff("s=gr2(t,y,yd)","s=y(1)")
133 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
134
135 // Same problem, with psol and pjac example routines
136
137 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
138 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,"gr2",info,"psol1","pjac1");
139 assert_checkalmostequal(nn(1),81.163512,0.009);
140 deff("s=gr2(t,y,yd)","s=y(1)")
141 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,"psol1","pjac1");
142 assert_checkalmostequal(nn(1),81.163512,0.009);
143
144 // Same problem, with psol and pjac macros
145
146 // Redefine pjac to use res2
147 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
148     ires = 0;
149     SQuround = 1.490D-08;
150     tx = t;
151     nrow = 0;
152     e = zeros(1, neq);
153     wp = zeros(neq*neq, 1);
154     iwp = zeros(neq*neq, 2);
155     for i=1:neq
156         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
157         if h*ydot(i) < 0 then del = -del; end
158         ysave = y(i);
159         ypsave = ydot(i);
160         y(i) = y(i) + del;
161         ydot(i) = ydot(i) + cj*del;
162         [e ires]=res2(tx, y, ydot);
163         if ires < 0 then return; end
164         delinv = 1/del;
165         for j=1:neq
166             wp(nrow+j) = delinv*(e(j)-savr(j));
167             iwp(nrow+j,1) = i;
168             iwp(nrow+j,2) = j;
169         end
170         nrow = nrow + neq;
171         y(i) = ysave;
172         ydot(i) = ypsave;
173     end
174 endfunction
175 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,"gr2",info,psol,pjac);
176 assert_checkalmostequal(nn(1),81.163512,0.003);
177 deff("s=gr2(t,y,yd)","s=y(1)")
178 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,psol,pjac);
179 assert_checkalmostequal(nn(1),81.163512,0.003);
180 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
181
182 //           Hot Restart
183
184 [yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,"res2","jac2",ng,"gr2",info);
185 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(3:4,qq);y0d1=yy(4:5,qq);
186 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,"res2","jac2",ng,"gr2",info,hotd);
187 assert_checkalmostequal(nn(1),162.57763,0.004);
188
189 // Same with C code
190 ilib_verbose(0);
191
192 cd TMPDIR;
193 mkdir("daskr_test1");
194 cd("daskr_test1");
195
196 code=["#include <math.h>"
197 "void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)"
198 "{res[0] = yd[0] - y[1];"
199 " res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}"
200 " "
201 "void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)"
202 "{pd[0]=*cj - 0.0;"
203 " pd[1]=    - (-200.0*y[0]*y[1] - 1.0);"
204 " pd[2]=    - 1.0;"
205 " pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}"
206 " "
207 "void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)"
208 "{ groot[0] = y[0];}"];
209 mputl(code,"t22.c") ;
210 ilib_for_link(["res22" "jac22" "gr22"],"t22.c","","c");
211 exec("loader.sce");
212
213 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
214 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
215 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
216 // Hot restart
217 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(3:4,qq);y0d1=yy(4:5,qq);
218 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,"res22","jac22",ng,"gr22",info,hotd);
219
220 rtol=[1.d-6;1.d-6];
221 atol=[1.d-6;1.d-4];
222 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
223 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,"res22","jac22",ng,"gr22",info);
224 // Hot restart
225 [yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,"res22","jac22",ng,"gr22",info);
226 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(3:4,qq);y0d1=yy(4:5,qq);
227 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,"res22","jac22",ng,"gr22",info,hotd);