# Convert interpolating function to a Bezier curve

I've been wondering, is there a way to convert an interpolation into a sequence of BezierCurve objects that form each step in the interpolation?

There is an easy way to plot an InterpolatingFunction as shown in Easy way to plot ODE solutions from NDSolve?:

ifn = NDSolveValue[{y''[x] + y[x]^3 == 0, y == 1, y' == 0},
y, {x, 0, 10}, Method -> "Extrapolation"];
ifn // ListLinePlot But as can be seen, when a method takes large steps, the path plotted is not smooth. This can be easily fixed with ListLinePlot[ifn, InterpolationOrder -> 3] or with Plot[ifn[t], {t,...}]. So my question is perhaps academic. However, both methods end up sampling the function at points between the steps. It is a trivial objection, I know, since I'm pretty sure any curve (e.g. Circle[] as well as BezierCurve[]) is produced by connecting sample points.

The function hermiteToBezierPoints[] computes the Bezier control points for BezierCurve from values of the function and its derivatives up to some order at each endpoint of an interval. The function ifnToBezierCurve[ifn] passes the relevant data for an interpolating function ifn to hermiteToBezierPoints[] and wraps the results in BezierCurve or JoinedCurve[BezierCurve[...],...]. I originally started down this line of development when I read this Q&A: Plot with NDSolve for a range of initial values. So the motivation was to plot NDSolve solutions via Bezier curves.

The function ifnToBezierCurve[] might convert any univariate scalar-valued InterpolatingFunction, although it complains about those that do not have "Hermite" as the "InterpolationMethod". I had to guess what InterpolatingFunction does when the function data at the end points is insufficient to determine an interpolant of degree InterpolationOrder. These are the ones that can have kinks. It uses data from adjacent points, but I had to guess exactly how. So far it has worked on the random examples I threw at it. It handles ones that have jump discontinuities with a simple application of Split[..., Unequal] and InternalPartitionRagged.

### Examples

OP's.

ifn = NDSolveValue[{y''[x] + y[x]^3 == 0, y == 1, y' == 0},
y, {x, 0, 10}, Method -> "Extrapolation"];
Show[
ListLinePlot[ifn, PlotStyle -> AbsoluteThickness[4.6]],
Graphics[{
ColorData, AbsoluteThickness[2.6],
ifnToBezierCurve[ifn]
}]
] DeveloperPackedArrayForm of order 3

ifnToBezierCurve[
NDSolveValue[{y''[x] + y[x]^3 == 0, y == 1, y' == 0},
y, {x, 0, 1}]] // Graphics Variable order, non-Hermite method. Graphics output same as above.

ifnToBezierCurve[
NDSolveValue[{y''[x] + y[x]^3 == 0, y == 1, y' == 0},
y, {x, 0, 1}, InterpolationOrder -> All]] // Graphics


ifnToBezierCurve::meth: Currently converting interpolations using interpolation method Local Taylor series are unimplemented; returning a BezierCurve based on Hermite method with SplineDegree -> 3.

You can override the warning by explicitly setting SplineDegree, which then creates Bezier curves of uniform degree from the interpolating function and its derivatives. (Graphics same as above.)

ifnToBezierCurve[
NDSolveValue[{y''[x] + y[x]^3 == 0, y == 1, y' == 0},
y, {x, 0, 1}, InterpolationOrder -> All],
SplineDegree -> 5] // Graphics


Discontinuities

Here's an ODE that has jumps:

Graphics[{
ifnToBezierCurve[
NDSolveValue[{y'[x] + y[x] ==
2 DiracDelta[x - 1] + 2 DiracDelta[x - 2] + 2 DiracDelta[x - 3],
y == 1}, y, {x, 0, 4}]]
},
PlotRange -> All,
Options@Plot] Irregular, user-defined interpolation with kinks

idata = {{{0}, 1, 2, -1}, {{1}, -1}, {{3/2}, 1}, {{2}, 3}, {{5/2},
0}, {{3}, -1, 1, 2, -3, 1}};
foo = Interpolation[N@idata, InterpolationOrder -> 4];
Show[
Plot[foo[t], {t, 0, 3}, PlotStyle -> AbsoluteThickness[5.6],
PlotRange -> All],
Graphics[{ColorData, AbsoluteThickness[2.6],
ifnToBezierCurve@foo
}]
] ### Code dump

(* Bezier pts from derivatives at x0, x1 *)
hermiteToBezierPoints // ClearAll;
hermiteToBezierPoints[{{x0_}, f0__}, {{x1_}, f1__}] :=
hermiteToBezierPoints[{x0, x1}, {f0}, {f1}];
hermiteToBezierPoints[{x0_, x1_}, f0_?VectorQ, f1_?VectorQ] :=
Module[
{dx, order0, rng0, order1, rng1, n, yy0, yy1},
dx = x1 - x0;
order0 = Length@f0 - 1;
rng0 = Range[0, order0];
order1 = Length@f1 - 1;
rng1 = Range[0, order1];
n = order0 + order1 + 1;(* Bezier/Hermite order *)
yy0 = f0*dx^rng0/Pochhammer[n + 1 - rng0, rng0];
yy1 = f1*(-dx)^rng1/Pochhammer[n + 1 - rng1, rng1];
Transpose@{
Subdivide[x0, x1, n],
Join[
Outer[Binomial, rng0, rng0] . yy0,
Outer[Binomial, rng1, rng1] . yy1 // Reverse]
}
];

(*"
* Single BezierCurve per segment alternative
"*)

(* unneeded as yet
(*bit field positions-inferred,perhaps mistaken*)
$$bitFlagsPos=2;$$extrapBit=0;(*whether to warn about extrapolation*)
$$rectArrayBit=1; (*whether the data (f,f',...) is a rectangular array \ (not ragged) *)$$machPrecBit=2; (*whether the data (f,f',...) is MachinePrecision *)
$repeatedBit=4; (*whether repeated abscissae are permitted *) *) ifnGetData // ClearAll; iIfnGetData // ClearAll; (*** internal function ***) iIfnGetData[ if : Verbatim[InterpolatingFunction][_, type_, coords_, {DeveloperPackedArrayForm, split_, a_}, ___], Automatic] /; if["InterpolationMethod"] === "Hermite"(*/;BitAnd[ type[[$bitFlagsPos]],2^\$rectArrayBit]>0*):=
<|"Coordinates" -> coords,
"FunctionData" ->
InternalPartitionRagged[a, Differences@split]|>;
iIfnGetData[
if : Verbatim[InterpolatingFunction][_, type_, coords_,
a : {__List}, ___], Automatic] /;
if["InterpolationMethod"] === "Hermite" :=
<|"Coordinates" -> coords,
"FunctionData" -> a|>;
(* non-Hermite *)
iIfnGetData[if_InterpolatingFunction, Automatic] /;
if["InterpolationMethod"] =!= "Hermite" := (
Message[ifnToBezierCurve::meth, if["InterpolationMethod"],
First@if["InterpolationOrder"]];
iIfnGetData[if, First@if["InterpolationOrder"]]);
(* non-Automatic SplineDegree *)
iIfnGetData[if_InterpolatingFunction, sdeg_Integer?Positive] :=
(* SplineDegree deg will be effectively be odd *)
With[{coords = if@"Coordinates",
a = Transpose@
Through[NestList[Derivative, if, Floor[sdeg/2]][
"ValuesOnGrid"]]},
<|"Coordinates" -> coords,
"FunctionData" -> a|>
];
(*** UI function ***)
ifnGetData[if_InterpolatingFunction, sdeg_] :=
Module[{data},
data = iIfnGetData[if, sdeg];
With[{coordsegs = Split[First@data@"Coordinates", Unequal]},
{coordsegs,
InternalPartitionRagged[data@"FunctionData",
Length /@ coordsegs]}
]];

ifnToBezierCurve // ClearAll;
ifnToBezierCurve::meth =
"Currently converting interpolations using interpolation method  \
are unimplemented; returning a BezierCurve based on Hermite method \
with SplineDegree -> .";
(* If SplineDegree is not Automatic then Hermite interpolation model \
is used *)
ifnToBezierCurve // Options = {SplineDegree -> Automatic};
ifnToBezierCurve[ifn_InterpolatingFunction, OptionsPattern[]] /;
Length@ifn@"Domain" == 1 && ifn@"OutputDimensions" === {} :=
With[{data = ifnGetData[ifn, OptionValue@SplineDegree]},
Function[{t, x},(* time steps, function step data *)
If[
MatrixQ[x] &&
2 Length@x[] > First@ifn@"InterpolationOrder",
(* all components the same sufficient degree *)
With[{sd = 2 Length@x[] - 1},(* SplineDegree *)
BezierCurve[ (* composite BezierCurve *)
Join[
{{t[], x[[1, 1]]}},
Flatten[
Rest[hermiteToBezierPoints[##]] &,
{Partition[t, 2, 1], Most[x], Rest[x]}
],
1]
],
SplineDegree -> sd]
],
(* degrees of components vary or
are less than InterpolatingOrder *)
Module[{delayedRest, iorder, ncoeffs, i},
delayedRest := (delayedRest = Rest; Identity);
ncoeffs = Length /@ x;
JoinedCurve@ (* joined BezierCurve pieces *)
Table[
iorder = ncoeffs[[j]] + ncoeffs[[j + 1]] - 1;
i = 0;
While[iorder < First@ifn@"InterpolationOrder",
++i;
With[{dj = 3/4 - 1/4 (-1)^i (3 + 2 i)},
If[1 <= j + dj <= Length@x,
iorder += ncoeffs[[j + dj]]
]]
];
With[{bpts = hermiteToBezierPoints[
{t[[j]], t[[j + 1]]},
x[[j]],
Join[
x[[j + 1]],
Table[Derivative[k][ifn][t[[j + 1]]],
{k, ncoeffs[[j + 1]], iorder - ncoeffs[[j]]}]]
]},
BezierCurve[delayedRest@bpts,
SplineDegree -> Length[bpts] - 1]
],
{j, Length@t - 1}
]]
] (* end If[] *)
],
data
] /; FreeQ[data, ifnGetData]
];
`