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