NDSolve, InterpolatingFunction objects and Splines
- To: mathgroup at smc.vnet.net
- Subject: [mg54983] NDSolve, InterpolatingFunction objects and Splines
- From: zosi <zosi at to.infn.it>
- Date: Wed, 9 Mar 2005 06:34:12 -0500 (EST)
- Organization: INFN sezione di Torino
- Sender: owner-wri-mathgroup at wolfram.com
Gentlemen,
Question:
How can we extract the **coefficients** of the polynomials
which constitute the InterpolatingFunction Object (IFO) representing
the solution of NDSolve ?
We have been able to extract only the abscissae (x_i) - ordinates {y_i)
of the numerical output; we suspect they are hidden.
Our curiosity arises from the fact that we have tried
(rather naively) to improve the representation of the output
of NDSolve applied, to begin, to a very simple ODE
tfin= 3 Pi;
solnNDS = NDSolve[{x''[t]+x[t]==0, x[0] == 1,x'[0]==1}, x, {t, 0, tfin }];
solnum[t_] := x[t] /. First[solnNDS];
To begin, we have extracted the couples
nodes = (x /. First[solnNDS])[[3, 1]];
solnum[nodes];
Length[nodes];
then
residualfunction[t_] := Evaluate[x''[t] + x[t] /. First[solnNDS]];
plotresidualfunction=Plot[residualfunction[t], {t, 0, tfin },ImageSize
-> 500];
the residual ranges approx. between +/- 1.5*10^(-6), if we exclude
the very first points.
To have a comparison with Splines, we have
a) applied the command SplineFit (cubic)
to the same set of data (nodes and solnum[nodes]),
b) extracted the polynomials (in each interval)
c) calculated again the residual;
To our surprise the residual ranges between +/- 1*10^(-3).
A loss by a factor of 1000.
In this mail we do not shown the details to calculate the 1st and 2nd
derivative of the spline (see Note 1)
We said "surprise" because the polynomials contained in the IFO are
not so smooth as the splines.
About the behaviour of the IFO, a very simple example by Stan Wagon
(and slightly modified) is illuminating
data = {{-1, -1}, {-0.96, -0.4}, {-0.88, 0.3},
{-0.62, 0.78}, {0.13, 0.91}, {1, 1},
{1.2, 1.5}, {1.5, 1.6}, {2.2, 1.2}};
f = Interpolation[data];
Plot[f[t], {t, -1, 2.2},
Epilog -> {PointSize[0.025], Point /@ data}];
The two bumps clearly indicate that the first derivative
is not continuous in at least two-three points.
Needs["NumericalMath`SplineFit`"];
cub = SplineFit[data, Cubic]; (* see note 1 *)
ParametricPlot[cub[t], {t, 0, 8},
Compiled -> False, PlotRange->All, AspectRatio ->
Automatic,
Epilog -> {PointSize[0.02], Map[Point, data]}];
Conclusion: If we are using the same WorkingPrecision (=Default),
shouldn't be **better** the residual obtained
through the splines ?
Note 1. As we have not been able to calculate directly the 1st and 2nd
derivative
of cub, we have extracted the coefficients of each polynomial
contained in cub,
then made the calculation of the derivative. Is there any
quickier method
to calculate directly the **second** derivative of the
SplineFit output ?
Thanks for your patience and help.
Gianfranco Zosi
Dip Fisica Generale "A. Avogadro"
Universita di Torino
v. P. Giuria 1 - 10125 Torino - Italy