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.

Output from example problem

`pp1` and `pp1` are structures containing fields `breaks` and `coefficients`, that define the polynomial pieces. So you can use `mkpp` to build a new piecewise polynomial combining the pieces from `pp1` and `pp2`. You just need to pass the sorted breaks and coefficients to `mkpp`. The only cumbersome part is building the composite coefficient matrix with an order equal to the maximum of the orders of `pp1` and `pp2`, sorting the breaks (removing duplicates), and applying that same sorting to the rows of the composite coefficient matrix. — Luis Mendo, Jun 23 '15 at 15:47
@LuisMendo as you say it is cumbersome to sort and build the new piecewise polynomial. That's why I'm wondering if there is an easy way to do this, e.g. with a built-in function or some small combination of functions. — Delyle, Jun 23 '15 at 15:54
@LuisMendo (Unless I'm missing something) The work can be greatly reduced if both `pp1` and `pp2` use the same number of control points, either 5 or 7 or whatever. Then isn't the correct resultant coefficient matrix just the addition of the other 2 with `breaks = unique([pp1.breaks,pp2.breaks])`? — TroyHaskin, Jun 23 '15 at 16:02
Ah, I missed the arbitrary part and assumed the MWE was the use case. Well then, very interesting indeed! — TroyHaskin, Jun 23 '15 at 16:06
@TroyHaskin If I understand correctly what you say: Yes, if the breaks are the same it's just a matter of adding the coefficients (and then `breaks = pp1.breaks`, no need for `unique `). in the general case, it's an interesting problem! — Luis Mendo, Jun 23 '15 at 16:12
@TroyHaskin Yes, the problem is almost trivial if the same control points are used between `pp1` and `pp2`. However, in my case I use arbitrary control points that do not always line up, as in the example above. So I'm wondering if there is an easy(ish) way to generate `pp3` when `pp1` and `pp2` do not have the same control points. — Delyle, Jun 23 '15 at 16:18
Can you make do with an anonymous function? `pp12 = @(x) ppval(pp1,x)+ppval(pp2,x);`. Then you would call it normally: `plot(tnew, pp12(tnew))` — Luis Mendo, Jun 23 '15 at 16:18
@TroyHaskin Your implementation with the anonymous function is nicer but I would still like the piecewise-polynomial form of pp12. — Delyle, Jun 23 '15 at 16:38

Delyle · Accepted Answer · 2015-06-25T00:48:01.043

3

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:

Piecewise polynomials

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.

edited Jun 25 '15 at 00:48

answered Jun 23 '15 at 19:33

Delyle

529
3
14

This does give _exact_ results. The problem is not exactness, but rather _type_: you get an anonymous function, which is a different "object" than a piecewise polynomial (which is actually a struct with certain fields) – Luis Mendo Jun 24 '15 at 13:09
@LuisMendo I'm not sure I agree. In the solution here, I get a piecewise polynomial (`pp3`, as a struct with fields) that is derived from the anonymous function (more precisely, from the vector produced from evaluating the anonymous function `pp12` at the breakpoints). When I say `pp3` is not _exactly_ `pp1`+`pp2`, what I mean is that the coefficients in `pp3` are not what you would get if you did all the algebra to expand and collect the terms in `pp1` and `pp2` within each segment. My claim is validated by finding `max(abs(ppval(pp3,tnew) - pp12(tnew))) ~= 0`, is it not? – Delyle Jun 24 '15 at 18:36
2

But those are just _numerical errors_ because the computation is done differently in either case. Note that the difference is of the order of `eps`. For example, `sqrt(1e7)^2` and `1e7` _are_ equal, but try `sqrt(1e7)^2-1e7` and you won't get exactly zero. – Luis Mendo Jun 24 '15 at 22:42
@LuisMendo I see. You've convinced me. I'll amend the answer. – Delyle Jun 25 '15 at 00:45

score 0 · Answer 2 · edited Sep 02 '19 at 07:31

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 ."

Adding piecewise polynomials in MATLAB

2 Answers2