# Merging piecewise polynomial (pp) structures

10 views (last 30 days)

Show older comments

##### 0 Comments

### Answers (3)

Bruno Luong
on 28 Mar 2023

Edited: Bruno Luong
on 28 Mar 2023

Using BSFK from this FEX https://www.mathworks.com/matlabcentral/answers/1810010-free-knot-spline-approximation-bsfk-problem?s_tid=ans_lp_feed_leaf

% Example as in John's answer

breaks = [4 10]; coefs = [ 0 0 1 -2 53];

PP1 = mkpp(breaks,coefs);

breaks = [8 15]; coefs = [ -1 6 1 4 77];

PP2 = mkpp(breaks,coefs);

% overlap interval

oI = [8 10];

li = diff(oI); % its length

f = @(x) ppval(PP1,x);

g = @(x) ppval(PP2,x);

% define smooth transition weight function

sigfun = @(x) x.^2.*(3-2*x);

w2 = @(x) sigfun((x-oI(1))/li);

w1 = @(x) 1-w2(x);

% weighted sum of f and g, applied only in oI

h = @(x) w1(x).*f(x) + w2(x).*g(x);

%% Fit on 3 intervals

arg = {5,[],[],struct('KnotRemoval', 'yes','Animation',0)}; % EDIT first argument 4 => 5

xo = linspace(oI(1),oI(2),100);

ppo = BSFK(xo,h(xo),arg{:});

x1 = linspace(PP1.breaks(1),oI(1));

ppf = BSFK(x1,f(x1),arg{:}); % we can also use PP1 instead

x2 = linspace(oI(end),PP2.breaks(end));

ppg = BSFK(x2,g(x2),arg{:});

% Construct ppmerge structure

ppmerge = ppf;

ppmerge.breaks = [ppf.breaks(1:end-1),ppo.breaks(1:end-1),ppg.breaks(1:end)];

ppmerge.coefs = [ppf.coefs;

ppo.coefs;

ppg.coefs];

ppmerge.pieces = size(ppmerge.coefs,1);

% Graphical check

figure

hold on

xf = linspace(PP1.breaks(1),PP1.breaks(end));

xg = linspace(PP2.breaks(1),PP2.breaks(end));

plot(xf,f(xf),'b+');

plot(xg,g(xg),'rx');

xm = linspace(ppmerge.breaks(1),ppmerge.breaks(end));

plot(xm,ppval(ppmerge,xm),'g','LineWidth',2);

legend('f','g','merge')

A closer look in the overlap interval, it seems the transition is smooth as expected

##### 11 Comments

Bruno Luong
on 30 Mar 2023

Edited: Bruno Luong
on 30 Mar 2023

What I means in "subdivisions of f/g subintervals inside overlap region" is that

in the region of [8,10] in the example is the "overlap region".

In this region f has multiple break points, they define a set of "subintervals of f inside overlap region". Each of the subinterval of f is a polynomial of order k-1, since f is pp. So if you select k sample points in each sub intervals (eg. equidistance including the two end break points, thus the word "subdivision") then f is fully characterized by the values computed at this subdivision samples.

That define the sample points for f. Then do the same procedure for g. Take the union of two sets of sampling points, this gives us a way to select the points in the overlap region, and I'm sure it is large enough to recover accurately the approximation of

h = w1*f + w2*g.

on the overlap region by fitting. It is better than selecting emprirical 100 points as I did in my code and it can capture the sharp corner where it needs, etc...

Matthieu
on 24 Mar 2023

Are you trying to 'concatenate' your polynomial structures ?

I am not an expert, although I know those are defined by 1 vector and a matrix, breaks & coefs. break defines the range of action of the coefficients of your polynomial, stored in coefs.

You can get those vectors with:

breaks = [0 4 10 15];

coefs = [0 1 -1 1 1; 0 0 1 -2 53; -1 6 1 4 77];

pp = mkpp(breaks,coefs) ;

pp.breaks

pp.coefs

A way of merging those would then be to concatenate the breaks and coefs of all piecewise polynomials you want to merge, as below :

breaks2 = [15 20 22 28];

