Re: Problems with NDSolve and Splines

• To: mathgroup at smc.vnet.net
• Subject: [mg55418] Re: Problems with NDSolve and Splines
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Wed, 23 Mar 2005 05:34:31 -0500 (EST)
• References: <423816FD.1060804@to.infn.it>
• Sender: owner-wri-mathgroup at wolfram.com

```Dear Gianfranco:

>Recently (9 March), I have posted the mail mg54983, Wed 12:34 PM
>about NDSolve, InterpolatingFunction objects and Splines.
>
>Unfortunately, no one has replied,
>notwithstanding the detailed steps. Does it mean
>that my question was trivial ?
>
>As I have seen that in the past you have given (to the mathgroup)
>useful and clear hints, I would be grateful to you if you could
>
>               ---------------mg54983-------------
>
>Gentlemen,
>
>Question:
>How can we extract the **coefficients** of the polynomials
>which constitute the InterpolatingFunction Object (IFO) representing
>the solution of NDSolve ?

InterpolatingFunction uses divided differences to
construct Lagrange or Hermite interpolating
polynomials (Newton interpolation), see

http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html

As I understand it, Newton interpolation is used
because it is simpler to implement and more
efficient than using spline functions.

In 1994 Hon-Wa Tam, formerly from WRI sent me the
following notes on the structure of 1-D
InterpolatingFunction (but note that
InterpolatingFunction was modified in version 3):

The following discussion pertains only to 1-D interpolating functions.
See the description on interp_multidimen() for multi-dimensional
interpolating functions.

Interpolation uses Newton polynomials to construct an interpolation
table.  The object InterpolatingFunction[ range, table ] consists of
the range of the interval and the table.  Each element of table
is responsible for a subinterval (t_{n-1}, t_n].  A table element
has the structure

{ t_n, t_{n-1}, divided differences, coefficients }.

The order of the particular polynomial for (t_{n-1}, t_n] is equal to
Length( coefficients ).

Interpolation within (t_{n-1}, t_n] is done by

dvdf[1] + ( x - (t_n-coeff[1]) ) * dvdf[2] +
( x - (t_n-coeff[1]) ) * ( x - (t_n-coeff[2]) ) * dvdf[3] + ...

Each table has a degenerate element whose subinterval is of width 0.
This is for the purpose of finding the right subinterval by
binary search.

InterpolatingFunction objects can be constructed by Interpolation[]
or by NDSolve[].

If an InterpolatingFunction object is constructed by Interpolation[],
we have t_{n-1} < t_n in each table element.  This is also true
for NDSolve[] when the given initial point is at the left of the
domain, as in:

NDSolve[ { y'[x] == y[x], y[0] == 1 }, y, { x, 0, 1 } ].

In this normal case, the degenerate subinterval is on the L.H.S. and
is the first element of the interpolation table.

However, NDSolve can also integrate in the negative direction:

NDSolve[ { y'[x] == y[x], y[0] == 1 }, y, { x, 0, -1 } ].

In this case all subintervals are processed "backward" and we have
t_{n-1}  t_n.  The degenerate subinterval is now on the R.H.S. and
is the last element of the interpolation table.

It is also possible that some subintervals are in the +ve, and some
in the -ve direction:

NDSolve[ { y'[x] == y[x], y[0] == 1 }, y, { x, -1, 1 } ].

Here the degenerate subinterval is somewhere in the middle of the
interpolation table.

>We have been able to extract only the abscissae (x_i) - ordinates {y_i)
>of the numerical output; we suspect they are hidden.

You can extract other parts of the Interpolation
function. The fourth part is the one you are
after.

>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; see Note 1.

To improve the representation, have you tried
changing NDSolve options? For example, you could
increase the WorkingPrecision.

>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 show the details to calculate the 1st and 2nd
>derivative of the spline (see Note 2).
>
>We said "surprise" because the polynomials contained in the IFO are
>not so smooth as the splines.

I think that you are missing something here. If
you look up Interpolation in the HelpBrowser you
will see that:

Data can be given in the form {..., {xi, {fi, dfi, ddfi, ...}}, ...} to
specify derivatives as well as values of the function at the points xi. You
can specify different numbers of derivatives at different points.

Now, since NDSolve is solving the differential
equation, you should expect that it would take

>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.

Correct -- but if you had supplied the required
information you could have produced a smoother
InterpolatingFunction object.

>    Needs["NumericalMath`SplineFit`"];
>    cub = SplineFit[data, Cubic];    (* see note 2 *)
>    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 ?

No.

>Note 1. The main reason is that we have to solve
>in cascade three second order hyperbolic PDE
>(Takagi equations), i.e., the output from the
>first PDE becomes the input to the second PDE
>and so on. The output from the first PDE is
>**highly** oscillating and, even though  we have
>anyway obtained reasonable results by exploiting
>directly the IFO, why couldn't we exploit the
>Splines to represent the output from the PDE?
>Splines garantee a smoother behaviour.

I do not think that this is correct.
Nevertheless, it is obviously important for your
cascade to have reliable output from one stage
feeding into the next. I would have a close read
of the advanced documentation for NDSolve in the
HelpBrowser.

>Note 2. 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 ?

There have been several MathGroup postings on this topic. Do a search on

derivative splinefit

Cheers,
Paul

--
____________________________________________________________________
Paul Abbott                                   Phone: +61 8 6488 2734
School of Physics, M013                         Fax: +61 8 6488 1014
The University of Western Australia      (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009                      mailto:paul at physics.uwa.edu.au
AUSTRALIA                            http://physics.uwa.edu.au/~paul

Conference Chair for International Mathematica Symposium 2005
http://InternationalMathematicaSymposium.org/IMS2005/

____________________________________________________________________

```

• Prev by Date: Re: Q: Varying colors in Plot3D according to Height
• Next by Date: point in convex hull
• Previous by thread: Re: using Mathematica to plot in visual studio 2005
• Next by thread: point in convex hull