MathGroup Archive 2005

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

Search the Archive

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
    
















  • Prev by Date: Re: Table function
  • Next by Date: Re: Question to the world: Do you want an computer algebra program with a true GUI?
  • Previous by thread: Re: Google's aptitude test equation
  • Next by thread: Re: Question to the world: Do you want an computer algebra program with a true GUI?