# How can I get a discrete result from NDSolve?

When I use NDsolve to solve a differential equations, it returns an interpolation function of the discrete result. Because I haven't found another way to do Fourier analysis with an interpolation function, it seems I have to do the Fourier transform with Fourier. In my opinion, the process

$\qquad$ discrete result$\rightarrow$ iterpolation function$\rightarrow$ discrete result

will cause severe loss or distortion of the information. So I want to find a way to directly get the discrete result of a differential equations.

My differential equations is like this:

H = {{1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 3, 0}, {0, 0, 0, 4}};
Rho[t_] := Array[Subscript[rho, #1] &, {4}]
t0 = Table[j/10, {j, 1, 10}];
tf[t_] = t + 0.5;
solutionvec =
Table[s =
NDSolve[LogicalExpand[
D[Rho[t], t] == H.Rho[t] && Rho[t0[[j + 1]]] == {1, 0, 0, 0}],
Flatten[Table[Subscript[rho, i], {i, 1, 4}]], {t, t0[[j + 1]],
tf[t0[[j + 1]]]}, MaxStepSize -> 1/100, MaxSteps -> Infinity,
PrecisionGoal -> Infinity, WorkingPrecision -> 30];
{Evaluate[Rho[tf[t0[[j + 1]]]] /. s]}, {j, 0, 9}]

This equations doesn't work, but it can show the thing I'm deal with.

• You are much more likely to get help if you give a simple example of the sort of differential equation you are using. – mikado Jul 9 '17 at 6:13
• Your example doesn't work properly, please recheck it. Also, notice NDSolve doesn't usually use constant step size for solving ODE, while Fourier assumes uniform sampling, so I don't think directly extracting discrete from NDSolve is suitable for solving your problem. – xzczd Jul 9 '17 at 11:39
• Bettertomo, welcome to mma.se. We suggest, as you receive help, try to give it too, by answering questions in your area of expertise. Please take the tour and check the faqs! Please remember, when you see good questions and answers, to vote them up by clicking the gray triangles, and to accept the answer, if any, that solves your problem by clicking the checkmark sign. – kglr Jul 9 '17 at 11:44
• @xzczd, thank you for reminding me about the nonuniform sampling behavior of the NDSolve or it will make me suffer! – Bettertomo Jul 10 '17 at 1:14
• @kglr, thank you for reminding me, I am not very familiar about this excellent forum, I'll try my best! – Bettertomo Jul 10 '17 at 1:16

The problem the OP outlines is to get accurate discrete sample values of the solution (on a uniform grid). The OP is correct that there is error in the discretization by NDSolve (truncation error) and in interpolation (though perhaps not "severe"). The question is whether we can avoid the interpolation error. (1) Just taking the "ValuesOnGrid" either from the interpolating function constructed by NDSolve or by using EvaluationMonitor will get the values of the solution on the grid constructed by NDSolve, which will not be uniform unless a fixed-step method is employed. (2) Using Table to get the values on a regular grid from the solution will introduce interpolation error.

One way to get accurate, regular steps over an interval is to use NDSolveIterate[] and its helpers. (See also Subdivide for subdividing an interval into a specified number of steps.) The code below creates steps at specified regular intervals. (NDSolve will create others, but they are not used.)

ivp = {y'[x] == y[x] Cos[x + y[x]], y[0] == 1};
{state} = NDSolveProcessEquations[ivp, y, {x, 0, 30}];           (* initialize NDSolve *)
samples = 512;                                                    (* number of steps)
Do[NDSolveIterate[state, x], {x, Subdivide[0., 30., samples]}];  (* make steps *)
sol = NDSolveProcessSolutions[state];                            (* get solution *)
yvalues = Table[y[x] /. sol, {x, Subdivide[0., 30., samples]}];   (* get values at steps *)
Fourier[yvalues]
(* etc. *)

The other way, to use a fixed-step method and use if["ValuesOnGrid"] for the input to Fourier[], is also possible, although it can be difficult to get an error estimate to verify that the NDSolve[] solution is accurate. One could solve the system with a high WorkingPrecision and use that as a reference for the exact solution (but that's more than double the work). In any case, let's say you want 137 sample values (or 136 steps). But it turns out that you need at least ~250 steps to get a result accurate to machine precision. So we double the number of steps and downsample the NDSolve[] result:

ivp = {y'[x] == y[x] Cos[x + y[x]], y[0] == 1};
steps = 136;  (* desired number *)
{solIMP} =
NDSolve[ivp, y, {x, 0, 30},
Method -> {"FixedStep",
Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 15,
"ImplicitSolver" -> {"Newton",
AccuracyGoal -> MachinePrecision,
PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}},
StartingStepSize -> (30 - 0)/(2*steps)];          (* double the number of steps *)
fvalues = (y["ValuesOnGrid"] /. solIMP)[[;; ;; 2]];  (* downsample *)
Fourier[fvalues]

You can check that Length[fvalues] is 137.

• Thanks, it's very helpful. – Bettertomo Jul 10 '17 at 1:36
• @Bettertomo I was adding some explanation while you upvoted & commented (thanks!). I think my interpretation of what you needed was slightly different than the others. – Michael E2 Jul 10 '17 at 1:40
• Thanks for your summary ! – Bettertomo Jul 10 '17 at 2:30

Since you don't give any particular differential equation, I am using the 1st one from ref/NDSolve > Basic Examples.

xyPairs = {};
NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30},
EvaluationMonitor :> (xyPairs = Join[xyPairs, {{x, y[x]}}])];
ListPlot[xyPairs]

plot

• Thanks, I think your method is great. But I fail to suit it into my problem for a simple reason: I can't suit $xyPairs = Join[xyPairs, {{x, y[x]}}]$ into my problem... – Bettertomo Jul 9 '17 at 9:11
• I don't know how to express y[x] in your example when I'm dealing my problem, i.e. I don't know the specific way to suit $Evaluate[Rho[tf[t0[[j + 1]]]] /. s]$ in your method, I have give out my case in the question. – Bettertomo Jul 9 '17 at 9:19
• A nice additional trick is to use {x,30,30} for the x-range in these cases: NDSolve it will then only return a result for the endpoint in the form of a (degenerate) interpolating function, which might be much faster and need less memory then building the full interpolating function. As it has to solve for the whole x-range you will still get the same data in xyPairs. The effect is probably neglectable for a single ode but can be quite relevant for a larger system of odes or pdes... – Albert Retey Jul 9 '17 at 9:21
• It's a good way and I will try it . But I still think your first method is really great, so can you help me about it? – Bettertomo Jul 9 '17 at 9:25
• I'm not just curious about it. I just find a severe problem of my program, it's for calculation of absorption spectrum of a two level system. The correction of it will need the first method! – Bettertomo Jul 9 '17 at 9:40

Using the same example as in @m_goldberg's and @zhk's answers:

iF = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}][[1, 1, 2]];

ListPlot @ iF["ValuesOnGrid"]

• Thank you very much, I've made it! – Bettertomo Jul 9 '17 at 11:15
• @Bettertomo, my pleasure. If you haven't already, please consider to check the links in the welcome comment below your post. – kglr Jul 9 '17 at 11:46

Another way is to use Table

ode = y'[x] == y[x] Cos[x + y[x]];

sol = NDSolve[{ode, y[0] == 1}, {y}, {x, 0, 30}];

xy = Table[{x, y[x]} /. First@sol, {x, 0, 30, 2}];

% // TableForm

ListPlot[xy, PlotRange -> All];
• Thank you very much, I know your method is great, but I just can't suit it into my problem cause I know too little about mathematica. I don't know how to express y[x] in your example when I'm dealing my problem, I have give out my case in the question. – Bettertomo Jul 9 '17 at 9:17
• @Bettertomo Your code doesn't work. – zhk Jul 9 '17 at 10:20
• Yes, I fail to find a good simplification, so it just can show what I'm dealing with. The different with the original code is the Matrix H. – Bettertomo Jul 9 '17 at 10:28

A lot of the answers here are reaching into the InterpolatingFunction, to get at the raw output of NDSolve. There are certainly times where that is useful -- but this does not seem to be one of them. I think the original problem of the poster has to do with syntax in NDSolve. In particular, how to deal with vectors.

Here is a stripped down version of the problem in the original post:

H = {{1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 3, 0}, {0, 0, 0, 4}};
eqns = D[rho[t], t] == H.rho[t]
boundary = rho[0] == {1, 0, 0, 0}
sol = NDSolve[{eqns, boundary}, rho, {t, 0, 1}]

To get discrete points, one would write

discrete = Table[{t, #} & /@ (rho[t] /. sol[[1]]), {t, 0, 1, 0.1}]
ListPlot[Transpose@discrete, PlotRange -> All]

Alternatively, he/she could use spectral methods (solving the ODE in Fourier space).

• Thanks, my original problem is trying to avoid the process "raw discrete data in NDsolve $\rightarrow$ interpolation function from NDsolve $\rightarrow$ discrete to suit Fourier". My problem and target change with my knowledge about mathematica in the past ten hours thanks for help from people like you. – Bettertomo Jul 10 '17 at 1:26
• How to deal with vectors is actually my another problems! – Bettertomo Jul 10 '17 at 3:20