Sylvestre Ledru [Wed, 19 Sep 2012 08:38:21 +0000 (10:38 +0200)]
Change-Id: I23b12d92e1c30c412f86368494fb17f31c92bbab

12 files changed:

index 7a3ee17..03641ae 100644 (file)
@@ -79,6 +79,25 @@ xtitle("Cubic Shepard Interpolation of cos(x)cos(y) with randomly chosen interpo
legends("interpolation points",-9,1)
show_window()
]]></programlisting>
+        <scilab:image>
+            n = 150; // nb of interpolation points
+            xy = grand(n,2,"unf",0,2*%pi);
+            z = cos(xy(:,1)).*cos(xy(:,2));
+            xyz = [xy z];
+            tl_coef = cshep2d(xyz);
+
+            // evaluation on a grid
+            m = 30;
+            xx = linspace(0,2*%pi,m);
+            [X,Y] = ndgrid(xx,xx);
+            Z = eval_cshep2d(X,Y, tl_coef);
+            clf()
+            plot3d(xx,xx,Z,flag=[2 6 4])
+            param3d1(xy(:,1),xy(:,2),list(z,-9), flag=[0 0])
+            xtitle("Cubic Shepard Interpolation of cos(x)cos(y) with randomly chosen interpolation points")
+            legends("interpolation points",-9,1)
+            show_window()
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
index 5582c44..7156ce1 100644 (file)
@@ -103,15 +103,38 @@ zp = eval_cshep2d(xp,yp,S);
// compute facet (to draw one color for extrapolation region
// and another one for the interpolation region)
[xf,yf,zf] = genfac3d(xx,xx,zp);
-color = 2*ones(1,size(zf,2));
+colors = 2*ones(1,size(zf,2));
// indices corresponding to facet in the interpolation region
ind=find( mean(xf,"r")>0 & mean(xf,"r")<1 & mean(yf,"r")>0 & mean(yf,"r")<1 );
-color(ind)=3;
+colors(ind)=3;
clf();
-plot3d(xf,yf,list(zf,color), flag=[2 6 4])
+plot3d(xf,yf,list(zf,colors), flag=[2 6 4])
legends(["extrapolation region","interpolation region"],[2 3],1)
show_window()
]]></programlisting>
+        <scilab:image><![CDATA[
+deff("z=f(x,y)","z = 1+ 50*(x.*(1-x).*y.*(1-y)).^2")
+x = linspace(0,1,10);
+[X,Y] = ndgrid(x,x);
+X = X(:); Y = Y(:); Z = f(X,Y);
+S = cshep2d([X Y Z]);
+// evaluation inside and outside the square [0,1]x[0,1]
+m = 40;
+xx = linspace(-1.5,0.5,m);
+[xp,yp] = ndgrid(xx,xx);
+zp = eval_cshep2d(xp,yp,S);
+// compute facet (to draw one color for extrapolation region
+// and another one for the interpolation region)
+[xf,yf,zf] = genfac3d(xx,xx,zp);
+colors = 2*ones(1,size(zf,2));
+// indices corresponding to facet in the interpolation region
+ind=find( mean(xf,"r")>0 & mean(xf,"r")<1 & mean(yf,"r")>0 & mean(yf,"r")<1 );
+colors(ind)=3;
+clf();
+plot3d(xf,yf,list(zf,colors), flag=[2 6 4])
+legends(["extrapolation region","interpolation region"],[2 3],1)
+show_window()
+]]></scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
index df6c2f4..568cd88 100644 (file)
@@ -171,7 +171,33 @@ subplot(3,1,3)
plot2d(xx, [yy2k yy2f])
legends(["not_a_knot spline","fast sub-spline"], [1 2], "lr",%f)
xtitle("spline interpolation (second derivatives)")
-
+ ]]></programlisting>
+        <scilab:image>
+            a = -8; b = 8;
+            x = linspace(a,b,20)';
+            y = sinc(x);
+            dk = splin(x,y);  // not_a_knot
+            df = splin(x,y, "fast");
+            xx = linspace(a,b,800)';
+            [yyk, yy1k, yy2k] = interp(xx, x, y, dk);
+            [yyf, yy1f, yy2f] = interp(xx, x, y, df);
+            clf()
+            subplot(3,1,1)
+            plot2d(xx, [yyk yyf])
+            plot2d(x, y, style=-9)
+            legends(["not_a_knot spline","fast sub-spline","interpolation points"],...
+            [1 2 -9], "ur",%f)
+            xtitle("spline interpolation")
+            subplot(3,1,2)
+            plot2d(xx, [yy1k yy1f])
+            legends(["not_a_knot spline","fast sub-spline"], [1 2], "ur",%f)
+            xtitle("spline interpolation (derivatives)")
+            subplot(3,1,3)
+            plot2d(xx, [yy2k yy2f])
+            legends(["not_a_knot spline","fast sub-spline"], [1 2], "lr",%f)
+            xtitle("spline interpolation (second derivatives)")
+        </scilab:image>
+        <programlisting role="example"><![CDATA[
// here is an example showing the different extrapolation possibilities
x = linspace(0,1,11)';
y = cosh(x-0.5);
@@ -185,6 +211,19 @@ clf()
plot2d(xx,[yy0 yy1 yy2 yy3],style=2:5,frameflag=2,leg="C0@linear@natural@periodic")
xtitle(" different way to evaluate a spline outside its domain")
]]></programlisting>
+        <scilab:image>
+            x = linspace(0,1,11)';
+            y = cosh(x-0.5);
+            d = splin(x,y);
+            xx = linspace(-0.5,1.5,401)';
+            yy0 = interp(xx,x,y,d,"C0");
+            yy1 = interp(xx,x,y,d,"linear");
+            yy2 = interp(xx,x,y,d,"natural");
+            yy3 = interp(xx,x,y,d,"periodic");
+            clf()
+            plot2d(xx,[yy0 yy1 yy2 yy3],style=2:5,frameflag=2,leg="C0@linear@natural@periodic")
+            xtitle(" different way to evaluate a spline outside its domain")
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
@@ -152,6 +152,17 @@ plot(xx,[yy1;yy2;yy3],x,y,'*')
xtitle('interpolation of square function')
legend(['linear','spline','nearest'],a=2)
]]></programlisting>
+        <scilab:image>
+            x=linspace(0,3,20);
+            y=x^2;
+            xx=linspace(0,3,100);
+            yy1=interp1(x,y,xx,'linear');
+            yy2=interp1(x,y,xx,'spline');
+            yy3=interp1(x,y,xx,'nearest');
+            plot(xx,[yy1;yy2;yy3],x,y,'*')
+            xtitle('interpolation of square function')
+            legend(['linear','spline','nearest'],a=2)
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
index 41aeb4e..4524c0f 100644 (file)
@@ -166,6 +166,35 @@ plot3d(xx, yy, zz4, flag=[2 6 4])
xtitle("extrapolation with the natural outmode")
show_window()
]]></programlisting>
+        <scilab:image>
+            n = 7;  // a n x n interpolation grid
+            x = linspace(0,2*%pi,n); y = x;
+            z = cos(x')*cos(y);
+            C = splin2d(x, y, z, "periodic");
+
+            // now evaluate on a bigger domain than [0,2pi]x [0,2pi]
+            m = 80; // discretisation parameter of the evaluation grid
+            xx = linspace(-0.5*%pi,2.5*%pi,m); yy = xx;
+            [XX,YY] = ndgrid(xx,yy);
+            zz1 = interp2d(XX,YY, x, y, C, "C0");
+            zz2 = interp2d(XX,YY, x, y, C, "by_zero");
+            zz3 = interp2d(XX,YY, x, y, C, "periodic");
+            zz4 = interp2d(XX,YY, x, y, C, "natural");
+            clf()
+            subplot(2,2,1)
+            plot3d(xx, yy, zz1, flag=[2 6 4])
+            xtitle("extrapolation with the C0 outmode")
+            subplot(2,2,2)
+            plot3d(xx, yy, zz2, flag=[2 6 4])
+            xtitle("extrapolation with the by_zero outmode")
+            subplot(2,2,3)
+            plot3d(xx, yy, zz3, flag=[2 6 4])
+            xtitle("extrapolation with the periodic outmode")
+            subplot(2,2,4)
+            plot3d(xx, yy, zz4, flag=[2 6 4])
+            xtitle("extrapolation with the natural outmode")
+            show_window()
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
index 63e8e1c..7e15953 100644 (file)
@@ -49,6 +49,13 @@ plot2d(x',y',[-3],"011"," ",[-10,-40,50,50]);
yi=interpln([x;y],-4:45);
plot2d((-4:45)',yi',[3],"000");
]]></programlisting>
+        <scilab:image>
+            x=[1 10 20 30 40];
+            y=[1 30 -10 20 40];
+            plot2d(x',y',[-3],"011"," ",[-10,-40,50,50]);
+            yi=interpln([x;y],-4:45);
+            plot2d((-4:45)',yi',[3],"000");
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
index 0bc4b65..23d8ffc 100644 (file)
@@ -149,7 +149,18 @@ clf()
plot2d(xx,yy,style=2)
plot2d(x,y,style=-9, strf="000")
xtitle("linear interpolation of sin(x) with 11 interpolation points")
-
+ ]]></programlisting>
+        <scilab:image>
+            x = linspace(0,2*%pi,11);
+            y = sin(x);
+            xx = linspace(-2*%pi,4*%pi,400)';
+            yy = linear_interpn(xx, x, y, "periodic");
+            clf()
+            plot2d(xx,yy,style=2)
+            plot2d(x,y,style=-9, strf="000")
+            xtitle("linear interpolation of sin(x) with 11 interpolation points")
+        </scilab:image>
+        <programlisting role="example"><![CDATA[
// example 2 : bilinear interpolation
n = 8;
x = linspace(0,2*%pi,n); y = x;
@@ -164,7 +175,23 @@ param3d1(xg,yg, list(z,-9*ones(1,n)), flag=[0 0])
xtitle("Bilinear interpolation of 2sin(x)sin(y)")
legends("interpolation points",-9,1)
show_window()
-
+ ]]></programlisting>
+        <scilab:image>
+            n = 8;
+            x = linspace(0,2*%pi,n); y = x;
+            z = 2*sin(x')*sin(y);
+            xx = linspace(0,2*%pi, 40);
+            [xp,yp] = ndgrid(xx,xx);
+            zp = linear_interpn(xp,yp, x, y, z);
+            clf()
+            plot3d(xx, xx, zp, flag=[2 6 4])
+            [xg,yg] = ndgrid(x,x);
+            param3d1(xg,yg, list(z,-9*ones(1,n)), flag=[0 0])
+            xtitle("Bilinear interpolation of 2sin(x)sin(y)")
+            legends("interpolation points",-9,1)
+            show_window()
+        </scilab:image>
+        <programlisting role="example"><![CDATA[
// example 3 : bilinear interpolation and experimentation
//             with all the outmode features
nx = 20; ny = 30;
@@ -201,11 +228,12 @@ subplot(2,3,6)
plot3d(xp, yp, zp5, leg="x@y@z", flag = [2 4 4])
xtitle("by_nan")
show_window()
-
+ ]]></programlisting>
+        <programlisting role="example"><![CDATA[
// example 4 : trilinear interpolation (see splin3d help
//             page which have the same example with
//             tricubic spline interpolation)
-exec("SCI/modules/interpolation/demos/interp_demo.sci")
+exec("SCI/modules/interpolation/demos/interp_demo.sci");
func =  "v=(x-0.5).^2 + (y-0.5).^3 + (z-0.5).^2";
deff("v=f(x,y,z)",func);
n = 5;
index 1f283d4..e33cb9e 100644 (file)
@@ -133,6 +133,29 @@ plot2d(xd,[ye yd ys],style=[2 -2 3], ...
xtitle("a least square spline")
show_window()
]]></programlisting>
+        <scilab:image>
+            a = 0; b = 2*%pi;
+            sigma = 0.1;  // standard deviation of the gaussian noise
+            m = 200;       // number of experimental points
+            xd = linspace(a,b,m)';
+            yd = sin(xd) + grand(xd,"nor",0,sigma);
+
+            n = 6; // number of breakpoints
+            x = linspace(a,b,n)';
+
+            // compute the spline
+            [y, d] = lsq_splin(xd, yd, x);  // use equal weights
+
+            // plotting
+            ye = sin(xd);
+            ys = interp(xd, x, y, d);
+            clf()
+            plot2d(xd,[ye yd ys],style=[2 -2 3], ...
+            leg="exact function@experimental measures (gaussian perturbation)@fitted spline")
+            xtitle("a least square spline")
+            show_window()
+        </scilab:image>
+
</refsection>
<refsection role="see also">
<title>See Also</title>
index 1caa748..cc90efe 100644 (file)
@@ -50,6 +50,13 @@ plot2d(x',y',[3],"011"," ",[-10,-40,50,50]);
yi=smooth([x;y],0.1);
plot2d(yi(1,:)',yi(2,:)',[1],"000");
]]></programlisting>
+        <scilab:image>
+            x=[1 10 20 30 40];
+            y=[1 30 -10 20 40];
+            plot2d(x',y',[3],"011"," ",[-10,-40,50,50]);
+            yi=smooth([x;y],0.1);
+            plot2d(yi(1,:)',yi(2,:)',[1],"000");
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
@@ -212,7 +212,22 @@ clf()
plot2d(xx, [yyi yye], style=[2 5], leg="interpolation spline@exact function")
plot2d(x, y, -9)
xtitle("interpolation of the Runge function")
-
+ ]]></programlisting>
+        <scilab:image>
+            deff("y=runge(x)","y=1 ./(1 + x.^2)")
+            a = -5; b = 5; n = 11; m = 400;
+            x = linspace(a, b, n)';
+            y = runge(x);
+            d = splin(x, y);
+            xx = linspace(a, b, m)';
+            yyi = interp(xx, x, y, d);
+            yye = runge(xx);
+            clf()
+            plot2d(xx, [yyi yye], style=[2 5], leg="interpolation spline@exact function")
+            plot2d(x, y, -9)
+            xtitle("interpolation of the Runge function")
+        </scilab:image>
+        <programlisting role="example"><![CDATA[
// example 2 : show behavior of different splines on random data
a = 0; b = 1;        // interval of interpolation
n = 10;              // nb of interpolation  points
@@ -230,6 +245,23 @@ plot2d(x,y,-9, strf="000")  // to show interpolation points
xtitle("Various spline and sub-splines on random data")
show_window()
]]></programlisting>
+        <scilab:image>
+            a = 0; b = 1;        // interval of interpolation
+            n = 10;              // nb of interpolation  points
+            m = 800;             // discretisation for evaluation
+            x = linspace(a,b,n)'; // abscissae of interpolation points
+            y = rand(x);          // ordinates of interpolation points
+            xx = linspace(a,b,m)';
+            yk = interp(xx, x, y, splin(x,y,"not_a_knot"));
+            yf = interp(xx, x, y, splin(x,y,"fast"));
+            ym = interp(xx, x, y, splin(x,y,"monotone"));
+            clf()
+            plot2d(xx, [yf ym yk], style=[5 2 3], strf="121", ...
+            leg="fast@monotone@not a knot spline")
+            plot2d(x,y,-9, strf="000")  // to show interpolation points
+            xtitle("Various spline and sub-splines on random data")
+            show_window()
+        </scilab:image>
</refsection>
<refsection role="see also">
<title>See Also</title>
index c8ea5be..2eeac63 100644 (file)
@@ -144,7 +144,27 @@ plot3d(xx, yy, zz, flag=[2 4 4])
param3d1(X,Y,list(z,-9*ones(1,n)), flag=[0 0])
str = msprintf(" with %d x %d interpolation points. ermax = %g",n,n,emax)
xtitle("spline interpolation of cos(x)cos(y)"+str)
-
+ ]]></programlisting>
+        <scilab:image>
+            n = 7;  // a regular grid with n x n interpolation points
+            // will be used
+            x = linspace(0,2*%pi,n); y = x;
+            z = cos(x')*cos(y);
+            C = splin2d(x, y, z, "periodic");
+            m = 50; // discretisation parameter of the evaluation grid
+            xx = linspace(0,2*%pi,m); yy = xx;
+            [XX,YY] = ndgrid(xx,yy);
+            zz = interp2d(XX,YY, x, y, C);
+            emax = max(abs(zz - cos(xx')*cos(yy)));
+            clf()
+            plot3d(xx, yy, zz, flag=[2 4 4])
+            [X,Y] = ndgrid(x,y);
+            param3d1(X,Y,list(z,-9*ones(1,n)), flag=[0 0])
+            str = msprintf(" with %d x %d interpolation points. ermax = %g",n,n,emax)
+            xtitle("spline interpolation of cos(x)cos(y)"+str)
+        </scilab:image>
+
+        <programlisting role="example"><![CDATA[
// example 2 : different interpolation functions on random data
n = 6;
x = linspace(0,1,n); y = x;
@@ -171,7 +191,36 @@ subplot(2,2,4)
plot3d1(xp, yp, ZP4, flag=[2 2 4])
xtitle("monotone")
show_window()
-
+ ]]></programlisting>
+        <scilab:image>
+            // example 2 : different interpolation functions on random data
+            n = 6;
+            x = linspace(0,1,n); y = x;
+            z = rand(n,n);
+            np = 50;
+            xp = linspace(0,1,np); yp = xp;
+            [XP, YP] = ndgrid(xp,yp);
+            ZP1 = interp2d(XP, YP, x, y, splin2d(x, y, z, "not_a_knot"));
+            ZP2 = linear_interpn(XP, YP, x, y, z);
+            ZP3 = interp2d(XP, YP, x, y, splin2d(x, y, z, "natural"));
+            ZP4 = interp2d(XP, YP, x, y, splin2d(x, y, z, "monotone"));
+            xset("colormap", jetcolormap(64))
+            clf()
+            subplot(2,2,1)
+            plot3d1(xp, yp, ZP1, flag=[2 2 4])
+            xtitle("not_a_knot")
+            subplot(2,2,2)
+            plot3d1(xp, yp, ZP2, flag=[2 2 4])
+            xtitle("bilinear interpolation")
+            subplot(2,2,3)
+            plot3d1(xp, yp, ZP3, flag=[2 2 4])
+            xtitle("natural")
+            subplot(2,2,4)
+            plot3d1(xp, yp, ZP4, flag=[2 2 4])
+            xtitle("monotone")
+            show_window()
+        </scilab:image>
+        <programlisting role="example"><![CDATA[
// example 3 : not_a_knot spline and monotone sub-spline
//             on a step function
a = 0; b = 1; c = 0.25; d = 0.75;
@@ -196,6 +245,7 @@ subplot(1,2,2)
plot3d1(xp, xp, zp2, flag=[-2 6 4])
xtitle("subspline (monotone)")
]]></programlisting>
+
</refsection>
<refsection role="see also">
<title>See Also</title>
index 0052aa8..1b13e34 100644 (file)
@@ -95,6 +95,9 @@ vp_exact = f(xp,yp,zp);
vp_interp = interp3d(xp,yp,zp, tl);
er = max(abs(vp_exact - vp_interp))
// now retry with n=20 and see the error
+ ]]></programlisting>
+
+        <programlisting role="example"><![CDATA[

// example 2 (see linear_interpn help page which have the
//            same example with trilinear interpolation)