BCs with few knots incorrect

[asparc]

Hi!

Thanks for making this very useful tool. It's very easy to use. However, I think I stumbled upon some unexpected behaviour.

I wanted to make a spline with boundary conditions on the second derivative. So I expect the spline to evaluate to zero at both boundaries, like so:

julia> B1 = RecombinedBSplineBasis(Derivative(2), BSplineBasis(4, [0,0.25,0.5,0.75,1]));
julia> C1 = collocation_matrix(B1, [0,1], Derivative(2), Matrix{Float64});
julia> coefs1 = [1,1,1,1,1];
julia> println("2nd derivative at boundaries: ", C1*coefs1)
2nd derivative at boundaries: [-1.4210854715202004e-14, -1.4210854715202004e-14]  

So far, this is expected behaviour. However, when I reduce the number of knots from 5 to 3, which should still be valid, I get non-zero values for the second derivative at the right boundary:

julia> B2 = RecombinedBSplineBasis(Derivative(2), BSplineBasis(4, [0,0.5,1]));
julia> C2 = collocation_matrix(B2, [0,1], Derivative(2), Matrix{Float64});
julia> coefs2 = [1,1,1];
julia> println("2nd derivative at boundaries: ", C2*coefs2)
2nd derivative at boundaries: [-3.552713678800501e-15, 17.999999999999996]

Could this be a bug?

[jipolanco]

Hi, I'm glad you find this package useful!

This definitely looks like a bug, thanks for reporting it. I'll try to get it fixed when I find some time within the next few days.

[asparc]

Awesome! Thanks for the very quick reply!

[jipolanco]

After thinking about this issue, I believe this is a limitation of the recombination procedure, which defines basis functions as linear combinations of B-splines near the boundaries. This can be represented using a "recombination" matrix.

In your first example, the recombination matrix looks as follows:

julia> recombination_matrix(RecombinedBSplineBasis(Derivative(2), BSplineBasis(4, [0,0.25,0.5,0.75,1])))
7×5 RecombineMatrix{Float64, Tuple{Derivative{2}}, 2, 3, StaticArraysCore.SMatrix{3, 2, Float64, 6}}:
 1.2   ⋅    ⋅    ⋅    ⋅
 0.8  0.5   ⋅    ⋅    ⋅
  ⋅   1.5   ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅   1.5   ⋅
  ⋅    ⋅    ⋅   0.5  0.8
  ⋅    ⋅    ⋅    ⋅   1.2

which basically means that the first (last) two functions of the recombined basis are obtained as recombinations of the first (last) three functions of the original B-spline basis. In this case, the original basis has 7 functions and the recombined one has 5. Note that we need to recombine a total of 6 different B-spline functions (B[1], B[2], B[3], B[5], B[6], B[7]).

When you reduce the size of the original basis to 5, as in your second example, it is of course no longer possible to recombine 6 different B-spline functions. Right now the recombination matrix gives this:

julia> recombination_matrix(RecombinedBSplineBasis(Derivative(2), BSplineBasis(4, [0,0.5,1])))
5×3 RecombineMatrix{Float64, Tuple{Derivative{2}}, 2, 3, StaticArraysCore.SMatrix{3, 2, Float64, 6}}:
 1.2   ⋅    ⋅
 0.8  0.5   ⋅
  ⋅   1.5   ⋅
  ⋅    ⋅   0.8
  ⋅    ⋅   1.2

which is wrong. I think the only solution here is to just throw an error when the original basis is not large enough for the wanted recombination.

[asparc]

I agree that you can't treat the left/right boundary conditions independently any more; their recombinations will overlap. However, I'm still convinced that you have enough degrees of freedom to find a recombination that satisfies the Derivative(2) condition with a 2-segment cubic spline. I'd have to sit down and write it out at some point.

[jipolanco]

I'll have to give it some more thought as well. I'll be happy to know if you come up with something.