Dense Output, Chapter 17, ODE, specific y at x

[dpr]

Methods Forum Members,
I am interested in adapting the example of section 17.0.4 to generate output at values that I specify, while still using the adaptive stepsize algorithm. The discussion of section 17.2.2, Dense Output, suggests that the NR interpolation method, Dormand-Prince, as implemented in StepperDopr5 should be able to provide a solution to this problem. The interface used in the example to generate output is via the Output declaration. Is there a way to specify a vector of particular x values at which y values will be returned through this declaration? Is there an example on how to generate output at particular x values that might be useful?
Thank you

[kutta]

Hello NRMember,
Though this method below will help you to go further in your stated method,I shall be glad if you forward your view on this C progrm.GoodLuck
Thanks
C.R.Muthukumar(kutta)

#include <conio.h>


int main()
{
int i,n,j,k,l;
float xo,y[20],f[10][10],X[10],Y[10],h,u,p;
printf("Enter the value of n(No of data pairs-1):\n");
scanf("%d",&n);
printf("Enter the initial value of x: \n");
scanf("%f",&xo);
printf("Enterthe step size h:\n");
scanf("%f",&h);
printf("Enterr the values of y:\n");
for(i=0;i<n+1;i++)
scanf("%f",&y[i]);
printf("Enter the required no of interpolated values of y : \n");
scanf("%d",&l);
printf("Enter the values %d values of x for which values of y are required: \n",l);
for (k=0;k<l;k++)
scanf("%f",&X[k]);
for(j=0;j<n+1;j++)
f[0][j]=y[j];
for(i=1;i<n+1;i++)
for(j=0;j<n+1-i;j++)

f[i][j]=f[i-l][j+l]-f[i-l][j];
for(k=0;k<l;k++)
{u= ( X[k]-xo)/h;
Y[k]=y[0];
p=1;
for(i=1;i<n+1;i++)
{p=p*(u-i+1)/i;
Y[k]=Y[k]+p*f[i][0];
}
printf("The values of X and Y are : %f\t%\n",X[k],Y[k]);
}
getch();
}

[gjm]

kutta, I don't think your proposal is what the original poster is asking for. Your code takes (say) N equally spaced samples computed by the DE solver, fits the unique polynomial of degree <= N-1 that goes exactly through those points, and uses that for interpolation. (At least, I think that's what it does. I haven't checked in detail; in particular, I haven't looked for bugs.) The original poster is asking for an interpolation method that interpolates on each interval between samples using function and derivative data from the integration process. This is not at all the same thing. Global high-order polynomial fitting is likely to be horribly unstable when there are many points, and the local method is likely to be more accurate even without such instability, because it has more information than just the function values. Code for doing this interpolation is already in the NR routines.

dpr, if you want dense output at a nonuniformly spaced set of points then you need a modified version of Output. Something along the following lines should do, though it could be streamlined a bit if you didn't mind rewriting more of Output. (Warning: completely untested code; may consist entirely of bugs.)

struct ArbitraryDenseOutput : public Output {
  VecDoub desired_xsave;
  // xs is a vector of values after the first to save at. Must be nonempty and monotonic in same direction as integration.
  ArbitraryDenseOutput(VecDoub_I &xs)
     : Output(xs.size()), desired_xsave(xs) {}
  void out(const Int nstp, const Doub x, VecDoub_I &y, Stepper &s, donst Doub h) {
    if (!dense) throw("dense output not set in Output!");
    if (nstp == -1) { save(x,y); }
    else {
      while (count < desired_xsave.size() && (desired_xsave[count]-xout)*(x2-x1) > 0) {
        save_dense(s,desired_xsave[count],h);
        // save_dense already updates count
      }
    }
  }
};

[tae]

gjm, thanks for the warning about possible bugs in that code. I think I found a few:

if nstp==-1, rather than call save(x,y), do nothing

the 2nd condition for the while loop should read
(x-xout)*(x2-x1) > 0.0

the code in the while loop should read
save_dense(s, xout, h);
xout = desired_xsave[count];

odeint.h with these modifications seems to give the right dense output for an arbitrary array of points to save at. That is, if I choose some of the arbitrary points to coincide with the regularly-spaced points in the original dense output scheme, I get the same solution at those points. It doesn't appear to matter if the arbitrary points are monotonically increasing, decreasing, or neither.