find roots of a spline
in [Matlab]
|
From: Andres on 22 Apr 2010 08:49 Hi, I'd like to speed up my code for the following problem: input: two vectors x, y (numel(x) = numel(y) = ~5e5) output: roots of the cubic spline through (x,y) I've written a function 'rootspline' to get the roots of the piecewise polynomial form of the cubic spline 'pp' (see below), so I run pp = spline(x,y); xo = rootspline(pp); compare example 3 in the rootspline doc. Any ideas how to improve rootspline, its subfunction rootsbnd, or how to do it in a different way? Looking at the profiler results, it seems that calling rootsbnd for each piece of the spline hurts most... Thanks for your thoughts! Andres >> function zo = rootspline(pp) % rootspline find roots of a spline % % xo = rootspline(pp) returns the unique roots of the piecewise polynomial % contained in pp. Typically, pp was constructed using the functions % interp1, pchip, spline, or the spline utility mkpp. % % rootspline uses the roots function, so the input must not contain NaN or % Inf. % % EXAMPLES: % % % choose source data: % % ----------------- 1 ----------------- % x = -2:2; % y = [ 4 1 0 -1 -4 ]; % % ----------------- 2 ----------------- % x = [-2,-1,1,2]; % y = [ 4 1 1 4]; % % ----------------- 3 ----------------- % x = unique(rand(21,1)*20-10); % y = rand(size(x))-0.5; % % ------------------------------------- % % % calculate spline % pp = spline(x,y); % xo = rootspline(pp); % % plot % xx = linspace(min(x),max(x),50*numel(x)); % figure % plot(xx,ppval(pp,xx)) % hold on, grid on % plot(x,y,'k*') % plot(xo,0*xo,'ro') % legend({'spline','source data','roots'},'location','best') % xlabel('x'), ylabel('y'), title('rootspline example') % % Compare examples 1), 2) and note the roundoff susceptibility! % % % See also % roots, spline, interp1, pchip, mkpp x = pp.breaks(:); y = [pp.coefs(:,end);ppval(pp,x(end))]; % preallocate zo = zeros(pp.pieces,1); % initialize roots counter cnt = 0; % step through all pieces of the spline and find roots for idx = 1:pp.pieces r = rootsbnd(pp.coefs(idx,:),pp.breaks([idx,idx+1])); zo(cnt+1:cnt+numel(r)) = r; cnt = cnt + numel(r); end zo(cnt+1:end) = []; % get unique roots (zo is already sorted) if numel(zo)>1 zo = zo([true;diff(zo)~=0]); end end function r = rootsbnd(p,bnd) % rootsbnd find roots of a polynomial piece inside an interval % % r = rootsbnd(p,bnd) returns the roots of the polynomial represented by % the coefficient vector p (e.g. p = pp.coefs, pp produced by spline) % inside the interval bnd = [x1,x2]. r = roots(p); % sort out non-real and outlying roots idc = imag(r)==0; idc = idc & (r>=0 & r<=diff(bnd)); r = r(idc); % sort multiple roots if numel(r)>1 r = sort(r); end % shift to interval r = r + bnd(1); end <<
From: John D'Errico on 22 Apr 2010 09:27 "Andres" <rantore(a)werb.deNoRs> wrote in message <hqpgk4$l4q$1(a)fred.mathworks.com>... > Hi, > > I'd like to speed up my code for the following problem: > > input: two vectors x, y (numel(x) = numel(y) = ~5e5) > output: roots of the cubic spline through (x,y) > > I've written a function 'rootspline' to get the roots of the piecewise polynomial form of the cubic spline 'pp' (see below), so I run > > pp = spline(x,y); > xo = rootspline(pp); Some years ago, I wrote an inverse tool for cubic splines. Rather than using roots, use the explicit formula for the roots of a cubic polynomial, given the coefficients. This can be done in a fully vectorized form, so all the intervals can be processed at once. Other tricks you can do are to break up knot intervals into intervals such that the spline is strictly monotone over each interval. That will allow you to only apply a roots operation to the intervals where you KNOW a root exists already. John
|
Pages: 1 Prev: Normal direction of faceNormals Next: Wavelet toolbox and linear system |