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 Internal`PartitionRagged
.
Examples
OP's.
ifn = NDSolveValue[{y''[x] + y[x]^3 == 0, y[0] == 1, y'[0] == 0},
y, {x, 0, 10}, Method -> "Extrapolation"];
Show[
ListLinePlot[ifn, PlotStyle -> AbsoluteThickness[4.6]],
Graphics[{
ColorData[97][2], AbsoluteThickness[2.6],
ifnToBezierCurve[ifn]
}]
]
Developer`PackedArrayForm
of order 3
ifnToBezierCurve[
NDSolveValue[{y''[x] + y[x]^3 == 0, y[0] == 1, y'[0] == 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[0] == 1, y'[0] == 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[0] == 1, y'[0] == 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[0] == 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[97][2], 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_, {Developer`PackedArrayForm, split_, a_}, ___],
Automatic] /; if["InterpolationMethod"] === "Hermite"(*/;BitAnd[
type[[$bitFlagsPos]],2^$rectArrayBit]>0*):=
<|"Coordinates" -> coords,
"FunctionData" ->
Internal`PartitionRagged[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[1], 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,
Internal`PartitionRagged[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]},
MapThread[
Function[{t, x},(* time steps, function step data *)
If[
MatrixQ[x] &&
2 Length@x[[1]] > First@ifn@"InterpolationOrder",
(* all components the same sufficient degree *)
With[{sd = 2 Length@x[[1]] - 1},(* SplineDegree *)
BezierCurve[ (* composite BezierCurve *)
Join[
{{t[[1]], x[[1, 1]]}},
Flatten[
MapThread[
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]
];