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