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