ODEPACK solvers (daskr): adding Scilab macro example for pjac and psol

[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / daskr.tst

diff --git a/scilab/modules/differential_equations/tests/unit_tests/daskr.tst b/scilab/modules/differential_equations/tests/unit_tests/daskr.tst
index dbbe1cd..686ce31 100644 (file)
--- a/scilab/modules/differential_equations/tests/unit_tests/daskr.tst
+++ b/scilab/modules/differential_equations/tests/unit_tests/daskr.tst
@@ -1,4 +1,3 @@
-
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
@@ -10,8 +9,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 // =============================================================================
 
-// Execute with exec("SCI/modules/differential_equations/tests/unit_tests/daskr.tst");
-//  or test_run('differential_equations', 'daskr', ['no_check_error_output']);
+// Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
 
 //C-----------------------------------------------------------------------
 //C First problem.
@@ -35,6 +33,7 @@ info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
 if abs(nn(1)-2.53)>0.001 then pause,end
 
+// Same problem, but using macro for the derivative evaluation function 'res1'
 deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
 deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
 
@@ -49,6 +48,59 @@ y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
 if abs(nn(1)-2.53)>0.001 then pause,end
 
+// Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
+// pjac uses the macro res1 defined above.
+function [r, ier] = psol(wp, iwp, b)
+    ier = 0;
+    //Compute the LU factorization of R.
+    sp = sparse(iwp, wp);
+    [h, rk] = lufact(sp);
+    //Solve the system LU*X = b
+    r = lusolve(h, b);
+    ludel(h);
+endfunction
+function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
+    ires = 0;
+    SQuround = 1.490D-08;
+    tx = t;
+    nrow = 0;
+    e = zeros(1, neq);
+    wp = zeros(neq*neq, 1);
+    iwp = zeros(neq*neq, 2);
+    for i=1:neq
+        del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
+        if h*ydot(i) < 0 then del = -del; end
+        ysave = y(i);
+        ypsave = ydot(i);
+        y(i) = y(i) + del;
+        ydot(i) = ydot(i) + cj*del;
+        [e ires]=res1(tx, y, ydot);
+        if ires < 0 then return; end
+        delinv = 1/del;
+        for j=1:neq
+            wp(nrow+j) = delinv*(e(j)-savr(j));
+            iwp(nrow+j,1) = i;
+            iwp(nrow+j,2) = j;
+        end
+        nrow = nrow + neq;
+        y(i) = ysave;
+        ydot(i) = ypsave;
+    end
+endfunction
+
+y0=1;t=2:6;t0=1;y0d=3;
+info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
+if abs(nn(1)-2.47)>0.001 then pause,end
+y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
+if abs(nn(1)-2.5)>0.001 then pause,end
+y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
+if abs(nn(1)-2.53)>0.001 then pause,end
+info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
+
 //C
 //C-----------------------------------------------------------------------
 //C Second problem (Van Der Pol oscillator).
@@ -134,7 +186,7 @@ gr22=[...
 '      RETURN';
 '      END';]
 
-//Uncomment lines below: link may be machine dependent if some f77 libs are 
+//Uncomment lines below: link may be machine dependent if some f77 libs are
 //needed for linking
 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
 //unix_g('cd /tmp;rm -f /tmp/jac22.f');