coefs2 = [0 0 1 2 3; 0 0 0 0 5; -0.1 1 0 0 8] ;

breaks_merged = horzcat(breaks,breaks2(2:end));

coefs_merged = vertcat(coefs,coefs2);

pp_merged = mkpp(breaks_merged,coefs_merged) ;

t = -5:0.1:33 ;

plot(t,ppval(pp_merged,t))

Is this what you needed ?

##### 4 Comments

John D'Errico
on 24 Mar 2023

Note that a weighted sum will implicitly increase the order of the polynomial segment in the merged region.

A linear weight will result in a function that is one degree higher in the overlapped interval, but it will also result in a derivative discontinuity at the break. So continuous, but not differentiable. So it should also be possible to employ a nonlinear (actually just a higher polynomial order) weighting scheme, if it were important to acheive differentiability on the overlapped interval too.

John D'Errico
on 24 Mar 2023

Edited: John D'Errico
on 24 Mar 2023

Is there any method to "merge" two PP functions? Well, no. Should there be one? I've been using and writing splines tools for 40+ years, in MATLAB for 35 years or so, and answering questions about them for as long. And I've never seen that request, certainly not as you are asking. You don't write code to do something that nobody will ever want to use.

Could there be a way to do so,? Well, yes. The obvious one is to just expand the list of breaks, then append the pp segments as suggested already.

However, your question is even less, let me say, "expected", because you have an overlap. What would you do in the overlap region? Arbitrarily choose one function of the other? Take the mean? Flip a coin to decide?

Sigh. COULD you do something? Well, I suppose you could do something. I could think of several choices even here. The problem is in the overlap region, you need to resolve the conflct. That is, suppose we have intervals at the ends of domain 1 as [A,B], and at the beginning of domain 2 as [C,D]. I'll assume from what you have said is that A < C < B < D. Otherwise there is no problem. Actually, I can see another issue. What if B < C? That is, the problem is simple, if B == C. Now the direct merger already suggested works. But if B < C, then there is an undefined region. I suppose you could extrapolate the curves, over the hole in the middle between domains. But extrapolating splines is about the most obscene thing you can do. Avoid it at all costs.

Anyway, suppose we have the sceneario where A < C < B < D? We could now adjust the breaks to be exactly that sequence. Where on the interval [A,B), we use PP1 from the lower domain. On the interval [B,D), we use PP2, so the upper segment.

But on the interval [C,B) What choice do we have? Perhaps we might just use the average of the two functions on that interval. But that would now almost certainly introduce a point of discontinuity at the breaks, a bad thing.

One idea might be to use an average that varies along the interval. Essentially, this would force the curve to be a higher order segment. So instead of a cubic segment there, it would now be, a degree 4 segment at least. That would allow the resulting function to be continuous, though not necessarily differentiable. If we are willing to make it a 5th degree segment, that should allow one to construct a "merged" segment that would be both continous AND differentiable. (Just thinking off the cuff there.)

So, doable. But, for example, what would you do here?

PP1 = mkpp([0 2],[1 0]);

PP2 = mkpp([1 3],[-1 2]);

fnplt(PP1)

hold on

fnplt(PP2)

hold off

Lets see. The new pp form would have 4 breaks.

allbreaks = sort([PP1.breaks,PP2.breaks]);

S1 = [0,PP1.coefs];

S3 = [0,PP2.coefs(1),fnval(PP2,allbreaks(3))];

S2a = [PP1.coefs(1),fnval(PP1,allbreaks(2))];

S2b = PP2.coefs;

S2 = conv([-1,1],S2a) + conv([1 0],S2b);

PPmerged = mkpp(allbreaks,[S1;S2;S3]);

fnplt(PPmerged,'g')

hold on

fnplt(PP1,'r')

fnplt(PP2,'b')

hold off

So the "merged function is the same as the old ones, in the intervals [0,1], and [2,3]. nd it is made to be continuous at the overlap points, but differentiability was not enforced. On the sub-interval [1 2], the segment is now quadratic, so one degree higher order. That was necessary to enforce continuity.

As I said, doable.

##### 26 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!