Parametric Differential Equations

Nina Dokeva

Wolfram Research

Introduction

Motivation

Differential system models in many application areas usually depend on a variety of parameters.

Goals

to be able to compute sensitivities, which are derivatives of the solutions with respect to parameters. This helps for:

analyzing the influence of each parameter on the model behavior

finding parameters that satisfy certain requirements for the solution

to be able to repeatedly solve numerically a differential system for various parameters in an efficient way

integration with other numerical Mathematica functions

Scope

ParametricNDSolve, just like NDSolve, works with:

ordinary differential equations

partial differential equations

differential algebraic equations

delay differential equations

Example: Mackey-Glass model

Consider the Mackey–Glass equation for the density P[t] of mature circulating white blood cells.The parameter τ is the delay between production and release into the bloodstream of the cells.

ParametricNDSolve[{P'[t] == β (θ^n) P[t - τ]/ ((θ^n) + ((P[t - τ])^n)) - γ P[t], P[t /; t ≤ 0]==A},P[t],{t,0,1000},{β ,θ ,γ ,n,τ,A}]

pfun1=ParametricNDSolveValue[{P'[t] == β (θ^n) P[t - τ]/ ((θ^n) + ((P[t - τ])^n)) - γ P[t], P[t /; t ≤ 0]==A},P[t],{t,0,1000},{β ,θ ,γ ,n,τ,A}]

sol1=pfun1[0.2,  8, .1 , 10,15,0.1]

Plot[sol1,{t,0,1000}, AxesLabel->{t, P[t]}]

ParametricNDSolveValue can take for an output expression of the solution function(s):

pfun2 = ParametricNDSolveValue[{P'[t] == β (θ^n) P[t - τ]/ ((θ^n) + ((P[t - τ])^n)) - γ P[t], P[t /; t ≤ 0]==A},{P[t], P[t-τ]},{t,-100,1000},{β ,θ ,γ ,n,τ,A} ]

sol2=pfun2[0.2,  8, .1 , 10,15,0.1]

ParametricPlot[sol2,{t,0,1000},AxesLabel->{P[t], P[t-τ]}]

pfun3 = ParametricNDSolveValue[{P'[t] == β (θ^n) P[t - τ]/ ((θ^n) + ((P[t - τ])^n)) - γ P[t], P[t /; t ≤ 0]==A},{P[t], P[t-τ]},{t,-100,1000},{β ,θ ,γ ,n,τ,A} , Method->{"ParametricSensitivity"->None}];

Manipulate[ParametricPlot[Evaluate[pfun3[β ,θ ,γ ,n,τ,A]],{t,0,1000}, AxesLabel->{P[t], P[t-τ]}],{{β,0.2},0.1, 1},{{θ, 8},2,20},{{γ,.1},0.02,0.15},{{n,10},5,100},{{τ,15},10,30},{{A,0.1},0,1},ControlPlacement->Left]

Scope of the parameters and the output

Scope of the parameters

Parameters can appear almost anywhere, even in the options:

err = ParametricNDSolveValue[{x'[t] == c x[t], x[0]==a},x[1] -a Exp[c], {t,0,T}, {a,c,T,wp}, Method->"Extrapolation", WorkingPrecision->wp]

Table[err[1,1,1,2^k], {k,4,6}]

Scope of the output expression for ParametricNDSolveValue

When using ParametricNDSolveValue the target output need not be the list of dependent variables but can be some expression depending on them.
In the case of multiple solutions, ParametricNDSolveValue will return only one branch.

Sample the solution at discrete times

eqns={y''[t]+y[t]==3a Sin[y[t]],y[0]==y'[0]==1};
pfun=ParametricNDSolveValue[eqns,Table[y[t],{t,0,5,1}],{t,0,5},{a}]

pfun[1]

Return the integral of the solution function

eqns={y''[t]+y[t]==3a Sin[y[t]],y[0]==y'[0]==1};
pfun=ParametricNDSolveValue[eqns,(∫_0^5)y[s]\[DifferentialD]s,{t,0,5},{a}]

Plot[pfun[a],{a,-2,2}, AxesLabel->{a,(∫_0^5)y[s,a]\[DifferentialD]s}]

Compute the total energy of a damped driven harmonic oscillator starting from rest

eqns = {x''[t] + γ x'[t] + x[t] == ϵ Cos[ω t], x[0] == x'[0] == 0};
pf2 = ParametricNDSolveValue[eqns, (1/2)(x[#]^2+x'[#]^2)&, {t, T}, {γ, ϵ, ω, T}]

Plot[Evaluate[Table[pf2[.1, .2, ω,10][t],{ω,.6,1.4,.1}]],{t,0,10}, AxesLabel->{t,(x[t]^2+x'[t]^2)/2}]

Sensitivity solutions

The main purpose of ParametricNDSolve is to be able to compute sensitivities, which are derivatives of the solutions with respect to parameters. This is achieved by automatically setting up additional equations, called (forward) sensitivity equations, and solving the augmented system.

For ODEs:

Original system:

((∂Y(t,P))/(∂t))= f(t,Y(t,P), P),
Y((Subscript[t, 0]),P) = (Subscript[Y, 0])(P),    where   Y=((Subscript[y, 1]),. . . ,(Subscript[y, N])),    P=((Subscript[p, 1]),. . . ,(Subscript[p, L]))

Definition of sensitivities:

(Subscript[s, i ])(t,P)= ((∂Y(t,P))/(∂(Subscript[p, i]))) =  (((Subscript[y, 1])/(∂(Subscript[p, i]))) ,. . . ,((Subscript[y, N])/(∂(Subscript[p, i])))),    for    i=(1,. . . ,L)

Sensitivity equations:

((∂(Subscript[s, i ])(t,P))/(∂t))= ((∂ f(t,Y, P))/(∂Y)) (Subscript[s, i ])(t,P) + ((∂ f(t,Y, P))/(Subscript[p, i])),
(Subscript[s, i])((Subscript[t, 0]),P) =((∂(Subscript[Y, 0])(P))/(∂(Subscript[p, i]))),   for  i=(1,. . . ,L)

Derivation (Differentiation + Clairaut’s Theorem + Chain Rule):

((∂(Subscript[s, i ])(t,P))/(∂t)) =(∂/(∂t))((∂Y(t,P))/(∂(Subscript[p, i])))=(∂/(∂(Subscript[p, i])))((∂Y(t,P))/(∂t))= (∂/(∂(Subscript[p, i])))f(t,Y(t,P), P) =
= ((∂f)/(∂Y)) ((∂Y)/(∂(Subscript[p, i])))+ ((∂f)/(∂(Subscript[p, i])))= ((∂f)/(∂Y))  (Subscript[s, i ])+ ((∂f)/(∂(Subscript[p, i]))),
(Subscript[s, i ])((Subscript[t, 0]),P)  = ((∂Y((Subscript[t, 0]),P))/(∂(Subscript[p, i])))= ((∂(Subscript[Y, 0])(P))/(∂(Subscript[p, i]))),      i=(1,. . . ,L)

Dimension of the augmented system:

N (L+1)

Note that the equations for the sensitivity with respect to each parameter have the same Jacobian ((∂f)/(∂Y)) as does the original system, so it is possible to reuse the Jacobian. Unless the method is otherwise specified, NDSolve uses the codes CVODES and IDAS from the Sundials library for ODEs and DAEs respectively.

Numerical Sensitivity Analysis examples

Van der Pol oscillator

Van der Pol oscillator equation (known to appear in electrical engineering, physics, biology, seismology):

x''[t]-μ(1-x[t]^2)x'[t]+x[t]==0

ParametricNDSolve will turn this into a first order system of equations and add sensitivity equations to the system:

y1'[t] = y2[t]
y2'[t] = μ( 1- y1[t]^2)y2[t] - y1[t]

s11'[t]=s12[t]
s12'[t]=(−2μ y1[t]y2[t]−1)s11[t]+μ(1−y1[t]^2)s12[t]+(1−y1[t]^2)y2[t]

pfun=ParametricNDSolveValue[{x''[t]-μ(1-x[t]^2)x'[t]+x[t]==0,x[0]==2,x'[0]==0},x,{t,0,20},{μ}]

Find sensitivities at μ = 2.341:

With[{μ=2.341},GraphicsRow[{Plot[pfun[μ][t],{t,0,20}, ImageSize->270,AxesLabel->{t,x[t,μ]}],Plot[pfun'[μ][t],{t,0,20}, PlotRange->All,AxesLabel->{t,Derivative[0,1][x][t,μ]}],Plot[pfun''[μ][t],{t,0,20}, PlotRange->All,AxesLabel->{t,Derivative[0,2][x][t,μ]}]}]]

GraphicsRow[{Plot[pfun[μ][20],{μ,2.3,2.4}, ImageSize->280,AxesLabel->{μ,x[20,μ]}],Plot[Evaluate[pfun[2.341][t]+.5{0,-1,1}pfun'[2.341][t]],{t,0,20}, Filling->{2->{3}}],Plot[Evaluate[Table[pfun[μ][t],{μ,{2.341,2.5,2.1}}]],{t,0,20}]}]

Duffing equation

Models certain damped and driven oscillators.

xpfun =ParametricNDSolveValue[{x''[t]+δ x'[t]+α x[t]+βx[t]^3==γ Cos[ω t],x[0]==a,x'[0]==b},x,{t,0,100},{α,β,γ,δ, ω ,a,b}]

Parameters:
  δ controls the size of the damping (friction).
  α controls the size of the restoring force.
  β controls the amount of non-linearity in the restoring force.If β = 0 the Duffing equation describes a damped and driven simple harmonic oscillator.
  γ controls the amplitude of the periodic driving force.If γ = 0 we have a system without driving force.
  ω controls the frequency of the periodic driving force.

Sensitivity to initial conditions:

Sensitivity to b when α = -1:

GraphicsRow[{Plot[Evaluate[xpfun[-1,1,0.3,0.15, 1 ,-1,1][t]+.5{0,-1,1}(D[xpfun[-1,1,0.3,0.15, 1 ,-1,b],b]/.b->1) [t]], {t,0,50}, Filling->{2->{3}}, PlotRange->All,ImageSize->350],Plot[Evaluate[Table[xpfun[-1,1,0.3,0.15, 1 ,-1,b][t],{b,{1,1.001}}]], {t, 0,50}] }]

Sensitivity to b when α = 1:

GraphicsRow[{Plot[Evaluate[xpfun[1,1,0.3,0.15, 1 ,-1,1][t]+.5{0,-1,1}(D[xpfun[1,1,0.3,0.15, 1 ,-1,b],b]/.b->1) [t]], {t,0,50}, Filling->{2->{3}}, PlotRange->All,ImageSize->350],Plot[Evaluate[Table[xpfun[1,1,0.3,0.15, 1 ,-1,b][t],{b,{1,1.001}}]], {t, 0,50}] }]

Sensitivity to the size of the damping δ:

Sensitivity to δ when α = -1:

GraphicsRow[{Plot[Evaluate[xpfun[-1,1,0.3,0.15, 1 ,-1,1][t]+.5{0,-1,1}(D[xpfun[-1,1,0.3,δ, 1 ,-1,1],δ]/.δ->0.15) [t]], {t,0,50}, Filling->{2->{3}}, PlotRange->All,ImageSize->350],Plot[Evaluate[Table[xpfun[-1,1,0.3,δ, 1 ,-1,1][t],{δ,{0.15,0.151}}]], {t, 0,50}] }]

Sensitivity to δ when α = 1:

GraphicsRow[{Plot[Evaluate[xpfun[1,1,0.3,0.15, 1 ,-1,1][t]+.5{0,-1,1}(D[xpfun[1,1,0.3,δ, 1 ,-1,1],δ]/.δ->0.15) [t]], {t,0,100}, Filling->{2->{3}}, PlotRange->All,ImageSize->350],Plot[Evaluate[Table[xpfun[1,1,0.3,δ, 1 ,-1,1][t],{δ,{0.15,0.151}}]], {t, 0,100}] }]

Mackey-Glass equation

pfun = ParametricNDSolveValue[{P'[t] == β (θ^n) P[t - τ]/ ((θ^n) + ((P[t - τ])^n)) - γ P[t], P[t /; t ≤ 0]==A},P[t],{t,-100,1000},{β ,θ ,γ ,n,τ,A} ]

Sensitivity with respect to the initial condition A:

GraphicsRow[{Plot[Evaluate[{pfun[0.2, 8, .1 , 10, 15, 0.1],pfun[0.2, 8, .1 , 10, 15, 0.12]}],{t,0,500},PlotRange->All],Plot[Evaluate[D[pfun[0.2,8, .1 , 10,15,A],A]/.A->0.1],{t,0,500},PlotRange->All, ImageSize->350]}]

Sensitivity with respect to τ (white blood cells maturation time):

At τ = 15:

GraphicsRow[{Plot[Evaluate[{pfun[0.2, 8, .1 , 10, 15, 0.1],pfun[0.2, 8, .1 , 10, 15.2, 0.1]}],{t,0,500},PlotRange->All],Plot[Evaluate[D[pfun[0.2,8, .1 , 10,τ,.1],τ]/.τ->15],{t,0,500},PlotRange->All, ImageSize->350]}]

At τ = 4.4:

GraphicsRow[{Plot[Evaluate[Table[pfun[0.2, 8, .1 , 10, τ, .1],{τ,{4.4,4.5}}]],{t,0,500}, PlotRange->All],Plot[Evaluate[D[pfun[0.2,8, .1 , 10,τ,.1],τ]/.τ->4.4],{t,0,500}, PlotRange->All, ImageSize->350]}]

FindRoot with ParametricFunction

For a given parametric equation, find parameter values that satisfy certain equality constraints for the solution. 
FindRoot is using the Jacobian derived from numerically computing the solution to the sensitivity equations.

For the Van der Pol oscillator find the value of μ, such that x[20]=0:

pfun=ParametricNDSolveValue[{x''[t]-μ(1-x[t]^2)x'[t]+x[t]==0,x[0]==2,x'[0]==0},x,{t,0,20},{μ}]

sol=FindRoot[pfun[μ][20]==0,{μ,3}]

Plot[pfun[μ][t]/.sol,{t,0,20}]

For the Van der Pol oscillator find the value of μ, such that ((∂x(20, μ))/(∂μ))=0:

pfun=ParametricNDSolveValue[{x''[t]-μ(1-x[t]^2)x'[t]+x[t]==0,x[0]==2,x'[0]==0},x,{t,0,20},{μ}]

sol=FindRoot[pfun'[μ][20]==0,{μ,2}]

Plot[pfun'[μ][t]/.sol,{t,0,20}, PlotRange->All]

Find initial condition X'[0] = (Subscript[X, 0]) to satisfy a boundary condition X[T]=C:

We have x[0] = 0, y[0] = 1 and  x[T] = 5, y[T] = 2, want to find {x'[0], y'[0]}.

Assumptions:
We will stretch the sling all the way to the ground, a meters behind the sling base. 
We shoot an object of weight 1kg and the square root of the elasticity constant of the sling is 5.
The initial velocity is 5*{a,1}, need to find a:

pfun = ParametricNDSolve[{y''[t]==-9.81, x''[t]==0, x[0]==0,y[0]==1,x'[0]==5a,y'[0]==5 1},{x,y},{t,0,T},{a,T}]

sol=FindRoot[{x[a,T][T]-5==0,y[a,T][T]-2==0}/.pfun,{{a,2},{T,3}}]

image=;

ParametricPlot[Evaluate[{x[a,T][t],y[a,T][t]}/.pfun/.sol],{t,0,T/.sol}, AxesOrigin->{0,0}, Epilog->Inset[image,{x[a,T][T],y[a,T][T]}/.pfun/.sol], PlotRange->Full, PlotRangePadding->Scaled[0.05]]

Fitting data to a differential equation model

Find the parameter μ that gives the best fit for a given data to the Van der Pol equation. FindFit takes advantage of the sensitivity solutions to speed up the minimization process.

Create data to fit:

peqns={x''[t]-μ(1-x[t]^2)x'[t]+x[t]==0,x[0]==2, x'[0]==0};
data=NDSolveValue[peqns/.μ->1.23,Table[{t,x[t]+RandomReal[.1]},{t,0,10,.1}],{t,10}];

lp=ListPlot[data,PlotRange->All]

Create ParametricFunction and use FindFit:

pfun=ParametricNDSolveValue[peqns,x,{t,0,10},{μ}]

fit = FindFit[data, pfun[μ][t],{{μ,2}},t]//Quiet

Show[Plot[Evaluate[pfun[μ][t] /. fit], {t, 0, 10}], ListPlot[data]]

Preprocessing and caching

Preprocessing

When the parameters only appear in the equations, equation processing is only done once. A function is set up so that the preprocessed equations can take different parameter values.

Solve Duffing equation with 1000 sets of parameters and a single preprocessing:

peqns= {x''[t]+δ x'[t]+α x[t]+β x[t]^3==γ Cos[ω t],x[0]==a,x'[0]==b};

pfun =ParametricNDSolveValue[peqns,x,{t,0,10},{α,β,γ,δ, ω ,a,b}]

pardata = Table[RandomReal[],{1000},{7}];

AbsoluteTiming[Apply[pfun,pardata, {1}];]

With the preprocessing repeated every time:

fun = Function[{α,β,γ,δ,ω,a,b}, NDSolveValue[eqns,x,{t,10}]]/.eqns->peqns

AbsoluteTiming[Apply[fun, pardata, {1}];]

Caching

ParametricFunction automatically caches the most recent solution computed.

In the fitting example from the previous slide:

peqns={x''[t]-μ(1-x[t]^2)x'[t]+x[t]==0,x[0]==2, x'[0]==0};
data=NDSolveValue[peqns/.μ->1.23,Table[{t,x[t]+RandomReal[.1]},{t,0,10,.1}],{t,10}];

pfun=ParametricNDSolveValue[peqns,x,{t,0,10},{μ}];
AbsoluteTiming[fit = FindFit[data, pfun[μ][t],{{μ,2}},t]//Quiet]

With the caching turned off:

pfun=ParametricNDSolveValue[peqns,x,{t,0,10},{μ},Method->{"ParametricCaching"->None}];
AbsoluteTiming[fit = FindFit[data, pfun[μ][t],{{μ,2}},t]//Quiet]

Parametric dependence in a PDE

Wave equation:

pf = ParametricNDSolveValue[{D[u[t,x],t,t] == c^2 D[u[t,x],x,x], u[0,x] == Exp[-(a x)^2],Derivative[1,0][u][0,x] ==0,  u[t, -10] == u[t,10]}, u, {t,0,10},{x,-10,10},{a, c}]

Plot the solution:

{Plot3D[Evaluate[pf[1,2][t,x]],{t,0,10},{x,-10,10}, PlotRange->All,PlotPoints->100,ImageSize->300],DensityPlot[Evaluate[pf[1,2][t,x]], {t,0,10},{x,-10,10},PlotRange->All, PlotPoints->100,ImageSize->300]}//Quiet

Plot the sensitivity of the solution to the wave speed:

{DensityPlot[Evaluate[Derivative[0,1][pf][1,2][t,x]], {t,0,10},{x,-10,10},PlotRange->All,PlotPoints->100, ImageSize->230],Plot[Evaluate[{pf[1,2.][t,x],pf[1,2.1][t,x]}/.t->2],{x,-10,10}, PlotRange->All, ImageSize->230], Plot[Evaluate[{pf[1,2.][t,x],pf[1,2.1][t,x]}/.t->10],{x,-10,10}, PlotRange->All, ImageSize->230]}//Quiet

Manipulate[Plot[Evaluate[pf[1,2.1][t,x]],{x,-10,10}, PlotRange->All]//Quiet,{t,0,10}]

Plot the sensitivity of the solution to the width of the initial condition (when a increases, the wave becomes more narrow):

{DensityPlot[Evaluate[Derivative[1,0][pf][1,2][t,x]], {t,0,10},{x,-10,10}, PlotRange->All, PlotPoints->100,ImageSize->300],Plot[Evaluate@Table[Exp[-(((a x))^2)],{a,{0.5,1,2}}],{x,-10,10}, PlotRange->All,ImageSize->300]}//Quiet

Parametric equations with events

Forced impact oscillator - a forced oscillator that reverses direction upon impact with a wall at x[t]=0:

pfun=ParametricNDSolveValue[{x''[t]+x[t]==Cos[ω t],x[0]==0,x'[0]==0,WhenEvent[x[t]==0,x'[t]->-r x'[t]]},x,{t,0,300},{ω,r}];

Compare with a forced oscillator without impacts:

pfun1=ParametricNDSolveValue[{x''[t]+x[t]==Cos[ω t],x[0]==0,x'[0]==0},x,{t,0,300},{ω,r}];

Plot[Evaluate[{pfun[2.76,0.95][t],pfun1[2.76,0.95][t]}],{t,0,15}, AxesLabel->{t,x[t]}]

Compare the behavior for different parameters:

GraphicsGrid[{Apply[ParametricPlot[Evaluate[{pfun[#1,#2] [t],pfun[#1,#2]'[t]}],{t,200,300}, AxesLabel->{x[t],x'[t]},PlotPoints->200,PlotRange->{{-.1,1},{-1.1,1}},ImageSize->120]&,{{3,0.95},{2.76,0.52},{2.76,0.95},{2.76,0.82},{2.9,0.95}},{1}]}]

Conclusions

ParametricNDSolve[Value] serves to make solving and analyzing differential equation with parameters easier and more efficient and has the following features:

Parameters can appear almost anywhere in the input

The target output of ParametricNDSolveValue can be a general expression of the solution functions

Preprocessing of the equations is done only once when possible

Caching of the most recent solution is done automatically

Derivatives with respect to parameters are computed by automatically setting up sensitivity equations that augment the original system

ParametricNDSolve[Value] works well with other numerical functions that can take advantage of the sensitivity computations.

Thank you!