</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
- <synopsis>[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, jac]] [, hd])</synopsis>
+ <synopsis>[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, pjac]] [, hd])</synopsis>
</refsynopsisdiv>
<refsection>
<title>Arguments</title>
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
- <literal>[r,ier] = psol(neq,R,iR,b)</literal> and the
+ <literal>[r, ier] = psol(wp, iwp, b)</literal> and the
<literal>psol</literal> function must return the solution of that system in
<literal>r</literal> and an error flag <literal>ier</literal>.
</para>
<literal>psol</literal> is now
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-psol(R,iR,b,x1,x2,...)
+psol(wp,iwp,b,x1,x2,...)
]]></programlisting>
<para>
<literal>psol</literal> still returns the solution in <literal>r</literal>.
<listitem>
<para>A Scilab function.</para>
<para>Its calling sequence must be
- <literal>[R,iR,ier] = pjac(neq,t,y,ydot,h,cj)</literal> and in return,
- the arrays R and iR must contain all preconditioner information in compressed sparse row format.
+ <literal>[wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)</literal> and in return,
+ the arrays wp and iwp must contain all preconditioner information in compressed sparse row format.
</para>
</listitem>
<listitem>
<literal>pjac</literal> is now
</para>
<programlisting role="no-scilab-exec"><![CDATA[
-pjac(neq,t,y,ydot,h,cj,x1,x2,...)
+pjac(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)
]]></programlisting>
<para>
<literal>pjac</literal> still returns factorised
<literal>dg/dy + cj*dg/dydot</literal> as a function of
- <literal>(neq,t,y,ydot,h,cj,x1,x2,...)</literal>.
+ <literal>(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)</literal>.
</para>
</listitem>
<listitem>
double ysave = 0;
double ypsave = 0;
double* e = NULL;
-
+ int neq = *nequations;
double SQuround = sqrt(C2F(dlamch)("P"));
tx = *tOld;
- e = (double*)calloc(*nequations, sizeof(double));
+ e = (double*)calloc(neq, sizeof(double));
- for (i = 0 ; i < *nequations ; ++i)
+ for (i = 0 ; i < neq ; ++i)
{
del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
del *= (*h * actualP[i] >= 0) ? 1 : -1;
}
delinv = 1. / del;
- for (j = 0 ; j < *nequations ; j++)
+ for (j = 0 ; j < neq ; j++)
{
wp[nrow + j] = (e[j] - savr[j]) * delinv;
+ iwp[nrow + j] = nrow + j + 1;
+ iwp[nrow + j + neq * neq] = nrow + j + 1;
}
- nrow += *nequations;
+ nrow += neq;
actual[i] = ysave;
actualP[i] = ypsave;
- iwp[i] = i;
- iwp[i + *nequations] = i;
}
free(e);
int info = 0;
int *ipiv = NULL;
- ipiv = (int*)malloc(sizeof(int) * *nequations);
+ ipiv = (int *) malloc(*nequations * sizeof(int));
C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
// The terms are also available at
// 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.
//C The initial value problem is..
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 bugmes();quit;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]')
y0=1;t=2:6;t0=1;y0d=3;
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 bugmes();quit;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 bugmes();quit;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 bugmes();quit;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 bugmes();quit;end
+info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
//C
//C-----------------------------------------------------------------------
//C Second problem (Van Der Pol oscillator).
! 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');
-
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
// 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.
[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]')
[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).
' 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');
projects / scilab.git / commitdiff