# Using NDSolve with Interpolation

I don't think I full understand the nature of NDSolve or Interpolation. I can't understand why the second output doesn't match the first:

func1[{x0_, y0_}, T_] := Module[{}, X[t_] = t; Y[t_] = Sin[t];
nd = NDSolve[{Thread[{x'[t], y'[t]} == {X[t] - x[t], Y[t] - y[t]}/
Sqrt[(X[t] - x[t])^2 + (Y[t] - y[t])^2]], x == x0, y == y0}, {x, y},
{t, 0, T}];Table[Evaluate[{x[t], y[t]} /. nd[]], {t, 0, T}]];

func2[{x0_, y0_}, T_] := Quiet@Module[{}, tableB = N@Table[{t, Sin[t]}, {t, 1, 10, 1}];
X = Quiet[Interpolation /@ Table[Take[tableB[[All, 1]]], {k, 1, 10}]];
Y = Quiet[Interpolation /@ Table[Take[tableB[[All, 2]]], {k, 1, 10}]];
nd = NDSolve[{{Derivative[x][t] == (-x[t] + X[[t]][t])/
Sqrt[(-x[t] + X[[t]][t])^2 + (-y[t] + Y[[t]][t])^2],
Derivative[y][t] == (-y[t] + Y[[t]][t])/Sqrt[(-x[t] + X[[t]][t])^2 +
(-y[t] + Y[[t]][t])^2]}, x == x0, y == y0}, {x, y}, {t, 0, T}];
Quiet@Table[Evaluate[{x[t], y[t]} /. nd[]], {t, 0, T}]];

func1[{2, 2}, 5]
func2[{2, 2}, 5]


What am I doing wrong?

• It looks like there something odd going on with Take[tableB[[All, 1]], y] in your defs of X and Y? Interpolation[tableB[[All, 1]], y] gives you a somewhat sensible output (but not the same as the first, and a different set of error messages). Also, as far as I can tell you're defining tableB for t = 1,... 100, but then using its Interpolation for t = 0,..., T. – aardvark2012 Sep 19 '17 at 1:31
• @aardvark2012 updated, so should make sense, but still no meaningful output – martin Sep 19 '17 at 1:38
• Could you clarify why you want to replace t and Sin[t] with InterpolatingFunctions? Your X function seems to be a very complicated way of writing t. Interpolation[N@Table[{t, Sin[t]}, {t, 1, 100, 1}]] will give you an interpolated Sin function (over the domain [1, 100]) which you could use instead of Y. But could you explain why you'd want to interpolate X = t? – aardvark2012 Sep 19 '17 at 2:06
• @aardvark2012 the data is variable - Sin is just a toy function. The actual data will be far less predictable, and will be cumulative - hence the complicated X and Y functions. If there is a numeric FindRoot or similar alternative, I'd be happy with that. – martin Sep 19 '17 at 2:14
• X[[t]] makes no sense if t is real. – Michael E2 Sep 19 '17 at 3:16

If you can, I would recommend not using Interpolation here. A possible alternative using pure functions passed as arguments is

func4[{x0_, y0_}, T_, X_, Y_] :=
Module[{nd},
nd = First@
NDSolve[{Thread[{x'[t], y'[t]} == {X[t] - x[t], Y[t] - y[t]}/
Sqrt[(X[t] - x[t])^2 + (Y[t] - y[t])^2]], x == x0,
y == y0}, {x, y}, {t, 0, T}];
Table[Evaluate[{x[t], y[t]} /. nd], {t, 0, T}]];


which you can then use as:

ListLinePlot[func4[{2, 2}, 20, # &, Sin]]
ListLinePlot[func4[{2, 2}, 20, Cos, Sin]]
ListLinePlot[func4[{2, 2}, 20, # Cos[#] &, Tan[2 # - 1] &]]   If that's not possible, and you really need Interpolation for some reason, then there are two possibilities I can see, depending on what your goal is with X and Y. The first is that your input function is just Y[t] = Sin[t], in which case interpolating X[t] = t conceptually problematic. The second is that your input functions are parameterized curves in R^2, such as {X[t], Y[t]} = {Cos[t], Sin[t]}, for example.

Case I: Y(t) = Sin(t)

The easy way of doing this (ie, without interpolating anything) would be

func[{x0_, y0_}, T_] :=
Module[{nd},
nd = First@
NDSolve[{Thread[{x'[t], y'[t]} == {t - x[t], Sin[t] - y[t]}/
Sqrt[(t - x[t])^2 + (Sin[t] - y[t])^2]], x == x0,
y == y0}, {x, y}, {t, 0, T}];
Table[Evaluate[{x[t], y[t]} /. nd], {t, 0, T}]
];

ListLinePlot[func[{2, 2}, 20]] To build in Y[t] as an InterpolatingFunction (and leave X[t] = t) just set Y = Interpolation[tableB].

func2[{x0_, y0_}, T_] :=
Module[{tableB = N@Table[{t, Sin[t]}, {t, 0, T, 1}], nd, X, Y},
Y = Interpolation[tableB];
nd = First@
NDSolve[{Thread[{x'[t], y'[t]} == {t - x[t], Y[t] - y[t]}/
Sqrt[(t - x[t])^2 + (Y[t] - y[t])^2]], x == x0,
y == y0}, {x, y}, {t, 0, T}];
Table[Evaluate[{x[t], y[t]} /. nd], {t, 0, T}]];

ListLinePlot[func2[{2, 2}, 20]] which clearly agrees with the plot for the original func.

Case II: {X[t], Y[t]} = {Cos[t], Sin[t]}

For this case, the non-interpolating way would be

func[{x0_, y0_}, T_] :=
Module[{nd},
nd = First@
NDSolve[{Thread[{x'[t], y'[t]} == {Cos[t] - x[t], Sin[t] - y[t]}/
Sqrt[(Cos[t] - x[t])^2 + (Sin[t] - y[t])^2]], x == x0,
y == y0}, {x, y}, {t, 0, T}];
Table[Evaluate[{x[t], y[t]} /. nd], {t, 0, T}]
];

ListLinePlot[func[{2, 2}, 20]] The interpolating version is not much more complicated. Note how tableB is set up, and how X and Y are defined from it.

func3[{x0_, y0_}, T_] :=
Module[{X, Y, nd,
tableB = N@Table[{t, Cos[t], Sin[t]}, {t, 0, T, 1}]},
X = Interpolation[tableB[[;; , {1, 2}]]];
Y = Interpolation[tableB[[;; , {1, 3}]]];
nd = First@
NDSolve[{Thread[{x'[t], y'[t]} == {X[t] - x[t], Y[t] - y[t]}/
Sqrt[(X[t] - x[t])^2 + (Y[t] - y[t])^2]], x == x0,
y == y0}, {x, y}, {t, 0, T}];
Table[Evaluate[{x[t], y[t]} /. nd], {t, 0, T}]];

ListLinePlot[func3[{2, 2}, 17]] Which, again, agrees with the simple version. However, note that I've taken T = 17 and not T = 20. That's because this setup gets very slow very quickly.

• great - thanks for your help. Will try it out for a range of data - I had a feeling it might get slow quickly. – martin Sep 19 '17 at 6:39