How do I make a second order differential equation match a difference equation?

Question

I have fitted a difference equation to some data and would now like to find a differential equation that matches the difference equation. The difference equation and differential equation are

$$ y(n)=f(n) y(n-k)+g(n) y(n-2 k) $$ $$ y''(t)+p(t) y'(t)+q(t) y(t)=0 $$ These are linear difference and differential equations with variable coefficients. f[n] and g[n] are cubic interpolation functions and I would like p[t] and q[t] to also be interpolation functions.

Mathematica can solve symbolic difference equations and differential equations with cubic variable coefficients. Thus

ClearAll[y, f0, f1, f2, f3, g0, g1, g2, g3, n, p0, p1, p2, p3, q0, q1,
   q2, q3, t];
RSolve[y[n] == (f0 + f1 n + f2 n^2 + f3 n^3) y[
     n - 5] + (g0 + g1 n + g2 n^2 + g3 n^3) y[n - 10], y[n], n]
DSolve[y''[t] + (p0 + p1 t + p2 t^2 + p3 t^3) y'[
     t] + (q0 + q1 t + q2 t^2 + q3 t^3) y[t] == 0, y[t], t]

The output is in terms of DifferenceRoot and DifferentialRoot respectively.

For difference and differential equations with constant coefficients the solutions are just

RSolve[y[n] == f0 y[n - k] + g0 y[n - 2 k], y[n], n]
DSolve[y''[t] + p0 y'[t] + q0 y[t] == 0, y[t], t]

So knowing the solution to the difference equation enables the coefficients of the differential equation to be found by simple matching (also time t = (n-1) dt where dt is the sampling increment) how do I do this for my variable coefficients?

Here is some fitted data.

The vertical lines are the locations of knots in the fitting data. One interval of my difference equation is given by

k = 5; (* Backward step increment *)
t1 = 0.07894736842105263`; (* Starting time *)
dt = 0.0004`; (* Time increment *)
dd = {-0.9573507176866912`, -2.024394628072429`, -3.0177598379156`, \
-3.902923724734958`, -4.648984716450747`, -5.22977777434798`, \
-5.624836149234829`, -5.820162645659958`, -5.808779618675913`, \
-5.591034714271433`, -5.174647219824499`}; (* Initial conditions data \
*)
(* Variable coefficients data *)
f = Function[{t}, 
  1.8268003676664075` - 8.64235779514729` t - 
   5.046775951053236` t^2 + 104.88522863362422` t^3];
g = Function[{t}, -1.0031787122707823` + 0.3806275165536718` t - 
    1.0674501046346905` t^2 + 0.7316362694127528` t^3];
ClearAll[y];
ic = Table[y[n] == dd[[n]], {n, 1, 10}];
op = RecurrenceTable[
   Join[{y[n] == 
      f[t1 + (n - 1) dt] y[n - 5] + g[t1 + (n - 1) dt] y[n - 10]}, 
    ic], y[n], {n, 1, 67}];


ListPlot[op, Frame -> True]

How can I get a differential equation with variable coefficients that matches this difference equation? If I can do each interval in the interpolation function in turn as a piecewise function then I have a solution. Thanks for any help.

Answer 1

Fixing ideas for the case with $k = 1$ the correspondence between

$$ \cases{ y(n)=f(n) y(n-k)+g(n) y(n-2 k)\\ y''(t)+p(t) y'(t)+q(t) y(t)=0 } $$

can be obtained as a limit. The cases with $k>1$ would require higher order analogies. So following with $k = 1$ and considering the differential operator approximation $\dot y \approx \frac{y(t+dT)-y(t)}{dT}$ we have

$$ \ddot y + p\dot y+q y \approx \frac{y(t_k)-2y(t_k-dT)+y(t_k-2dT)}{dT^2}+p(t_k)\frac{y(t_k)-y(t_k-dT)}{dT}+q(t_k)y(t_k) $$

then making $j=\frac{t_k}{dT}$ and calling $Y(j) = y(j dT)$ we get the relationship

$$ \ddot y + p\dot y+q y \approx \frac{Y(j)-2Y(j-1)+Y(j-2)}{dT^2}+P(j)\frac{Y(j)-Y(j-1)}{dT}+Q(j)Y(j) $$

and also

$$ \lim_{dT\to 0}\left(\frac{Y(j)-2Y(j-1)+Y(j-2)}{dT^2}+P(j)\frac{Y(j)-Y(j-1)}{dT}+Q(j)Y(j)\right) = \ddot y + p\dot y+q y $$

Now considering the code

Clear["Global`*"]
p[t_] := p0 + p1 t + 0 p2 t^2 + 0 p3 t^3
q[t_] := q0 + q1 t + 0 q2 t^2 + 0 q3 t^3
parms = {p0 -> 1, p1 -> 1, p2 -> 0, p3 -> 0, q0 -> 0, q1 -> 1, q2 -> 1, q3 -> 0, dT -> 1/4}
tmax = 4;
soly = DSolve[y''[t] + p[t] y'[t] + q[t] y[t] == 0, y[t], t][[1]] /. parms // Simplify;
yt = y[t] /. soly;
cond1 = (yt /. {t -> 0}) - 1;
cond2 = (yt /. {t -> tmax});
solC = Solve[{cond1 == 0, cond2 == 0}, {C[1], C[2]}][[1]];
yt0 = yt /. solC;
gr1 = Plot[yt0, {t, 0, tmax}, PlotRange -> All];

equdY = (1 + dT (p0 + dT k p1 + dT (q0 + dT k q1))) Y[k] + Y[k - 2] - (2 + dT (p0 + dT k p1)) Y[k - 1] /. parms;
solY = RSolve[equdY == 0, Y, k][[1]] ;
Ynu = Y[nu] /. solY;
nu0 = Ceiling[tmax/dT /. parms];
cond1 = (Ynu /. {nu -> 0}) - 1;
cond2 = (Ynu /. {nu -> nu0});
solC = Solve[{cond1 == 0, cond2 == 0}, {C[1], C[2]}][[1]];
Ynu0 = Ynu /. solC;
Yk = Table[{nu (dT /. parms), Ynu0}, {nu, 0, nu0}];
gr2 = ListPlot[Yk, PlotStyle -> Red]
Show[gr1, gr2]

we can corroborate those ideas.

NOTE

We condidered that only the cases involving

p[t_] := p0 + p1 t 
q[t_] := q0 + q1 t

could be used due to the existence of solutions in closed form.