dbbe1cdd1bb1c40066970f4467c4a9594cce3fbe
[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / daskr.tst
1
2 // =============================================================================
3 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 // Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
5 //
6 // This file must be used under the terms of the CeCILL.
7 // This source file is licensed as described in the file COPYING,
8 // which you should have received as part of this distribution.
9 // The terms are also available at
10 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 // =============================================================================
12
13 // Execute with exec("SCI/modules/differential_equations/tests/unit_tests/daskr.tst");
14 //  or test_run('differential_equations', 'daskr', ['no_check_error_output']);
15
16 //C-----------------------------------------------------------------------
17 //C First problem.
18 //C The initial value problem is..
19 //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
20 //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
21 //C The two root functions are..
22 //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
23 //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
24 //C-----------------------------------------------------------------------
25 y0=1;t=2:6;t0=1;y0d=3;
26 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
27 atol=1.d-6;rtol=0;ng=2;
28 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
29 if abs(nn(1)-2.47)>0.001 then pause,end
30 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
31 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
32 if abs(nn(1)-2.5)>0.001 then pause,end
33 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
34 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
35 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
36 if abs(nn(1)-2.53)>0.001 then pause,end
37
38 deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
39 deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
40
41 y0=1;t=2:6;t0=1;y0d=3;
42 atol=1.d-6;rtol=0;ng=2;
43 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
44 if abs(nn(1)-2.47)>0.001 then pause,end
45 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
46 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
47 if abs(nn(1)-2.5)>0.001 then pause,end
48 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
49 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
50 if abs(nn(1)-2.53)>0.001 then pause,end
51
52 //C
53 //C-----------------------------------------------------------------------
54 //C Second problem (Van Der Pol oscillator).
55 //C The initial value problem is..
56 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
57 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
58 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
59 //C The root function is  G = Y1.
60 //C An analytic solution is not known, but the zeros of Y1 are known
61 //C to 15 figures for purposes of checking the accuracy.
62 //C-----------------------------------------------------------------------
63 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
64 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
65 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
66 if abs(nn(1)-81.163512)>0.001 then pause,end
67
68 deff('[delta,ires]=res2(t,y,ydot)',...
69 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
70 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
71 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)]')
72 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
73 deff('s=gr2(t,y,yd)','s=y(1)')
74 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
75
76 //           Hot Restart
77
78 [yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
79 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
80 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hotd);
81 if abs(nn(1)-162.57763)>0.004 then pause,end
82
83 //Test of Dynamic link (Require f77!)
84 //         1 making the routines
85 res22=[...
86 '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
87 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
88 '      INTEGER NEQ';
89 '      DIMENSION Y(*), YDOT(*), DELTA(*)';
90 '      NEQ=2';
91 '      CALL F2(NEQ,T,Y,DELTA)';
92 '      DO 10 I = 1,NEQ';
93 '         DELTA(I) = YDOT(I) - DELTA(I)';
94 ' 10   CONTINUE';
95 '      RETURN';
96 '      END';
97 '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
98 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
99 '      INTEGER NEQ';
100 '      DOUBLE PRECISION T, Y, YDOT';
101 '      DIMENSION Y(*), YDOT(*)';
102 '      YDOT(1) = Y(2)';
103 '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
104 '      RETURN';
105 '      END';]
106
107 jac22=[...
108 '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
109 ' ';
110 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
111 '      INTEGER  NROWPD';
112 '      DOUBLE PRECISION T, Y, PD';
113 '      PARAMETER (NROWPD=2)';
114 '      DIMENSION Y(2), PD(NROWPD,2)';
115 '      PD(1,1) = 0.0D0';
116 '      PD(1,2) = 1.0D0';
117 '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
118 '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
119 '      PD(1,1) = CJ - PD(1,1)';
120 '      PD(1,2) =    - PD(1,2)';
121 '      PD(2,1) =    - PD(2,1)';
122 '      PD(2,2) = CJ - PD(2,2)';
123 '      RETURN';
124 '      END';]
125
126
127 gr22=[...
128 '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
129 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
130 '      INTEGER NEQ, NG';
131 '      DOUBLE PRECISION T, Y, GROOT';
132 '      DIMENSION Y(*), GROOT(*)';
133 '      GROOT(1) = Y(1)';
134 '      RETURN';
135 '      END';]
136
137 //Uncomment lines below: link may be machine dependent if some f77 libs are 
138 //needed for linking
139 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
140 //unix_g('cd /tmp;rm -f /tmp/jac22.f');
141 //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
142 //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
143 //unix_g('cd /tmp;make /tmp/jac22.o');
144 //          2  Linking the routines
145 //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
146 //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
147 //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
148 //          3 Calling the routines by dasrt
149 //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
150 // hot restart
151 //[yy,nn,hotd]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
152 //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
153 //[yy,nn,hotd]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hotd);