--- title: "ParametricNDSolve" language: "en" type: "Symbol" summary: "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} \\[Element] \\[CapitalOmega], pars] solves the partial differential equations eqns over the region \\[CapitalOmega]. ParametricNDSolve[eqns, u, {t, tmin, tmax}, {x, y} \\[Element] \\[CapitalOmega], pars] solves the time-dependent partial differential equations eqns over the region \\[CapitalOmega]. ParametricNDSolve[eqns, {u1, u2, ...}, ...] solves for the functions ui." keywords: - sensitivity - sensitivity of differential equations - sensitivity solutions - differential sensitivity - parameter sweeps - system identification - model calibration - parametric differential equations - parametric dependence of solutions - parameters in differential equations - adjoint sensitivity - forward sensitivity - IDAS - CVODES canonical_url: "https://reference.wolfram.com/language/ref/ParametricNDSolve.html" source: "Wolfram Language Documentation" related_guides: - title: "Equation Solving" link: "https://reference.wolfram.com/language/guide/EquationSolving.en.md" - title: "Differential Equations" link: "https://reference.wolfram.com/language/guide/DifferentialEquations.en.md" - title: "Systems Modeling" link: "https://reference.wolfram.com/language/guide/SystemsModeling.en.md" related_functions: - title: "ParametricNDSolveValue" link: "https://reference.wolfram.com/language/ref/ParametricNDSolveValue.en.md" - title: "ParametricFunction" link: "https://reference.wolfram.com/language/ref/ParametricFunction.en.md" - title: "NDSolve" link: "https://reference.wolfram.com/language/ref/NDSolve.en.md" - title: "DSolve" link: "https://reference.wolfram.com/language/ref/DSolve.en.md" - title: "SystemModelParametricSimulate" link: "https://reference.wolfram.com/language/ref/SystemModelParametricSimulate.en.md" --- # 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 $p_{\min }\leq p\leq p_{\max }$ | * 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 | | | | | ------------------ | ---------------- | ---------------------------------------------------------- | | AccuracyGoal | Automatic | digits of absolute accuracy sought | | Compiled | Automatic | whether expressions should be compiled automatically | | DependentVariables | Automatic | the list of all dependent variables | | EvaluationMonitor | None | expression to evaluate whenever the function is evaluated | | InitialSeeding | {} | seeding equations for some algorithms | | InterpolationOrder | Automatic | the continuity degree of the final output | | MaxStepFraction | 1 / 10 | maximum fraction of range to cover in each step | | MaxSteps | Automatic | maximum number of steps to take | | MaxStepSize | Automatic | maximum size of each step | | Method | Automatic | method to use | | NormFunction | Automatic | the norm to use for error estimation | | PrecisionGoal | Automatic | digits of precision sought | | StartingStepSize | Automatic | initial step size used | | StepMonitor | None | expression to evaluate when a step is taken | | WorkingPrecision | MachinePrecision | precision to use in internal computations | ## Examples (22) ### Basic Examples (3) Get a parametric solution for ``y`` : ```wl In[1]:= sol = ParametricNDSolve[{y'[t] == a y[t], y[0] == 1}, y, {t, 0, 10}, {a}] Out[1]= {y -> ParametricFunction[<>]} ``` Evaluating with a numerical value of ``a`` gives an approximate function solution for ``y`` : ```wl In[2]:= y1 = y[1] /. sol Out[2]= InterpolatingFunction[{{0., 10.}}, {5, 7, 1, {86}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«1017»"], {Developer`PackedArrayForm, {0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 4 ... , 94, 96, 98, 100, 102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140, 142, 144, 146, 148, 150, 152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172}, CompressedData["«1608»"]}, {Automatic}] ``` Evaluate at a time $t=10$ : ```wl In[3]:= y1[10] /. sol Out[3]= 22026.5 ``` Plot the solutions for several different values of the parameter: ```wl In[4]:= Plot[Evaluate[Table[y[a][t] /. sol, {a, -1, 1, .1}]], {t, 0, 1}, PlotRange -> All] Out[4]= [image] ``` --- Get a function of the parameter ``a`` that gives the function ``f`` at $t=10$ : ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + a y[t] == 0, y[0] == 1, y'[0] == 0}, {y}, {t, 0, 10}, {a}] Out[1]= {y -> ParametricFunction[<>]} ``` This plots the value of ``f[10]`` as a function of the parameter ``a`` : ```wl In[2]:= Plot[Evaluate[y[a][10] /. sol], {a, 0, 2}] Out[2]= [image] ``` Find a value of ``a`` for which ``y[10] = 0`` : ```wl In[3]:= FindRoot[y[a][10] /. sol, {a, 1}] Out[3]= {a -> 1.20903} ``` --- Show the sensitivity of the solution of a differential equation to parameters: ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + a y[t] == 0, y[0] == b, y'[0] == 0}, y, {t, 0, 20}, {a, b}] Out[1]= {y -> ParametricFunction[<>]} ``` The sensitivity with respect to ``a`` increases with ``t`` : ```wl In[2]:= Plot[Evaluate[(y[1, 1][t] + {0, .1, -.1}D[y[a, 1], a][t] /. a -> 1) /. sol], {t, 0, 20}, Filling -> {2 -> {3}}] Out[2]= [image] ``` The sensitivity with respect to ``b`` does not increase with ``t``: ```wl In[3]:= Plot[Evaluate[(y[1, 1][t] + {0, .1, -.1}D[y[1, b], b][t] /. b -> 1) /. sol], {t, 0, 20}, Filling -> {2 -> {3}}] Out[3]= [image] ``` ### Scope (5) #### Parameter Dependence (3) ``ParametricNDSolve`` returns a substitution to a ``ParametricFunction`` object: ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + y[t] == 3a Sin[y[t]], y[0] == y'[0] == 1}, y, {t, 0, 5}, {a}] Out[1]= {y -> ParametricFunction[<>]} ``` Get the solution for $a=0$: ```wl In[2]:= if0 = y[0] /. sol Out[2]= InterpolatingFunction[{{0., 5.}}, {5, 7, 2, {57}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«693»"], {Developer`PackedArrayForm, {0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51, 54, 57, 60, 63, 66, 69, 72, 75, 78, 81, 84, 87, 90, 93, 96, 99, 102, 105, 108, 111, 114, 117, 120, 123, 126, 129, 132, 135, 138, 141, 144, 147, 150, 153, 156, 159, 162, 165, 168, 171}, CompressedData["«1750»"]}, {Automatic}] ``` Plot the solution for $a=0$ : ```wl In[3]:= Plot[if0[t], {t, 0, 5}] Out[3]= [image] ``` Plot solutions for values of $a$ ranging from $-0.5$ to $0.5$: ```wl In[4]:= Plot[Evaluate@Table[y[a][t] /. sol, {a, -.5, .5, .05}], {t, 0, 5}, PlotRange -> All] Out[4]= [image] ``` --- Initial conditions can be specified as parameters: ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + 2y[t] == Sin[.4t] + Cos[y[t]], y[0] == a, y'[0] == b}, y, {t, 0, 10}, {a, b}] Out[1]= {y -> ParametricFunction[<>]} ``` Plot solutions with $y(0)=a$ with $y'(0)=1$ for values of $a$ ranging from $-2$ to $2$ : ```wl In[2]:= Plot[Evaluate@Table[y[a, 1][t] /. sol, {a, -2, 2, .5}], {t, 0, 10}, PlotRange -> All] Out[2]= [image] ``` Plot solutions with $y(0)=1$ with $y'(0)=b$ for values of $b$ ranging from $-2$ to $2$ : ```wl In[3]:= Plot[Evaluate@Table[y[1, b][t] /. sol, {b, -2, 2, .5}], {t, 0, 10}, PlotRange -> All] Out[3]= [image] ``` --- Differential equation coefficients and boundary conditions can be specified as parameters: ```wl In[1]:= 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}] Out[1]= {u -> ParametricFunction[<>]} ``` Plot solutions with $u''(x)=a u(x)$ for values of $a$ ranging from $-2$ to $2$ and with $u(0)=1$ and $u(1)=1/2$ : ```wl In[2]:= Plot[Evaluate@Table[u[a, 1 / 2][t] /. sol, {a, -2, 2, .5}], {t, 0, 10}, PlotRange -> All] Out[2]= [image] ``` #### Parameter Sensitivity (2) Solve the classical harmonic oscillator with parametric amplitude $a$ : ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + y[t] == 0, y[0] == a, y'[0] == 0}, y, {t, 0, 10}, {a}]; ``` Plot the solution for $a=2$ and several nearby values of $a$ : ```wl In[2]:= 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}]] Out[2]= [image] ``` The sensitivity of $y$ with respect to $a$ is by definition $y_a(t,a)$. Plot the sensitivity at $a=2$ : ```wl In[3]:= dy2 = y'[2] /. sol; Plot[dy2[t], {t, 0, 10}] Out[3]= [image] ``` Plot the sensitivity $y_a(t,a)$ as a band around the solution $y(t,a)$ for $a=2$ : ```wl In[4]:= Plot[Evaluate[y2[t] + {-1, 0, 1}dy2[t]], {t, 0, 10}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}] Out[4]= [image] ``` --- Sensitivity analysis of a differential equation with multiple parameters: ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + 2y[t] == Sin[t] + Cos[y[t]], y[0] == a, y'[0] == b}, y, {t, 0, 10}, {a, b}]; In[2]:= ifun = y[1, 1] /. sol Out[2]= InterpolatingFunction[{{0., 10.}}, {5, 7, 2, {157}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«1769»"], {Developer`PackedArrayForm, CompressedData["«438»"], CompressedData["«5166»"]}, {Automatic}] ``` Plot the sensitivity with respect to the initial condition $y(0)=a$ at $a=1$, $b=1$ : ```wl In[3]:= ifuna = D[y[a, 1], a] /. a -> 1 /. sol Out[3]= InterpolatingFunction[{{0., 10.}}, {5, 7, 2, {308}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, CompressedData["«3393»"], {Developer`PackedArrayForm, CompressedData["«726»"], CompressedData["«10044»"]}, {Automatic}] In[4]:= Plot[Evaluate[ifun[t] + {-.5, 0, .5}ifuna[t]], {t, 0, 10}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}] Out[4]= [image] ``` Plot the sensitivity with respect to the initial condition $y'(0)=b$ at $a=1$, $b=1$ : ```wl In[5]:= ifunb = D[y[1, b], b] /. b -> 1 /. sol; In[6]:= Plot[Evaluate[ifun[t] + {-.5, 0, .5}ifunb[t]], {t, 0, 10}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}] Out[6]= [image] ``` ### Generalizations & Extensions (1) Solve $x'(t)=x(t)$, $x(0)=1$ for various values of ``WorkingPrecision`` and plot the error: ```wl In[1]:= sol = ParametricNDSolve[{x'[t] == x[t], x[0] == 1}, x[1], {t, 0, 1}, {wp}, Method -> "Extrapolation", WorkingPrecision -> wp ] Out[1]= {x[1] -> ParametricFunction[<>]} In[2]:= ListLogPlot[Table[{wp, Abs[x[1][wp] - E] /. sol}, {wp, 16, 60}]] Out[2]= [image] ``` ### Options (2) #### Method (2) ##### ParametricCaching (1) Prevent caching of solutions to save memory: ```wl In[1]:= 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: ```wl In[2]:= Do[Print[Module[{m = MemoryInUse[]}, AbsoluteTiming[{ByteCount[x[2, 2] /. sol], MemoryInUse[] - m}]]], {3}] During evaluation of In[2]:= {0.039815, {449832, 381888}} During evaluation of In[2]:= {0.035419, {449832, 208}} During evaluation of In[2]:= {0.035642, {449832, 152}} ``` The default is to cache the most recently computed solution: ```wl In[3]:= 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: ```wl In[4]:= Do[Print[Module[{m = MemoryInUse[]}, AbsoluteTiming[{ByteCount[x[2, 2] /. sol], MemoryInUse[] - m}]]], {3}] During evaluation of In[4]:= {0.035887, {449832, 497392}} During evaluation of In[4]:= {0.000019, {449832, 56}} During evaluation of In[4]:= {9.00000000000000022800771687370158`0.9748424227189488*^-6, {449832, 56}} ``` ##### ParametricSensitivity (1) Specify that no sensitivity should be computed: ```wl In[1]:= sol = ParametricNDSolve[{x'[t] == a x[t], x[0] == b}, x, {t, 0, 1}, {a, b}, Method -> {"ParametricSensitivity" -> None}]; ``` The function evaluates quickly: ```wl In[2]:= Plot3D[Evaluate[x[a, b][1] /. sol], {a, 1, 2}, {b, 0, 1}] Out[2]= [image] ``` The derivative is not computed: ```wl In[3]:= dfun = Derivative[1, 0][x /. sol]; dfun[1, 1] ``` ParametricNDSolve::scompute: The sensitivity derivative \!\(\*SubscriptBox[\(\[PartialD]\), \(a\)]x\) has total degree exceeding the specification in ParametricSensitivity->None. ```wl Out[4]= ParametricFunction[<>][1, 1] ``` ### Applications (10) #### Parameter Sweeps (4) Find an initial value $x'(0)$ for which the solution $x(t)$ of a differential equation will have $x(10)=0$ : ```wl In[1]:= sol = ParametricNDSolve[{x''[t] - x'[t] + x[t] == 0, x[0] == 1, x'[0] == s}, x, {t, 0, 30}, {s}]; ``` Find a solution with $x(10)=0$: ```wl In[2]:= root = FindRoot[Evaluate[x[s][10] /. sol], {s, 6}] Out[2]= {s -> 1.40296} ``` Check the resulting solution: ```wl In[3]:= NDSolve[{x''[t] - x'[t] + x[t] == 0, x[0] == 1, x'[0] == s} /. root, x, {t, 0, 10}]; In[4]:= Plot[x[t] /. %, {t, 0, 10}] Out[4]= [image] ``` Compare to nearby values of the parameter ``s``: ```wl In[5]:= 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]] Out[5]= [image] ``` --- Find several solutions to the boundary value problem $y''(t)=\sin (y(t))$, $y(0)=0$, $y(10)=0$. First consider several possible values for $y'(0)$: ```wl In[1]:= sol = ParametricNDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 0, y'[0] == s}, y, {t, 0, 10}, {s}]; In[2]:= Plot[Evaluate@Table[Tooltip[y[s][t] /. sol, s], {s, 0.3, 1.9, 0.1}], {t, 0, 10}] Out[2]= [image] ``` Run a parameter sweep to determine nontrivial solution values of $y'(0)$ : ```wl In[3]:= Plot[(y[s] /. sol)[10], {s, 0, 2.1}] Out[3]= [image] ``` Using approximate initial values from the graph above: ```wl In[4]:= val = Table[FindRoot[(y[s] /. sol)[10], {s, s0}], {s0, {0.9, 1.9, 2}}] Out[4]= {{s -> 0.924845}, {s -> 1.87817}, {s -> 1.99927}} ``` Plot the solutions that were found: ```wl In[5]:= Plot[Evaluate[y[s][t] /. sol /. val], {t, 0, 10}] Out[5]= [image] ``` --- Find all eigenvalues $\pi /2<\lambda _i<11 \pi /2$ and eigenfunctions $y_i(x)$ for $y''(x)=-\lambda ^2 y(x)$ with $y(0)=y(1)=0$. Start by exploring the possible parameter values: ```wl In[1]:= psol = y /. ParametricNDSolve[{y''[x] == -λ^2 y[x], y[0] == 0, y'[0] == λ}, y, {x, 0, 1}, {λ}]; In[2]:= Plot[Evaluate@Table[Tooltip[psol[λ][x], λ], {λ, π / 2, 11π / 2, π / 2}], {x, 0, 1}] Out[2]= [image] ``` Find the exact values for $\lambda _i$: ```wl In[3]:= pfun = y[1] /. ParametricNDSolve[{y''[x] == -λ^2 y[x], y[0] == 0, y'[0] == λ}, y[1], {x, 0, 1}, {λ}]; In[4]:= Plot[pfun[λ], {λ, π / 2, 11π / 2}] Out[4]= [image] In[5]:= sol = Table[FindRoot[pfun[λ] == 0, {λ, λ0}], {λ0, {3, 6, 9, 12, 16}}]//Quiet Out[5]= {{λ -> 3.14159}, {λ -> 6.28319}, {λ -> 9.42478}, {λ -> 12.5664}, {λ -> 15.708}} ``` Plot them: ```wl In[6]:= Plot[Evaluate[psol[λ][x] /. sol], {x, 0, 1}] Out[6]= [image] ``` --- Find the value of $a$ for which the solution of $=3 a \sin (y(t))$, $y(0)=y'(0)=1$ has minimal arc length from $t=0$ to $t=10$. Begin by plotting the solutions for values of $a$ ranging from 0 to 1: ```wl In[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}]; In[2]:= Plot[Evaluate@Table[y[a][t] /. psol, {a, 0, 1, .1}], {t, 0, 10}] Out[2]= [image] ``` Plot $a$ versus the arc length of the solution: ```wl In[3]:= Plot[Evaluate[α[a][10] /. psol], {a, 0, 1}, PlotPoints -> 100] Out[3]= [image] ``` The minimum arc length solution for $0\leq a\leq 1$ seems to occur at $a=0$ : ```wl In[4]:= 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]] Out[4]= [image] ``` Find the local minimum, which appears near $a=0.725$ : ```wl In[5]:= lmin = FindMinimum[α[a][10] /. psol, {a, 0.72}, AccuracyGoal -> 7] Out[5]= {14.9146, {a -> 0.725089}} ``` Plot the corresponding solution of (locally) minimal arc length together with some nearby solutions: ```wl In[6]:= 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]}] Out[6]= [image] ``` #### Parameter Sensitivities (4) Perturb a parameter in a differential equation and view several of the resulting perturbed solutions: ```wl In[1]:= pfun = y /. ParametricNDSolve[{y''[t] + y[t] == a Sin[2t ] + Cos[2y[t]], y[0] == 1, y'[0] == 0}, y, {t, 0, 10}, {a}]; In[2]:= 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}]] Out[2]= [image] ``` Plotting the sensitivity solutions gives qualitatively the same result: ```wl In[3]:= Plot[Evaluate[pfun[0][t] + {-.5, 0, .5}pfun'[0][t]], {t, 0, 5}, PlotStyle -> {Gray, {Blue}, Gray}, Filling -> {1 -> {3}}] Out[3]= [image] ``` --- Simulate an inverted pendulum stabilized by an oscillating base: ```wl In[1]:= 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 $\theta$ with respect to the amplitude ``a`` increases with time: ```wl In[2]:= if = pfun[50, 2]; ifa = D[pfun[w, a], a] /. {w -> 50, a -> 2}; In[3]:= Plot[Evaluate[if[t] + {-.1, 0, .1}ifa[t]], {t, 0, 5}, PlotStyle -> {Gray, Blue, Gray}, Filling -> {1 -> {3}}] Out[3]= [image] ``` --- Parametric dependence of the heat equation $u_t(t,x)=c^2 u_{\text{xx}}(t,x)$, $u(0,x)=\exp \left(-a x^2\right)$ : ```wl In[1]:= 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}]; ``` Find the sensitivities $\frac{\partial u(t,x,a,c)}{\partial a}$ and $\frac{\partial u(t,x,a,c)}{\partial c}$: ```wl In[2]:= 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: ```wl In[3]:= 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}}]]] Out[3]= [image] ``` Indicate sensitivity to ``a`` and ``c`` by changing the color of the solution surface: ```wl In[4]:= 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]] In[5]:= {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]} Out[5]= [image] ``` --- Find the sensitivity of the Lorenz equations to a parameter: ```wl In[1]:= 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 -> ∞]; In[2]:= Plot[Evaluate[RealExponent[D[ndf[a, 8], a]]], {a, 26, 27}, PlotRange -> All] Out[2]= [image] ``` #### Parameter Fitting (2) Sample the solution of a differential equation and add noise: ```wl In[1]:= {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]; In[2]:= ListPlot[noisysamples] Out[2]= [image] ``` Fit a trigonometric model to the noisy data: ```wl In[3]:= pfun = y /. ParametricNDSolve[{y''[t] == -w y[t], y[0] == a, y'[0] == b}, y, {t, 0, 10}, {w, a, b}]; In[4]:= fit = FindFit[noisysamples, pfun[w, a, b][t], {{w, .1}, {a, 1}, {b, 1}}, t] Out[4]= {w -> 0.156438, a -> -7.68318, b -> 35.6125} In[5]:= Plot[pfun[w, a, b][t] /. fit, {t, 0, 10}, Epilog -> {Point@noisysamples}] Out[5]= [image] ``` A quadratic model is a better fit: ```wl In[6]:= pfun = y /. ParametricNDSolve[{y''[t] == -g, y[0] == p, y'[0] == v}, y, {t, 0, 10}, {g, p, v}]; In[7]:= fit = FindFit[noisysamples, pfun[g, p, v][t], {{g, 5}, {p, .5}, {v, .5}}, t] Out[7]= {g -> 10.0381, p -> -5.66798, v -> 42.4953} In[8]:= Plot[pfun[g, p, v][t] /. fit, {t, 0, 10}, Epilog -> {Point@noisysamples}] Out[8]= [image] ``` --- Find the parameters that make the solution of a differential equation the best fit to data: ```wl In[1]:= 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: ```wl In[2]:= BoltzmannConstant = Quantity[1.3806504`*^-23, "Joules" / "Kelvins"]; In[3]:= tempK = UnitConvert[temp, "Kelvins"]; fitdata = QuantityMagnitude[Transpose[{time, tempK}]]; Subscript[tempK, 0] = QuantityMagnitude[First[tempK]]; Subscript[tempK, ∞] = QuantityMagnitude[UnitConvert[Quantity[79, "DegreesFahrenheit"], "Kelvins"]]; In[4]:= c = QuantityMagnitude[(2260 18/6.02 10^23 BoltzmannConstant)]; tend = QuantityMagnitude[Last[time]]; In[5]:= 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}]; In[6]:= fit = FindFit[fitdata//N, pfun[K1, K2][t], {{K1, 0.0001}, {K2, 10000}}, t] Out[6]= {K1 -> 0.0000322221, K2 -> 73856.8} In[7]:= Show[ListPlot[fitdata], Plot[Evaluate[pfun[K1, K2][t] /. fit], {t, 0, tend}, PlotStyle -> Red]] Out[7]= [image] ``` ### Properties & Relations (1) Use ``SystemModelParametricSimulate`` to simulate larger hierarchical system models: ```wl In[1]:= pf = SystemModelParametricSimulate[[image], "Inertia2.w", 40, {"Resistor1.R", "SpringDamper1.d"}] Out[1]= ParametricFunction["SystemModelSimulate", "DocumentationExamples.Simulation.HybridMotor", {"Resistor1.R", "SpringDamper1.d"}, WSMLink`WSMExperiment[ "DocumentationExamples.Simulation.HybridMotor", "/tmp/WolframSystemModeler-secavarat-14.3/wsml_ ... framSystemModeler-secavarat-1\ 4.3/WL/e_wsml_HybridMotor_lsim260803_104994_45794_337803023323419989901970908995_1.sim", Missing[], Missing[], Missing[]], "Inertia2.w", {0, 40}, Association[], Sequence @@ WSM`PackageScope`patchOpts[{}, False]] ``` Simulate with two sets of resistor and spring damper parameters: ```wl In[2]:= sims = {pf[0.1, 0.25], pf[0.4, 0.3]}; ``` Compare the resulting angular velocities: ```wl In[3]:= Plot[Evaluate@Comap[sims, t], {t, 0, 40}] Out[3]= [image] ``` ## See Also * [`ParametricNDSolveValue`](https://reference.wolfram.com/language/ref/ParametricNDSolveValue.en.md) * [`ParametricFunction`](https://reference.wolfram.com/language/ref/ParametricFunction.en.md) * [`NDSolve`](https://reference.wolfram.com/language/ref/NDSolve.en.md) * [`DSolve`](https://reference.wolfram.com/language/ref/DSolve.en.md) * [`SystemModelParametricSimulate`](https://reference.wolfram.com/language/ref/SystemModelParametricSimulate.en.md) ## Related Guides * [Equation Solving](https://reference.wolfram.com/language/guide/EquationSolving.en.md) * [Differential Equations](https://reference.wolfram.com/language/guide/DifferentialEquations.en.md) * [Systems Modeling](https://reference.wolfram.com/language/guide/SystemsModeling.en.md) ## History * [Introduced in 2012 (9.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn90.en.md) \| [Updated in 2014 (10.0)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn100.en.md)