I'm using ParametricNDSolve to find a solution to a set of first order ODEs which have one independent variable, z. The set of initial values contains a function of t.

G = 0.00589667;
B = 2066.84;

ode = {y1'[z] == 2/3 * G * y2[z] * y3[z],
       y2'[z] == -B * y3[z] - 2/3 * G * y1[z] * y3[z],
       y3'[z] == B * y2[z]};

iv = {y1[0] == 0., y2[0] == 0., y3[0] == -Exp[-(t - 5)^2]};

ParametricNDSolve treats t as a parameter and returns three parametric interpolating functions with parameter t, independent variable z, and dependent variables y1, y2, and y3:

sol = ParametricNDSolve[{ode, iv}, {y1, y2, y3}, {z, 0., 50.}, t]

{y1 -> ParametricFunction[ <> ],
 y2 -> ParametricFunction[ <> ], 
 y3 -> ParametricFunction[ <> ]}

Then I can evaluate these functions easily to obtain the y value at any t and z in range like this: y1[t][z] /. sol.

I want to plot the y values at a single value of z, over a range of t:

Plot[Evaluate@{y1[t][l] /. sol, y2[t][l] /. sol, y3[t][l] /. sol}, {t, 0, 10}, PlotRange -> All]

This does produce the plot I want, but it takes a very long time, over an hour.

I timed ParametricNDSolve[] and that is not the problem. I then tried a different method of plotting: I evaluated the functions for t and z first, then used ListPlot[]:

y1out = y1[#][50.] & /@ Range[0, 10, 0.01] /. sol;
ListPlot[y1out, PlotRange -> All]

But the improvement was negligible, so I believe it is the evaluation of the parametric functions that is inefficient.

I should note that using the Plot[] function is ideal because it changes its sampling rate with the slope of the curve. I can't predict the shape of the curves, and want to produce the shape of the curves reliably.

I am not a programmer and unfortunately don't have a Mathematica expert around to consult. I do have access to a high performance computer and so I can look into learning how to code in parallel if need be, but I want to know if there is a more efficient way to evaluate and plot a ParametricFunction over a range of values of the parameter?


It is faster to use NDSolve to obtain an interpolation over the parameter than solving the ODE over and over with different boundary conditions:

G = 0.00589667;
B = 2066.84;

ode = {Derivative[1, 0][y1][z, t] == 0.0039311133 y2[z, t] y3[z, t],
       Derivative[1, 0][y2][z, t] == -2066.84 y3[z, t] - 0.0039311133 y3[z, t] Derivative[1, 0][y1][z, t],
       Derivative[1, 0][y3][z, t] == 2066.84 y2[z, t]};

iv = {y1[0, t] == 0., y2[0, t] == 0., y3[0, t] == -Exp[-(t - 5)^2]};

sol = NDSolve[{ode, iv}, {y1, y2, y3}, {z, 0., 50.}, {t, 0., 10}];

Plot[Evaluate[y1[50, t] /. sol], {t, 0, 10}, PlotRange -> All]

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.