ParametricNDSolve
ParametricNDSolve[eqns,u,{x,xmin,xmax},pars]
finds a numerical solution to the ordinary differential equations eqns for the function u with the independent variable x in the range xmin to xmax with parameters pars.
ParametricNDSolve[eqns,u,{x,xmin,xmax},{y,ymin,ymax},pars]
solves the partial differential equations eqns over a rectangular region.
ParametricNDSolve[eqns,u,{x,y}∈Ω,pars]
solves the partial differential equations eqns over the region Ω.
ParametricNDSolve[eqns,u,{t,tmin,tmax},{x,y}∈Ω,pars]
solves the time-dependent partial differential equations eqns over the region Ω.
ParametricNDSolve[eqns,{u1,u2,…},…]
solves for the functions ui.
Details and Options
- ParametricNDSolve gives results in terms of ParametricFunction objects.
- A specification for the parameters pars of {pspec1,pspec2,…} can be used to specify ranges.
- Possible forms for pspeci are:
-
p p has range Reals or Complexes Element[p,Reals] p has range Reals Element[p,Complexes] p has range Complexes Element[p,{v1,…}] p has discrete range {v1,…} {p,pmin,pmax} p has range 
- In ParametricNDSolve[eqns,{u1,u2,…},…], ui can be any expression. Typically, ui will depend on the parameters indirectly through the solution of the differential equations but may depend explicitly on the parameters. A ParametricFunction object that will return a list can be obtained using ParametricNDSolve[eqns,{{u1,u2,…}},…] or by using ParametricNDSolveValue[eqns,{u1,u2,…},…].
- Derivatives of the resulting ParametricFunction objects with respect to the parameters are computed using a combination of symbolic and numerical sensitivity methods when possible.
- ParametricNDSolve takes the same options and settings as NDSolve.
- NDSolve and ParametricNDSolve typically solve differential equations by going through several different stages, depending on the type of equations. With Method->{s1->m1,s2->m2,…}, stage si is handled by method mi. The actual stages used and their order are determined by NDSolve, based on the problem to be solved.
- Possible solution stages are the same as for NDSolve, with the addition of:
-
"ParametricCaching" caching of computed solutions "ParametricSensitivity" computation of derivatives with respect to parameters
List of all options
Examples
open all close allBasic Examples (3)
Get a parametric solution for y:
sol = ParametricNDSolve[{y'[t] == a y[t], y[0] == 1}, y, {t, 0, 10}, {a}]Evaluating with a numerical value of a gives an approximate function solution for y:
y1 = y[1] /. sol
y1[10] /. solPlot the solutions for several different values of the parameter:
Plot[Evaluate[Table[y[a][t] /. sol, {a, -1, 1, .1}]], {t, 0, 1}, PlotRange -> All]Get a function of the parameter a that gives the function f at 
:
sol = ParametricNDSolve[{y''[t] + a y[t] == 0, y[0] == 1, y'[0] == 0}, {y}, {t, 0, 10}, {a}]This plots the value of f[10] as a function of the parameter a:
Plot[Evaluate[y[a][10] /. sol], {a, 0, 2}]Find a value of a for which y[10]=0:
FindRoot[y[a][10] /. sol, {a, 1}]Show the sensitivity of the solution of a differential equation to parameters:
sol = ParametricNDSolve[{y''[t] + a y[t] == 0, y[0] == b, y'[0] == 0}, y, {t, 0, 20}, {a, b}]The sensitivity with respect to a increases with t:
Plot[Evaluate[(y[1, 1][t] + {0, .1, -.1}D[y[a, 1], a][t] /. a -> 1) /. sol], {t, 0, 20}, Filling -> {2 -> {3}}]The sensitivity with respect to b does not increase with t:
Plot[Evaluate[(y[1, 1][t] + {0, .1, -.1}D[y[1, b], b][t] /. b -> 1) /. sol], {t, 0, 20}, Filling -> {2 -> {3}}]Scope (5)
Parameter Dependence (3)
ParametricNDSolve returns a substitution to a ParametricFunction object:
sol = ParametricNDSolve[{y''[t] + y[t] == 3a Sin[y[t]], y[0] == y'[0] == 1}, y, {t, 0, 5}, {a}]
if0 = y[0] /. sol
Plot[if0[t], {t, 0, 5}]Plot solutions for values of 
ranging from 
to 
:
Plot[Evaluate@Table[y[a][t] /. sol, {a, -.5, .5, .05}], {t, 0, 5}, PlotRange -> All]Initial conditions can be specified as parameters:
sol = ParametricNDSolve[{y''[t] + 2y[t] == Sin[.4t] + Cos[y[t]], y[0] == a, y'[0] == b}, y, {t, 0, 10}, {a, b}]Plot solutions with 
with 
for values of 
ranging from 
to 
:
Plot[Evaluate@Table[y[a, 1][t] /. sol, {a, -2, 2, .5}], {t, 0, 10}, PlotRange -> All]Plot solutions with 
with 
for values of 
ranging from 
to 
:
Plot[Evaluate@Table[y[1, b][t] /. sol, {b, -2, 2, .5}], {t, 0, 10}, PlotRange -> All]Differential equation coefficients and boundary conditions can be specified as parameters:
sol = ParametricNDSolve[{u''[x] == a u[x], DirichletCondition[u[x] == 1, x == 0], DirichletCondition[u[x] == b, x == 1]}, u, {x}∈Line[{{0}, {1}}], {a, b}]Plot solutions with 
for values of 
ranging from 
to 
and with 
and 
:
Plot[Evaluate@Table[u[a, 1 / 2][t] /. sol, {a, -2, 2, .5}], {t, 0, 10}, PlotRange -> All]Parameter Sensitivity (2)
Solve the classical harmonic oscillator with parametric amplitude 
:
sol = ParametricNDSolve[{y''[t] + y[t] == 0, y[0] == a, y'[0] == 0}, y, {t, 0, 10}, {a}];Plot the solution for 
and several nearby values of 
:
y2 = y[2] /. sol;
Show[Plot[Evaluate@Table[y[a][t] /. sol, {a, 1, 3, .2}], {t, 0, 10}, PlotStyle -> Thin], Plot[y2[t], {t, 0, 10}]]The sensitivity of 
with respect to 
is by definition 
. Plot the sensitivity at 
:
dy2 = y'[2] /. sol;
Plot[dy2[t], {t, 0, 10}]Plot the sensitivity 
as a band around the solution 
for 
:
Plot[Evaluate[y2[t] + {-1, 0, 1}dy2[t]], {t, 0, 10}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}]Sensitivity analysis of a differential equation with multiple parameters:
sol = ParametricNDSolve[{y''[t] + 2y[t] == Sin[t] + Cos[y[t]], y[0] == a, y'[0] == b}, y, {t, 0, 10}, {a, b}];ifun = y[1, 1] /. solPlot the sensitivity with respect to the initial condition 
at 
, 
:
ifuna = D[y[a, 1], a] /. a -> 1 /. solPlot[Evaluate[ifun[t] + {-.5, 0, .5}ifuna[t]], {t, 0, 10}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}]Plot the sensitivity with respect to the initial condition 
at 
, 
:
ifunb = D[y[1, b], b] /. b -> 1 /. sol;Plot[Evaluate[ifun[t] + {-.5, 0, .5}ifunb[t]], {t, 0, 10}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}]Generalizations & Extensions (1)
Solve 
, 
for various values of WorkingPrecision and plot the error:
sol = ParametricNDSolve[{x'[t] == x[t], x[0] == 1}, x[1], {t, 0, 1}, {wp}, Method -> "Extrapolation", WorkingPrecision -> wp ]ListLogPlot[Table[{wp, Abs[x[1][wp] - E] /. sol}, {wp, 16, 60}]]Options (2)
Method (2)
ParametricCaching (1)
Prevent caching of solutions to save memory:
sol = ParametricNDSolve[{x''[t] + a x[t] == 0, x[0] == b, x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞, Method -> {"ParametricCaching" -> None}];With no caching, the only extra memory required is for the processed equations:
Do[Print[Module[{m = MemoryInUse[]},
AbsoluteTiming[{ByteCount[x[2, 2] /. sol], MemoryInUse[] - m}]]], {3}]
The default is to cache the most recently computed solution:
sol = ParametricNDSolve[{x''[t] + a x[t] == 0, x[0] == b, x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞];With caching, the memory requirement is much greater:
Do[Print[Module[{m = MemoryInUse[]},
AbsoluteTiming[{ByteCount[x[2, 2] /. sol], MemoryInUse[] - m}]]], {3}]
ParametricSensitivity (1)
Specify that no sensitivity should be computed:
sol = ParametricNDSolve[{x'[t] == a x[t], x[0] == b}, x, {t, 0, 1}, {a, b}, Method -> {"ParametricSensitivity" -> None}];The function evaluates quickly:
Plot3D[Evaluate[x[a, b][1] /. sol], {a, 1, 2}, {b, 0, 1}]The derivative is not computed:
dfun = Derivative[1, 0][x /. sol];
dfun[1, 1]
Applications (11)
Parameter Sweeps (5)
Find an initial value 
for which the solution 
of a differential equation will have 
:
sol = ParametricNDSolve[{x''[t] - x'[t] + x[t] == 0, x[0] == 1, x'[0] == s}, x, {t, 0, 30}, {s}];
root = FindRoot[Evaluate[x[s][10] /. sol], {s, 6}]NDSolve[{x''[t] - x'[t] + x[t] == 0, x[0] == 1, x'[0] == s} /. root, x, {t, 0, 10}];Plot[x[t] /. %, {t, 0, 10}]Compare to nearby values of the parameter s:
Show[Plot[Evaluate@Table[Tooltip[x[s][t], s] /. sol, {s, 1, 1.8, 0.1}], {t, 0, 10}, PlotStyle -> Lighter[Gray, 0.5]],
Plot[Evaluate[x[s][t] /. sol /. root], {t, 0, 10}, PlotStyle -> Orange]]Find several solutions to the boundary value problem 
, 
, 
. First consider several possible values for 
:
sol = ParametricNDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 0, y'[0] == s}, y, {t, 0, 10}, {s}];Plot[Evaluate@Table[Tooltip[y[s][t] /. sol, s], {s, 0.3, 1.9, 0.1}], {t, 0, 10}]Run a parameter sweep to determine nontrivial solution values of 
:
Plot[(y[s] /. sol)[10], {s, 0, 2.1}]Using approximate initial values from the graph above:
val = Table[FindRoot[(y[s] /. sol)[10], {s, s0}], {s0, {0.9, 1.9, 2}}]Plot the solutions that were found:
Plot[Evaluate[y[s][t] /. sol /. val], {t, 0, 10}]Find all eigenvalues 
and eigenfunctions 
for 
with 
. Start by exploring the possible parameter values:
psol = y /. ParametricNDSolve[{y''[x] == -λ^2 y[x], y[0] == 0, y'[0] == λ}, y, {x, 0, 1}, {λ}];Plot[Evaluate@Table[Tooltip[psol[λ][x], λ], {λ, π / 2, 11π / 2, π / 2}], {x, 0, 1}]
pfun = y[1] /. ParametricNDSolve[{y''[x] == -λ^2 y[x], y[0] == 0, y'[0] == λ}, y[1], {x, 0, 1}, {λ}];Plot[pfun[λ], {λ, π / 2, 11π / 2}]sol = Table[FindRoot[pfun[λ] == 0, {λ, λ0}], {λ0, {3, 6, 9, 12, 16}}]//QuietPlot[Evaluate[psol[λ][x] /. sol], {x, 0, 1}]Find the value of 
for which the solution of 
, 
has minimal arc length from 
to 
. Begin by plotting the solutions for values of 
ranging from 0 to 1:
psol = ParametricNDSolve[{y''[t] + y[t] == 3a Sin[y[t]], y[0] == y'[0] == 1, α'[t] == Sqrt[1 + y[t] ^ 2], α[0] == 0}, {y, α}, {t, 0, 10}, {a}];Plot[Evaluate@Table[y[a][t] /. psol, {a, 0, 1, .1}], {t, 0, 10}]Plot 
versus the arc length of the solution:
Plot[Evaluate[α[a][10] /. psol], {a, 0, 1}, PlotPoints -> 100]The minimum arc length solution for 
seems to occur at 
:
Show[Plot[Evaluate@Table[y[a][t] /. psol, {a, 0, 1, .1}], {t, 0, 10}, PlotRange -> All, PlotPoints -> 500, PlotStyle -> Lighter[Gray, 0.3]], Plot[Evaluate[y[0][t] /. psol], {t, 0, 10}, PlotStyle -> Thick]]Find the local minimum, which appears near 
:
lmin = FindMinimum[α[a][10] /. psol, {a, 0.72}, AccuracyGoal -> 7]Plot the corresponding solution of (locally) minimal arc length together with some nearby solutions:
Show[{Plot[Evaluate@Table[y[a][t] /. psol, {a, .65, .8, .005}], {t, 0, 10}, PlotStyle -> Lighter[Gray, 0.5]], Plot[Evaluate[y[a][t] /. psol /. lmin[[2]]], {t, 0, 10}, PlotStyle -> Thick]}]Solve a radial form of the Klein-Gordon equation for a potential of the form 
:
sol = ParametricNDSolve[{Derivative[2][R][r] - ((1 - (UnitConvert[Quantity[1, "FineStructureConstant"]] E^-r Z + ϵ)^2) R[r]) == 0, R[0] == 0, R'[0] == 1}, R, {r, 0, 50}, {ϵ, Z}]Plot the values of the energy 
and the strength parameter 
that make the solution 
decay at the outer boundary:
ContourPlot[(R[ϵ, Z][50] /. sol) == 0, {Z, 50, 800}, {ϵ, -1, 1}, FrameLabel -> Automatic, PlotPoints -> 100]Parameter Sensitivities (4)
Perturb a parameter in a differential equation and view several of the resulting perturbed solutions:
pfun = y /. ParametricNDSolve[{y''[t] + y[t] == a Sin[2t ] + Cos[2y[t]], y[0] == 1, y'[0] == 0}, y, {t, 0, 10}, {a}];Show[Plot[Evaluate@Table[pfun[a][t], {a, -.6, .6, .1}], {t, 0, 5}, PlotStyle -> Lighter[Gray, 0.3], PlotRange -> All], Plot[pfun[0][t], {t, 0, 5}]]Plotting the sensitivity solutions gives qualitatively the same result:
Plot[Evaluate[pfun[0][t] + {-.5, 0, .5}pfun'[0][t]], {t, 0, 5}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}]Simulate an inverted pendulum stabilized by an oscillating base:
pfun = θ /. ParametricNDSolve[{θ''[t] + (9.8 / 10)Sin[θ[t]] == -(a / 10)w ^ 2 Sin[w t]Sin[θ[t]], θ[0] == 3.14158, θ'[0] == 0.01}, θ, {t, 0, 5}, {w, a}];Sensitivity of 
with respect to the amplitude a increases with time:
if = pfun[50, 2];
ifa = D[pfun[w, a], a] /. {w -> 50, a -> 2};Plot[Evaluate[if[t] + {-.1, 0, .1}ifa[t]], {t, 0, 5}, PlotStyle -> {Gray, Blue, Gray}, Filling -> {1 -> {3}}]Parametric dependence of the heat equation 
, 
:
pfun = u /. ParametricNDSolve[{D[u[t, x], t] == c^2 D[u[t, x], x, x], u[0, x] == Exp[-(a x) ^ 2], u[t, -10] == u[t, 10] == 0}, u, {t, 0, 5}, {x, -10, 10}, {a, c}];
ifun = pfun[1, 1];
ifda = D[pfun[a, 1], a] /. {a -> 1};
ifdc = D[pfun[1, c], c] /. {c -> 1};Plot the corresponding sensitivity bands:
Panel@Grid[Join[{{"", "a", "c"}}, Table[{Row[{"t=", τ}], Plot[Evaluate[(ifun[τ, x] + .5{0, 1, -1}ifda[τ, x])], {x, -10, 10}, Filling -> {2 -> {3}}, PlotRange -> All], Plot[Evaluate[(ifun[τ, x] + .5{0, 1, -1}ifdc[τ, x])], {x, -10, 10}, Filling -> {2 -> {3}}, PlotRange -> All]}, {τ, {1, 3, 5}}]]]Indicate sensitivity to a and c by changing the color of the solution surface:
ColorSensitivityPlot3D[f_, fsense_, opts___] :=
Block[{fs = fsense}, Plot3D[f, {t, 0, 5}, {x, -10, 10}, Mesh -> None, ColorFunction -> Function[{t, x}, ColorData["LakeColors"][fs]], PlotPoints -> 60, ColorFunctionScaling -> False, opts]]{ColorSensitivityPlot3D[pfun[1, 1][t, x], 1.5Abs[ifda[t, x]], PlotLabel -> a, PlotRange -> All], ColorSensitivityPlot3D[pfun[1, 1][t, x], Abs[ifdc[t, x]], PlotLabel -> c, PlotRange -> All]}Find the sensitivity of the Lorenz equations to a parameter:
ndf = x[T] /. ParametricNDSolve[{Derivative[1][x][t] == -3 (x[t] - y[t]), Derivative[1][y][t] == -x[t] z[t] + a x[t] - y[t], Derivative[1][z][t] == x[t] y[t] - z[t], x[0] == z[0] == 0, y[0] == 1}, x[T], {t, 0, T}, {a, T}, MaxSteps -> ∞];Plot[Evaluate[RealExponent[D[ndf[a, 8], a]]], {a, 26, 27}, PlotRange -> All]Parameter Fitting (2)
Sample the solution of a differential equation and add noise:
{sol, samples} = NDSolve[{y''[t] == -9.8, y[0] == 1, y'[0] == 40, WhenEvent[Mod[t, .5] == 0, Sow[{t, y[t]}]]}, y, {t, 0, 10}]//Reap;
noisysamples = Map[({#[[1]], #[[2]] + RandomReal[{-10, 10}]})&, First@samples];ListPlot[noisysamples]Fit a trigonometric model to the noisy data:
pfun = y /. ParametricNDSolve[{y''[t] == -w y[t], y[0] == a, y'[0] == b}, y, {t, 0, 10}, {w, a, b}];fit = FindFit[noisysamples, pfun[w, a, b][t], {{w, .1}, {a, 1}, {b, 1}}, t]Plot[pfun[w, a, b][t] /. fit, {t, 0, 10}, Epilog -> {Point@noisysamples}]A quadratic model is a better fit:
pfun = y /. ParametricNDSolve[{y''[t] == -g, y[0] == p, y'[0] == v}, y, {t, 0, 10}, {g, p, v}];fit = FindFit[noisysamples, pfun[g, p, v][t], {{g, 5}, {p, .5}, {v, .5}}, t]Plot[pfun[g, p, v][t] /. fit, {t, 0, 10}, Epilog -> {Point@noisysamples}]Find the parameters that make the solution of a differential equation the best fit to data:
time = Quantity[{0, 25, 34, 42, 54, 63, 80, 89, 103, 115, 131, 145, 158, 175, 195, 213, 225, 250, 274, 298, 315, 335, 353, 387, 411, 440, 475, 492, 520, 530, 550, 560, 570, 585, 600, 610, 620, 630, 640, 647, 660, 670, 680, 690, 700, 710, 720, 730, 740, 750, 760, 770, 778, 790, 800, 810, 820, 830, 840, 850, 860, 870, 880, 890, 900, 910, 920, 930, 940, 950, 960, 970, 980, 990, 1000, 1010, 1020, 1030, 1040, 1050, 1060, 1070, 1080, 1090, 1100, 1110, 1120, 1130, 1140, 1150, 1160, 1170, 1180, 1190, 1200, 1210, 1220, 1230, 1240, 1250, 1260, 1270, 1280, 1300, 1340, 1356, 1380, 1418, 1725, 1740, 1825, 1890, 1955, 2022, 2077, 2160, 2244, 2320, 2408, 2473, 2613}, "Seconds"];
temp = Quantity[{210, 204, 200, 198, 196, 194, 192, 190, 188, 186, 184, 182, 180, 178, 176, 174, 172, 170, 168, 166, 164, 162, 160, 158, 156, 154, 152, 150, 150, 148, 148, 147, 146, 145, 144, 144, 144, 143, 143, 142, 142, 141, 141, 140, 140, 139.5, 139, 139, 138, 138, 137, 137, 136, 136, 136, 135, 135, 134, 134, 133, 133, 133, 132.5, 132, 132, 131, 131, 131, 130, 130, 129.5, 129, 129, 128.5, 128, 128, 128, 127.5, 127, 127, 126, 126, 126, 125.5, 125, 125, 124.5, 124, 124, 123.5, 123, 123, 122.5, 122, 122, 122, 121.5, 121, 121, 121, 120.5, 120.5, 120, 120, 119, 118, 117, 116, 111, 110.5, 109, 108, 107, 106, 105, 104, 103, 102, 101, 100, 99}, "DegreesFahrenheit"];The data comes from a cooling body, so use Newton's law of cooling:
BoltzmannConstant = Quantity[1.3806504`*^-23, "Joules" / "Kelvins"];tempK = UnitConvert[temp, "Kelvins"];
fitdata = QuantityMagnitude[Transpose[{time, tempK}]];
Subscript[tempK, 0] = QuantityMagnitude[First[tempK]];
Subscript[tempK, ∞] = QuantityMagnitude[UnitConvert[Quantity[79, "DegreesFahrenheit"], "Kelvins"]];c = QuantityMagnitude[(2260 18/6.02 10^23 BoltzmannConstant)];
tend = QuantityMagnitude[Last[time]];pfun = T /. ParametricNDSolve[{(Derivative[1][T][t] == -K1 (T[t] - Subscript[tempK, ∞]) - K2 (E^-c / T[t] - E^-c / Subscript[tempK, ∞])), T[0] == Subscript[tempK, 0]}, T, {t, 0, tend}, {K1, K2}];fit = FindFit[fitdata//N, pfun[K1, K2][t], {{K1, 0.0001}, {K2, 10000}}, t]Show[ListPlot[fitdata], Plot[Evaluate[pfun[K1, K2][t] /. fit], {t, 0, tend}, PlotStyle -> Red]]Properties & Relations (1)
Use SystemModelParametricSimulate to simulate larger hierarchical system models:
pf = SystemModelParametricSimulate[\!\(\*GraphicsBox[«8»]\), "Inertia2.w", 40, {"Resistor1.R", "SpringDamper1.d"}]Simulate with two sets of resistor and spring damper parameters:
sims = {pf[0.1, 0.25], pf[0.4, 0.3]};Compare the resulting angular velocities:
Plot[Evaluate@Comap[sims, t], {t, 0, 40}]Text
Wolfram Research (2012), ParametricNDSolve, Wolfram Language function, https://reference.wolfram.com/language/ref/ParametricNDSolve.html (updated 2014).
CMS
Wolfram Language. 2012. "ParametricNDSolve." Wolfram Language & System Documentation Center. Wolfram Research. Last Modified 2014. https://reference.wolfram.com/language/ref/ParametricNDSolve.html.
APA
Wolfram Language. (2012). ParametricNDSolve. Wolfram Language & System Documentation Center. Retrieved from https://reference.wolfram.com/language/ref/ParametricNDSolve.html