add tests
Allan Cornet [Thu, 7 Feb 2008 13:46:25 +0000 (13:46 +0000)]
13 files changed:
scilab/modules/differential_equations/tests/unit_tests/bvode.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/bvode.tst [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/impl.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/impl.tst [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/int2d.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/int2d.tst [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/int3d.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/int3d.tst [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/intg.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/intg.tst [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/matode.dia.ref
scilab/modules/differential_equations/tests/unit_tests/ode.dia.ref [new file with mode: 0644]
scilab/modules/differential_equations/tests/unit_tests/ode.tst [new file with mode: 0644]

diff --git a/scilab/modules/differential_equations/tests/unit_tests/bvode.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/bvode.dia.ref
new file mode 100644 (file)
index 0000000..80cde27
--- /dev/null
@@ -0,0 +1,25 @@
+// to check that bvode works
+deff('df=dfsub(x,z)','df=[0,0,-6/x**2,-6/x]')
+deff('f=fsub(x,z)','f=(1 -6*x**2*z(4)-6*x*z(3))/x**3')
+deff('g=gsub(i,z)','g=[z(1),z(3),z(1),z(3)];g=g(i)')
+deff('dg=dgsub(i,z)',['dg=[1,0,0,0;0,0,1,0;1,0,0,0;0,0,1,0]';
+                      'dg=dg(i,:)'])
+deff('[z,mpar]=guess(x)','z=0;mpar=0')// unused here
+ //define trusol for testing purposes
+deff('u=trusol(x)',[
+   'u=0*ones(4,1)';
+   'u(1) =  0.25*(10*log(2)-3)*(1-x) + 0.5 *( 1/x   + (3+x)*log(x) - x)'
+   'u(2) = -0.25*(10*log(2)-3)       + 0.5 *(-1/x^2 + (3+x)/x      + log(x) - 1)'
+   'u(3) = 0.5*( 2/x^3 + 1/x   - 3/x^2)'
+   'u(4) = 0.5*(-6/x^4 - 1/x/x + 6/x^3)'])
+fixpnt=0;m=4;
+ncomp=1;aleft=1;aright=2;
+zeta=[1,1,2,2];
+ipar=zeros(1,11);
+ipar(3)=1;ipar(4)=2;ipar(5)=2000;ipar(6)=200;ipar(7)=1;
+ltol=[1,3];tol=[1.e-11,1.e-11];
+res=aleft:0.1:aright;
+z=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
+ fsub,dfsub,gsub,dgsub,guess);
+z1=[];for x=res,z1=[z1,trusol(x)]; end;
+if norm(z-z1) > 1D8*%eps then bugmes();quit;end;
diff --git a/scilab/modules/differential_equations/tests/unit_tests/bvode.tst b/scilab/modules/differential_equations/tests/unit_tests/bvode.tst
new file mode 100644 (file)
index 0000000..74cd052
--- /dev/null
@@ -0,0 +1,27 @@
+// to check that bvode works\r
+deff('df=dfsub(x,z)','df=[0,0,-6/x**2,-6/x]')\r
+deff('f=fsub(x,z)','f=(1 -6*x**2*z(4)-6*x*z(3))/x**3')\r
+deff('g=gsub(i,z)','g=[z(1),z(3),z(1),z(3)];g=g(i)')\r
+deff('dg=dgsub(i,z)',['dg=[1,0,0,0;0,0,1,0;1,0,0,0;0,0,1,0]';\r
+                      'dg=dg(i,:)'])\r
+deff('[z,mpar]=guess(x)','z=0;mpar=0')// unused here\r
+\r
+ //define trusol for testing purposes\r
+deff('u=trusol(x)',[ \r
+   'u=0*ones(4,1)';\r
+   'u(1) =  0.25*(10*log(2)-3)*(1-x) + 0.5 *( 1/x   + (3+x)*log(x) - x)'\r
+   'u(2) = -0.25*(10*log(2)-3)       + 0.5 *(-1/x^2 + (3+x)/x      + log(x) - 1)'\r
+   'u(3) = 0.5*( 2/x^3 + 1/x   - 3/x^2)'\r
+   'u(4) = 0.5*(-6/x^4 - 1/x/x + 6/x^3)'])\r
+\r
+fixpnt=0;m=4;\r
+ncomp=1;aleft=1;aright=2;\r
+zeta=[1,1,2,2];\r
+ipar=zeros(1,11);\r
+ipar(3)=1;ipar(4)=2;ipar(5)=2000;ipar(6)=200;ipar(7)=1;\r
+ltol=[1,3];tol=[1.e-11,1.e-11];\r
+res=aleft:0.1:aright;\r
+z=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...\r
+ fsub,dfsub,gsub,dgsub,guess);\r
+z1=[];for x=res,z1=[z1,trusol(x)]; end;  \r
+if norm(z-z1) > 1D8*%eps then pause,end;
\ No newline at end of file
diff --git a/scilab/modules/differential_equations/tests/unit_tests/impl.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/impl.dia.ref
new file mode 100644 (file)
index 0000000..a684c33
--- /dev/null
@@ -0,0 +1,3 @@
+y=impl([1;0;0],[-0.04;0.04;0],0,0.4,'resid','aplusp');
+if size(y) <> [3 1] then bugmes();quit;end
+if norm( y - [0.9851721;0.0000339;0.0147941]) < %eps then bugmes();quit;end
diff --git a/scilab/modules/differential_equations/tests/unit_tests/impl.tst b/scilab/modules/differential_equations/tests/unit_tests/impl.tst
new file mode 100644 (file)
index 0000000..cd9d36e
--- /dev/null
@@ -0,0 +1,3 @@
+y=impl([1;0;0],[-0.04;0.04;0],0,0.4,'resid','aplusp');\r
+if size(y) <> [3 1] then pause,end\r
+if norm( y - [0.9851721;0.0000339;0.0147941]) < %eps then pause,end
\ No newline at end of file
diff --git a/scilab/modules/differential_equations/tests/unit_tests/int2d.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/int2d.dia.ref
new file mode 100644 (file)
index 0000000..7967e10
--- /dev/null
@@ -0,0 +1,6 @@
+X=[0,0;1,1;1,0];
+Y=[0,0;0,1;1,1];
+deff('z=f(x,y)','z=cos(x+y)')
+[I,e]=int2d(X,Y,f);
+if e <%eps then bugmes();quit;end;
+if abs(I-0.5) < %eps then bugmes();quit;end;
diff --git a/scilab/modules/differential_equations/tests/unit_tests/int2d.tst b/scilab/modules/differential_equations/tests/unit_tests/int2d.tst
new file mode 100644 (file)
index 0000000..74c6690
--- /dev/null
@@ -0,0 +1,6 @@
+X=[0,0;1,1;1,0];\r
+Y=[0,0;0,1;1,1];\r
+deff('z=f(x,y)','z=cos(x+y)')\r
+[I,e]=int2d(X,Y,f);\r
+if e <%eps then pause,end;\r
+if abs(I-0.5) < %eps then pause,end;
\ No newline at end of file
diff --git a/scilab/modules/differential_equations/tests/unit_tests/int3d.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/int3d.dia.ref
new file mode 100644 (file)
index 0000000..9a5bb60
--- /dev/null
@@ -0,0 +1,28 @@
+X=[0;1;0;0];
+Y=[0;0;1;0];
+Z=[0;0;0;1];
+[RESULT,ERROR]=int3d(X,Y,Z,'int3dex');
+// computes the integrand exp(x*x+y*y+z*z) over the
+//tetrahedron (0.,0.,0.),(1.,0.,0.),(0.,1.,0.),(0.,0.,1.)
+//integration over a cube  -1<=x<=1;-1<=y<=1;-1<=z<=1
+//  bottom  -top-     right    -left-   front   -rear-
+X=[ 0, 0,    0, 0,    0, 0,    0, 0,    0, 0,    0, 0;
+   -1,-1,   -1,-1,    1, 1,   -1,-1,   -1,-1,   -1,-1;
+    1,-1,    1,-1,    1, 1,   -1,-1,    1,-1,    1,-1;
+    1, 1,    1, 1,    1, 1,   -1,-1,    1, 1,    1, 1];
+Y=[ 0, 0,    0, 0,    0, 0,    0, 0,    0, 0,    0, 0;
+   -1,-1,   -1,-1,   -1, 1,   -1, 1,   -1,-1,    1, 1;
+   -1, 1,   -1, 1,    1, 1,    1, 1,   -1,-1,    1, 1;
+    1, 1,    1, 1,   -1,-1,   -1,-1,   -1,-1,    1, 1];
+Z=[ 0, 0,    0, 0,    0, 0,    0, 0,    0, 0,    0, 0;
+   -1,-1,    1, 1,   -1, 1,   -1, 1,   -1,-1,   -1,-1;
+   -1,-1,    1, 1,   -1,-1,   -1,-1,   -1, 1,   -1, 1;
+   -1,-1,    1, 1,    1, 1,    1, 1,    1, 1,    1, 1];
+function v=f1(xyz,numfun),v=exp(xyz'*xyz),endfunction
+[result,err]=int3d(X,Y,Z,f1,1,[0,100000,1.d-5,1.d-7]);
+if result > 26 then bugmes();quit;end
+if err <%eps then bugmes();quit;end
+function v=f2(xyz,numfun),v=1,endfunction
+[result,err]=int3d(X,Y,Z,f2,1,[0,100000,1.d-5,1.d-7]);
+if result <> 8 then bugmes();quit;end
+if err <%eps then bugmes();quit;end
diff --git a/scilab/modules/differential_equations/tests/unit_tests/int3d.tst b/scilab/modules/differential_equations/tests/unit_tests/int3d.tst
new file mode 100644 (file)
index 0000000..5c4e0e4
--- /dev/null
@@ -0,0 +1,35 @@
+X=[0;1;0;0];\r
+Y=[0;0;1;0];\r
+Z=[0;0;0;1];\r
+[RESULT,ERROR]=int3d(X,Y,Z,'int3dex');\r
+// computes the integrand exp(x*x+y*y+z*z) over the \r
+//tetrahedron (0.,0.,0.),(1.,0.,0.),(0.,1.,0.),(0.,0.,1.)\r
+\r
+\r
+//integration over a cube  -1<=x<=1;-1<=y<=1;-1<=z<=1\r
+\r
+//  bottom  -top-     right    -left-   front   -rear- \r
+X=[ 0, 0,    0, 0,    0, 0,    0, 0,    0, 0,    0, 0;          \r
+   -1,-1,   -1,-1,    1, 1,   -1,-1,   -1,-1,   -1,-1; \r
+    1,-1,    1,-1,    1, 1,   -1,-1,    1,-1,    1,-1;     \r
+    1, 1,    1, 1,    1, 1,   -1,-1,    1, 1,    1, 1];         \r
+Y=[ 0, 0,    0, 0,    0, 0,    0, 0,    0, 0,    0, 0; \r
+   -1,-1,   -1,-1,   -1, 1,   -1, 1,   -1,-1,    1, 1;\r
+   -1, 1,   -1, 1,    1, 1,    1, 1,   -1,-1,    1, 1;   \r
+    1, 1,    1, 1,   -1,-1,   -1,-1,   -1,-1,    1, 1]; \r
+Z=[ 0, 0,    0, 0,    0, 0,    0, 0,    0, 0,    0, 0;\r
+   -1,-1,    1, 1,   -1, 1,   -1, 1,   -1,-1,   -1,-1; \r
+   -1,-1,    1, 1,   -1,-1,   -1,-1,   -1, 1,   -1, 1;  \r
+   -1,-1,    1, 1,    1, 1,    1, 1,    1, 1,    1, 1];      \r
+\r
+function v=f1(xyz,numfun),v=exp(xyz'*xyz),endfunction\r
+[result,err]=int3d(X,Y,Z,f1,1,[0,100000,1.d-5,1.d-7]);\r
+\r
+if result > 26 then pause,end\r
+if err <%eps then pause,end\r
+\r
+function v=f2(xyz,numfun),v=1,endfunction\r
+[result,err]=int3d(X,Y,Z,f2,1,[0,100000,1.d-5,1.d-7]);\r
+\r
+if result <> 8 then pause,end\r
+if err <%eps then pause,end\r
diff --git a/scilab/modules/differential_equations/tests/unit_tests/intg.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/intg.dia.ref
new file mode 100644 (file)
index 0000000..0a6772e
--- /dev/null
@@ -0,0 +1,5 @@
+// to check that intg works
+function y=f(x),y=x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)),endfunction
+exact=-2.5432596188;
+I=intg(0,2*%pi,f);
+if abs(exact-I) > 1 then bugmes();quit;end
diff --git a/scilab/modules/differential_equations/tests/unit_tests/intg.tst b/scilab/modules/differential_equations/tests/unit_tests/intg.tst
new file mode 100644 (file)
index 0000000..c5f0705
--- /dev/null
@@ -0,0 +1,5 @@
+// to check that intg works\r
+function y=f(x),y=x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)),endfunction\r
+exact=-2.5432596188;\r
+I=intg(0,2*%pi,f);\r
+if abs(exact-I) > 1 then pause,end
\ No newline at end of file
index 778e32e..4fb9241 100644 (file)
@@ -55,7 +55,7 @@ if (y6-yref) > 2*0.00001 then bugmes();quit;end
 //
 a=rand(3,3);ea=expm(a);
 deff('[ydot]=f(t,y)','ydot=a*y')
-Warning : redefining function : f                       
+Warning : redefining function: f                       
 
 t1=1;y=ode('adams',eye(a),t0,t1,f);
 if maxi(ea-y) > Leps then bugmes();quit;end
@@ -81,10 +81,10 @@ deff('[p]=dgbydy(t,y,s)','[-0.04,1.d4*y(3),1.d4*y(2);...
 rand('seed',0);rand('normal');
 nx=20;A=rand(nx,nx);A=A-4.5*eye();
 deff('y=f(t,x)','y=A*x')
-Warning : redefining function : f                       
+Warning : redefining function: f                       
 
 deff('J=j(t,x)','J=A')
-Warning : redefining function : j                       
+Warning : redefining function: j                       
 
 x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
 nt=size(t,'*');
@@ -99,7 +99,7 @@ xf=ode(x0,t0,t,f);   //lsoda
 if norm(xf(:,nt)-eAt*x0)>Leps then bugmes();quit;end
 xfj=ode(x0,t0,t,f,j);   //lsoda with jacobian
 Warning: Jacobian external is given, but
- not used!,  see %ODEOPTIONS(6).
+ not used!,  see ODEOPTIONS(6).
 
 if norm(xfj(:,nt)-eAt*x0)>Leps then bugmes();quit;end
 //itask=2;   --->  solution at mesh points  ---> t=tmax
@@ -123,7 +123,7 @@ deff('y=fcrit(t,x)',['if t<=tcrit then'
                       'else'
                       ' y=A*x;chk=resume(1);end'])
 x42=ode(x0,t0,t,fcrit);
-Warning : integration up to tcrit.
+Warning: integration up to tcrit.
 
 if chk==1 then bugmes();quit;end
 [p,q]=size(x42);
@@ -150,14 +150,14 @@ if norm(x61-xf,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=0; // jacobian nor called nor estimated
 x60=ode('st',x0,t0,t,f,j);   //Jacobian not used (warning)
 Warning: Jacobian external is given, but
- not used!,  see %ODEOPTIONS(6).
+ not used!,  see ODEOPTIONS(6).
 
 x60=ode('st',x0,t0,t,f);    //Jacobian not used
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=1;//Jacobian estimated
 x60=ode('st',x0,t0,t,f)  ;
 Warning: No Jacobian external given but
- one is required by %ODEOPTIONS(6) value !
+ one is required by ODEOPTIONS(6) value !
 
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 //test of %ODEOPTIONS(6)=jacflag   (adams)
@@ -167,7 +167,7 @@ if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=0;// jacobian nor called nor estimated
 x60=ode('ad',x0,t0,t,f,j);   //Jacobian not used (warning)
 Warning: Jacobian external is given, but
- not used!,  see %ODEOPTIONS(6).
+ not used!,  see ODEOPTIONS(6).
 
 x60=ode('ad',x0,t0,t,f);    //Jacobian not used
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
@@ -178,7 +178,7 @@ if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=1;//Jacobian estimated
 x60=ode(x0,t0,t,f);
 Warning: No Jacobian external given but
- one is required by %ODEOPTIONS(6) value !
+ one is required by ODEOPTIONS(6) value !
 
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 //   Banded Jacobian
@@ -195,7 +195,7 @@ for k=1:nx-1, A(k+1,k)=-1;end
 for k=1:nx-2, A(k+2,k)=2;end
 for k=1:nx-3, A(k+3,k)=-3;end
 deff('xd=f(t,x)','xd=A*x')
-Warning : redefining function : f                       
+Warning : redefining function: f                       
 
 ml=3;mu=2;
 %ODEOPTIONS(11:12)=[ml,mu];
@@ -203,7 +203,7 @@ for i=1:nx;
     for j=1:nx;
 if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
 end;end;
-Warning : redefining function : j                       
+Warning : redefining function: j                       
 
 // J is a ml+mu+1 x ny matrix.
 // Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
@@ -223,7 +223,7 @@ if norm(xnotband-xband,1)>Leps then bugmes();quit;end
 deff('jj=j(t,x)','jj=J')
 xband=ode('st',x0,t0,t,f,j);
 Warning: Jacobian external is given, but
- not used!,  see %ODEOPTIONS(6).
+ not used!,  see ODEOPTIONS(6).
 
 if norm(xnotband-xband,1)>Leps then bugmes();quit;end
 //            Test of %ODEOPTIONS(7)
@@ -269,7 +269,7 @@ A=diag([-10,-0.01,-1]);
 deff('uu=u(t)','uu=sin(t)');
 B=rand(3,1);
 deff('y=f(t,x)','y=A*x+B*u(t)')
-Warning : redefining function : f                       
+Warning : redefining function: f                       
 
 %ODEOPTIONS(1)=2;
 yy1=ode('stiff',[1;1;1],0,1,f);
@@ -281,10 +281,10 @@ clear %ODEOPTIONS;
 rand('seed',0);rand('normal');
 nx=20;A=rand(nx,nx);A=A-4.5*eye();
 deff('y=f(t,x)','y=A*x')
-Warning : redefining function : f                       
+Warning : redefining function: f                       
 
 deff('J=j(t,x)','J=A')
-Warning : redefining function : j                       
+Warning : redefining function: j                       
 
 //%ODEOPTIONS(1)=1;
 y2=ode('stiff',ones(nx,1),0,2,f,j);
@@ -297,7 +297,7 @@ yaj=ode('adams',ones(nx,1),0,2,f,j);
 ysf=ode('stiff',ones(nx,1),0,2,f,j);
 ysj=ode('stiff',ones(nx,1),0,2,f,j);
 deff('xd=f(t,x)','xd=A*x+B*sin(3*t)')
-Warning : redefining function : f                       
+Warning : redefining function: f                       
 
 A=rand(10,10)-4.5*eye();B=rand(10,1);
 x=ode(ones(10,1),0,[1,2,3],f);
diff --git a/scilab/modules/differential_equations/tests/unit_tests/ode.dia.ref b/scilab/modules/differential_equations/tests/unit_tests/ode.dia.ref
new file mode 100644 (file)
index 0000000..eabe735
--- /dev/null
@@ -0,0 +1,7 @@
+ // to check that ode works
+ // ---------- Simple one dimension ODE (Scilab function external)
+ // dy/dt=y^2-y sin(t)+cos(t), y(0)=0
+ function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
+ y0=0;t0=0;t=0:0.1:%pi;
+ y=ode(y0,t0,t,f);
+ if size(y) <> [1 32] then bugmes();quit;end
diff --git a/scilab/modules/differential_equations/tests/unit_tests/ode.tst b/scilab/modules/differential_equations/tests/unit_tests/ode.tst
new file mode 100644 (file)
index 0000000..0e9dd8a
--- /dev/null
@@ -0,0 +1,7 @@
+ // to check that ode works\r
+ // ---------- Simple one dimension ODE (Scilab function external)\r
+ // dy/dt=y^2-y sin(t)+cos(t), y(0)=0\r
+ function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction\r
+ y0=0;t0=0;t=0:0.1:%pi;\r
+ y=ode(y0,t0,t,f);\r
+ if size(y) <> [1 32] then pause,end
\ No newline at end of file