필터 지우기
필터 지우기

How to speed up code with ppval and integral

조회 수: 1 (최근 30일)
Martin Andersson
Martin Andersson 2017년 12월 5일
답변: David Goodmanson 2017년 12월 6일
Hi I want to speed up my code as much as possible. Everything that could help me reduce the time to run the code is appreciated. A simplification of my code is given below and it takes on my computer 120 seconds to run which is a little too long.
data=[0 1 2 3 4 5];% X data for first datapoints
X=repmat(Xdata,44000,1); % make all data points from first example
Ydata=[10 8.3 5.2 4.2 2.9 1.1];% X data for first datapoints
Y=repmat(Ydata,44000,1); % make all data points from first example
Xqdata=[2 2.1 2.3 2.6 3.0 3.6 3.9 4.0 4.2 4.4 4.8];% X data for second datapoints
xq=repmat(Xqdata,44000,1);% make all data points from second example
for i=1:44000 % loop over all datapoints
yi=pchip(X(i,:),Y(i,:), xq(i,:)); % make pchip over first data points and generate yi points
pp=pchip(xq(i,:),yi); % creates a piecewise polynomial structure for use with ppval
Result = integral(@(varible)ppval(pp,varible),xq(i,1),xq(i,end)); % integrates over xq (this is the time consuming part)
end
Is there any way to speed up the integration step?

답변 (1개)

David Goodmanson
David Goodmanson 2017년 12월 6일
Hi Martin,
Here's some code to integrate a piecewise polynomial on the basis of its polynomial coefficients, no actual numerical integration. On my pc with the 44000 repetitions, replacing integral(...) with intpp(pp) cut down the execution time by about a factor of three.
function defint = intpp(pp)
% definite integral of a piecwise polynomial
% c has highest power coefficient in first col, const coefficient in last col
brk = pp.breaks;
c = pp.coefs;
% [brk c] = unmkpp(pp); % still works, one less line of code
dbrk = diff(brk);
dbrk = dbrk(:);
[nrow ncol] = size(c);
intseg = zeros(nrow,1);
for k = 1:ncol
c(:,k) = c(:,k)/(ncol+1-k);
intseg = (intseg + c(:,k)).*dbrk;
end
defint = sum(intseg); % sum of integrations over the segments

카테고리

Help CenterFile Exchange에서 Polynomials에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by