Adding piecewise polynomials in MATLAB

Question

I need to add piecewise polynomials derived from multiple datasets. Is there an easy way to add piecewise polynomials together without interpolating? In other words, given PP1 and PP2, is there a way to generate PP3 (where PP3 remains in a piecewise polynomial form)? e.g...

    t1 = linspace(0,1,5);
    t2 = linspace(0,1,7);
    pp1 = spline(t1,sin(pi*t1));
    pp2 = spline(t2,t2.^2);
    close all
    hold on
    tnew = linspace(0,1,50);
    h(:,1) = plot(tnew,ppval(pp1,tnew));
    plot(t1,ppval(pp1,t1),'bs')
    h(:,2) = plot(tnew,ppval(pp2,tnew));
    plot(t2,ppval(pp2,t2),'rs')
    h(:,3) = plot(tnew,ppval(pp1,tnew)+ppval(pp2,tnew));
    legend(h,{'spline of sin(\pi t)','spline of t^2','sin(\pi t)+t^2'},...
                'location','northwest')
    xlabel('t')

But instead of specifying tnew explicitly, I would like a new piecewise polynomial pp3 that is effectively pp1+pp2.

Answer 1

This is probably the easiest way to get pp1 + pp2 Adding to the code in the question:

    pp12 = @(x) ppval(pp1,x)+ppval(pp2,x);

    breaks = unique([pp1.breaks,pp2.breaks]);
    pp3 = spline(breaks,pp12(breaks));

    plot(tnew,ppval(pp3,tnew),'k:');

Gives the dotted black line:

with pp3 being in piecewise-polynomial form with segments defined only by the breakpoints of pp1 and pp2. Running max(abs(ppval(pp3,tnew) - pp12(tnew))) yields 2.7756e-16, which is of the order of eps.

Thanks to @LuisMendo and @TroyHaskin for their suggestions.

Answer 2

x1  = pp1.breaks;
x2  = pp2.breaks;
x   = unique([x1,x2]);

idx1 = lookup(x1,x(1:end-1));
idx2 = lookup(x2,x(1:end-1));

coefs = pp1.coefs(idx1,:) + pp2.coefs(idx2,:);

pp3 = mkpp(x,coefs);
y   = ppval(pp3,x);
h(:,4)  = plot(tnew,ppval(pp3,tnew));
legend(h,{'spline of sin(\pi t)','spline of t^2','sin(\pi t)+t^2','pp3'},'location','northwest')
xlabel('t')

I thought this code should work, but was not true. I think because the mkpp resultant piecewise polynomial coefficients are local, adding them in the same interval can not solve the problem. But it can be a start point for a solution.

from Matlab help of mkpp: "Since the polynomial coefficients in coefs are local coefficients for each interval, you must subtract the lower endpoint of the corresponding knot interval to use the coefficients in a conventional polynomial equation. In other words, for the coefficients [a,b,c,d] on the interval [x1,x2], the corresponding polynomial is

f(x)=a(x−x1)3+b(x−x1)2+c(x−x1)+d ."