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/ ____________________________________________________________________