%MPTEST % This script contains various tests of the mp*.m files % Note that MPDER, MPSHIFT, MPNEXT are tested incidentally, % with the major stress on MPVAL, MPMAK, MPAPI. % cb aug97 % try trivial stuff first nn=0; nn=nn+1 % univariate polynomial p(x) = 1-x; mp = mpmak([1 -1],1); if(any(mpval(mp,[-1:2])-[2:-1:-1])) fprintf('error in test1 in MVTEST.\n'), end if any(mpbrk(mpder(mp,1),'co')-(-1)) fprintf('error in test1der in MVTEST.\n'), end nn=nn+1 % univariate polynomial p(x) = (1-x)^2; mp = mpmak([1 -2 1],1); if(any(mpval(mp,[-1:2])-(1-[-1:2]).^2)) fprintf('error in test2 in MVTEST.\n'), end nn=nn+1 % univariate polynomial p(x) = (1-x)^3; mp = mpmak([1 -3 3 -1],1); if(any(mpval(mp,[-1:2])-(1-[-1:2]).^3)) fprintf('error in test3 in MVTEST.\n'), end if any(mpbrk(mpder(mp,-2),'co')-[6 -12 6]) fprintf('error in test3der in MVTEST.\n'), end nn=nn+1 % a bivariate polynomial p(x,y) = (1-x)^2; mp = mpmak([1 0 -2 0 0 1],2); if(any(mpval(mp,dcube(2))-[0 0 4 4])) fprintf('error in test4 in MVTEST.\n'), end nn=nn+1 % a trivariate quadratic that vanishes at all the corners of the cube mp = mpmak([-3 0 0 0 1 0 1 0 0 1],3); if(any(mpval(mp,dcube(3)))) fprintf('error in test5 in MVTEST.\n'), end p = [1;-2;3]; xi = [2;1;-2]; h = 1.e-7; if abs(mpval(mp,[p+h*xi p-h*xi])*[1;-1]/(2*h) - mpval(mpder(mp,xi),p))>h fprintf('error in test5der in MVTEST.\n'), end nn=nn+1 % finally, check it out when there is less symmetry % i.e., the polynomial (x-1)(y-1)(z-1) which vanishes at all but one corner mp = mpmak([-1 1 1 1 0 -.5 0 -.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0],3); if(abs(mpval(mp,dcube(3))-[0 0 0 0 0 0 0 -8])>1.e-15) fprintf('error in test6 in MVTEST.\n'), end p = [1;-2;3]; xi = [2;1;-2]; h = 1.e-5; if abs(mpval(mp,[p+h*xi p p-h*xi])*[1;-2;1]/(h*h) - ... mpval(mpder(mp,[xi xi]),p))>10*h fprintf('error in test6der in MVTEST.\n'), end if mpbrk(mpder(mp,[xi -2*xi xi -xi 3*xi]),'co') fprintf('error in test6ader in MVTEST.\n'), end % Now check the vector-valued calculation: % (x-1)(y-1)(z-1) % (x+1)(y-1)(z-1) % (x-1)(y+1)(z-1) % (x-1)(y-1)(z+1) nn=nn+1 [-1 +1 +1 +1 0 -.5 0 -.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0 1 -1 -1 +1 0 +.5 0 -.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0 1 -1 +1 -1 0 -.5 0 +.5 -.5 0 0 0 0 0 0 1/6 0 0 0 0 1 +1 -1 -1 0 -.5 0 -.5 +.5 0 0 0 0 0 0 1/6 0 0 0 0]; mp = mpmak(ans,3); [0 0 0 0 0 0 0 -8 0 0 0 8 0 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8 0]; if(any(abs(mpval(mp,dcube(3))-ans)>1.e-15)) fprintf('error in test7 in MVTEST.\n'), end nn=nn+1 nn=nn+1 nn=nn+1 % Now we are ready to tackle the least interpolant... % first, the simplest case mp = mpapi(0,0); if(mpval(mp,0)) fprintf('error in test8 in MVTEST.\n'), end nn=nn+1 % try bivariate, linear interpolation points = [0 1 0;0 0 1]; values = [0 1 1]; mp = mpapi(points, values); if(mpval(mp,points) - [0 1 1]) fprintf('error in test9 in MVTEST.\n'), end nn=nn+1 % try bivariate, quadratic, still simple points = [0 1 2 0 1 0; 0 0 0 1 1 2]; coefs = [1 -1 1 -1 1 -1]; mp0 = mpmak(coefs,2); if any(abs(mpbrk(mpapi(points,mpval(mp0,points)),'co')-coefs)>1.e-14) fprintf('error in test10 in MVTEST.\n'), end nn=nn+1 % try bivariate, cubic, still simple points = [0 1 2 3 0 1 2 0 1 0;0 0 0 0 1 1 1 2 2 3]; coefs = [1 -1 1 -1 1 -1 1 -1 1 -1]; mp0 = mpmak(coefs,2); if any(abs(mpbrk(mpapi(points,mpval(mp0,points)),'co')-coefs)>1.e-14) fprintf('error in test11 in MVTEST.\n'), end nn=nn+1 % try trivariate interpolation to a linear function points = [0 1 0 0 1; 0 0 1 0 1; 0 0 0 1 1]; coefs = [1 -2 3 -1]; mp0 = mpmak(coefs,3); values = mpval(mp0, points); if any(abs(mpbrk(mpapi(points, values, 1.e-14, [0;0;0]),'co')- ... [coefs 0 0 0 0 0 0])>1.e-14) fprintf('error in test12 in MVTEST.\n'), end nn=nn+1 % try trivariate, quadratic, interpolation to a quadratic function points = [0 1 2 0 1 0 0 1 0 0;0 0 0 1 1 2 0 0 1 0;0 0 0 0 0 0 1 1 1 2]; coefs = [1 -1 1 -1 1 -1 1 -1 1 -1]; mp0 = mpmak(coefs,3); values = mpval(mp0, points); mp = mpapi(points, values); mps = mpshift(mp,[0;0;0]); if any(abs(mpbrk(mps,'co')-coefs)>1.e-14) fprintf('error in test13 in MVTEST.\n'), end if any(abs(mpbrk(mp,'co') - mpbrk(mpshift(mps,mpbrk(mp,'ce')),'co'))>1.e-14) fprintf('error in test13shift in MVTEST.\n'), end % try trivariate, cubic, interpolation to a quadratic function points = dcube(3); coefs = [1 -1 1 -1 1 -1 1 -1 1 -1]; mp0 = mpmak(coefs,3); values = mpval(mp0, points); % In this example, the coefficients are not reproduced, but the interpolant is % in Pi_2. To check it cc = mpbrk(mpapi(points, values),'co'); cc(1:10) = cc(1:10)-coefs; mpd = mpmak(cc,3); if any(abs(mpval(mpd,points))>1.e-14) fprintf('error in test14 in MVTEST.\n'), end % re-use not implemented for v4 %nn=nn+1 %% test the re-use of factorization information: %points = rand(3,10); values = [1 -2 1]*(points.*(points - 1)); %[mp, mpfactor] = mpapi(points,values,1.e-12,[0;0;0]); %mpf = mpapi(mpfactor,values); %if (max(abs(mpbrk(mp,'co')-mpbrk(mpf,'co')))) % fprintf('error in test15 in MVTEST.\n'), end nn=nn+1 % test the vector capability: points = rand(3,10); values = (points.*(points - 1)); mp = mpapi(points,values); if any(abs(values - mpval(mp,points))>1.e-14) fprintf('error in test16 in MVTEST.\n'), end nn=nn+1 % test the noise level fprintf(['The next test checks interpolation at 10 data sites along a',... ' straight line in R^3 to quartic values.\n',... 'The interpolant should formally be of degree 9 but its coefficients\n',... 'of degree > 4 should be zero. Also, the derivative in any direction\n',... 'perpendicular to that straight line should be zero.\n',... 'The next three numbers give the size of the supposedly zero coefficients\n',... 'relative to the size of the coefficients of the interpolant.\n']) t = 6*rand(1,10)-3; xi = [1;2;3]; points = xi*t; values = t.^4; mp = mpapi(points,values); % coefficients of degree > 4 should be zero: coefs = mpbrk(mp,'co'); dimpjd = mpbrk(mp,'p'); mc = max(abs(coefs)); mc5 = max(abs(coefs(dimpjd(6)+1:dimpjd(11)))); fprintf('noise ratio = %8g \n', mc5/mc) % Also, the interpolant should be constant in directions orthogonal to xi : noise1 = max(abs(mpbrk(mpder(mp,[1;-2;1]),'co'))); noise2 = max(abs(mpbrk(mpder(mp,[3;0;-1]),'co'))); fprintf('derivative noise = %8g, %8g \n',noise1/mc, noise2/mc) % nn=nn+1 % % Now add a speed test, to compare the old bivariate routine with the present % % one, on a large number of points: % nn = 8; % for n=1:4 % nn = nn*2 % points = rand(2,nn); values = mpval(mpmak([1 -1 1 -1 1 -1],2),points); % tic % mp = mpapi(points,values,1.e-8); % toc % tic % bp = bpapi(points,values); % toc % end