[Date Index] [Thread Index] [Author Index]

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
>send me your opinion.
>
>               ---------------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 
advantage of this information.

>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

at http://groups-beta.google.com/group/comp.soft-sys.math.mathematica.

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/

____________________________________________________________________