Introduction to Numerical Differential-Algebraic Equation Solvers
Differential Equations and Differential-Algebraic Equations (DAEs)
Differential Equations - ODE
Differential equations describe the behavior of a system through rates of change.
Simplest Example:
restart:
ⅆⅆtyt=0
ⅆⅆtyt=0
dsolve,y0=1
yt=1
This can extend to systems with more than one dependent variable.
dsys:=ⅆⅆtxt=yt,ⅆⅆtyt=−xt,x0=0,y0=1
dsys≔ⅆⅆtxt=yt,ⅆⅆtyt=−xt,x0=0,y0=1
dsolvedsys
xt=sint,yt=cost
It is not always possible to find the closed form solution of an ODE problem. This creates a need for solvers that compute a numerical approximation to the solution.
Differential Algebraic Equations - DAE
A DAE is a system that is a combination of differential equations and constraints.
One of the simplest examples is the pendulum system, where the m*g (mass times gravitational constant) term has been replaced by Pi^2.
dsys:=ⅆ2ⅆt2xt=−2λtxt,ⅆ2ⅆt2yt=−2λtyt−π2,xt2+yt2=1
dsys≔ⅆ2ⅆt2xt=−2λtxt,ⅆ2ⅆt2yt=−2λtyt−π2,xt2+yt2=1
To display derivatives and functions using a simpler notation, use the PDEtools[declare] command.
PDEtoolsdeclarext,yt,θt,rt,λt,prime=t,quiet;dsys
x''=−2λx,y''=−2λy−π2,x2+y2=1
The typical solution approach for this problem (as presented in textbooks) is to apply a change of variables to polar co-ordinates, which would proceed as follows.
tr:=xt=rtsinθt,yt=−rtcosθt; nsys≔PDEtoolsdchangetr,dsys,rt,θt
tr≔x=rsinθ,y=−rcosθ
nsys≔r''sinθ+2r'θ'cosθ+rθ''cosθ−rθ'2sinθ=−2λrsinθ,−r''cosθ+2r'θ'sinθ+rθ''sinθ+rθ'2cosθ=2λrcosθ−π2,r2sinθ2+r2cosθ2=1
The last equation implies r=1.
simplifynsys−1,trig
r2=1
Evaluate and remove lambda(t) as follows, giving an ODE in theta(t).
nsys:=nsysrt=1|nsysrt=1;isolatensys1,λt;
nsys≔θ''cosθ−θ'2sinθ=−2λsinθ,θ''sinθ+θ'2cosθ=2λcosθ−π2,sinθ2+cosθ2=1
λ=−θ''cosθ+θ'2sinθ2sinθ
evalnsys2,;
θ''sinθ+θ'2cosθ=−θ''cosθ+θ'2sinθcosθsinθ−π2
pode≔simplifyisolate,ⅆ2ⅆt2θt,trig
pode≔θ''=−sinθπ2
This is the more commonly known pendulum equation.
As an approximation for small displacements, we can approximate sin(theta) ~ theta giving us the simple harmonic oscillator equation:
ⅆ2ⅆt2θt+π2θt=0
θ''+π2θ=0
This is appropriate for the pendulum equation, but, in general, such a coordinate transform may not exist. Hence, there is the need for numerical solution methods for DAEs.
Structure of DAE
DAEs have an interesting property.
dsys
Consider our constraint:
c1≔dsys−1
c1≔x2+y2=1
Differentiate:
c2:=∂∂tc12
c2≔xx'+yy'=0
Differentiate, eliminate, and isolate:
c3:=eval∂∂tc2,dsys1..2;c3:=factorisolatec3,λt
c3≔x'2−2x2λ+y'2+y−2λy−π2=0
c3≔λ=−−x'2−y'2+yπ22x2+y2
One more differentiation and elimination yields an ODE for lambda(t).
d4:=eval∂∂tc3,dsys1..2
d4≔λ'=−4x'λx−2y'−2λy−π2+y'π22x2+y2+−x'2−y'2+yπ22xx'+2yy'2x2+y22
The above process of successive differentiation until we obtain an ODE for each variable is called index reduction. The number of required differentiations is called the gear index. It sometimes corresponds to the difficulty in obtaining a numerical solution of the DAE problem. (A greater index corresponds to a more difficult problem.)
Issues:
We obtained 3 algebraic relations for lambda, x, x', y, y': c1,c2,c3. Initial conditions must satisfy these—as a result, we really only have 2 degrees of freedom in the initial conditions for a problem that would normally have 5. This restricts the initial conditions. They must be checked for consistency. Note: Constraints are typically nonlinear, and cannot normally be solved in terms of the free initial conditions.
The algebraic constraints must hold along the entire solution.
Solution Methods
Method of Index Reduction
As previously shown, this involves successive differentiation of system until the DAE becomes an ODE system.
Advantages
Allows us to obtain all constraints directly
Allows use of a standard ODE solver to obtain the solution
Disadvantages
Does not enforce the constraints, other than at the initial point (resulting in constraint drift)
Example:
IRdsys:=dsys1,dsys2,d4g=π2|dsys1,dsys2,d4g=π2;IRcons≔c1,c2,c3g=π2|c1,c2,c3g=π2
IRdsys≔x''=−2λx,y''=−2λy−π2,λ'=−4x'λx−2y'−2λy−π2+y'π22x2+y2+−x'2−y'2+yπ22xx'+2yy'2x2+y22
IRcons≔x2+y2=1,xx'+yy'=0,λ=−−x'2−y'2+yπ22x2+y2
Find suitable initial conditions:
ICs:=y0=−1,Dx0=110
ICs≔y0=−1,Dx0=110
From the first constraint, we have:
IRcons1yt=−1|IRcons1yt=−1
x2+1=1
ICs:=ICs∪x0=0:
From the second constraint, we have:
eval(convertIRcons2,Dt=0|convertIRcons2,Dt=0,ICs)
−Dy0=0
ICs:=ICs∪Dy0=0:
And finally from the third constraint, we have:
eval(convertIRcons3,Dt=0|convertIRcons3,Dt=0,ICs)
λ0=1200+π22
ICs:=ICs∪
ICs≔λ0=1200+π22,x0=0,y0=−1,Dx0=110,Dy0=0
Now, we can obtain a numerical solution for the index reduced system ignoring the constraints.
(Note: We set maxfun=0 to remove the limit on the number of function evaluations.)
dsn:=dsolveIRdsys∪ICs,numeric,maxfun=0
dsn≔procx_rkf45...end proc
Consider the value at t=1000:
t1000a:=dsn1000
t1000a≔t=1000.,λ=4.72577908209305,x=0.0145466461773924,x'=0.0883768189980660,y=−1.04481889187418,y'=0.00110333926104828
Now consider the constraints at t=1000:
evalf15evalmapa→rhs−lhsa,convertIRcons,D,convertt1000a2..−1,D
−0.0918581217322063,−0.000132796612158576,1.22948907055331×10−7
This indicates that the constraints are not satisfied. This is one of the disadvantages of this approach.
Note: Problems arise when small changes to the problem dictate large changes in the solution at a later time (instability).
The Maple DAE extension solvers (rkf45_dae and rosenbrock_dae) are modifications to ODE solvers, and can operate in this mode:
dsn:=dsolveopdsys∪ICs,numeric,differential=true,projection=false,maxfun=0
dsn≔procx_rkf45_dae...end proc
t1000a2:=dsn1000
t1000a2≔t=1000.,λ=4.71144373570398,x=−0.0213094042255820,x'=0.0742773565891539,y=−1.04753334672762,y'=−0.00164448903842064
evalf15evalmapa→rhs−lhsa,convertIRcons,D,convertt1000a22..−1,D
−0.0977802032148257,−0.000139850889707689,1.07473340449360×10−7
Method of Index Reduction with Removal
This approach is similar to the Method of Index Reduction, but involves the removal of variables that can be solved for. Advantages and disadvantages are the same as those for the basic index reduction method, but this approach usually provides a better solution.
IRRdsys:=evaldsys1,dsys2,c3; IRRcons≔c1,c2
IRRdsys≔x''=−x'2−y'2+yπ2xx2+y2,y''=−x'2−y'2+yπ2yx2+y2−π2
IRRcons≔x2+y2=1,xx'+yy'=0
IRRICs:=removehas,ICs,λ
IRRICs≔x0=0,y0=−1,Dx0=110,Dy0=0
dsn:=dsolveIRRdsys∪IRRICs,numeric,maxfun=0
t1000b:=dsn1000
t1000b≔t=1000.,x=0.0126433163397750,x'=0.0919440494596145,y=−0.998759928311559,y'=0.00116624938372828
Examine the two remaining relevant constraints:
evalf15evalmapa→rhs−lhsa,convertIRRcons,D,convertt1000b2..−1,D
0.00231875215102262,2.32544800803330×10−6
We see that the O(1) error for the general index reduction method has been reduced to an O(10^(-2)) error by using the index reduction with removal method.
Note: The Maple DAE extension solvers can operate in this way, but also have additional error control capabilities. They do not remove the index 1 equation for lambda(t), but instead use it in the numerical integration, and monitor it for possible error. As a result, the solution is more accurate than previously.
dsn:=dsolveopdsys,opICs,numeric,maxfun=0,projection=false
t1000c:=dsn1000
t1000c≔t=1000.,λ=4.94437510816666,x=0.0126469052386010,x'=0.0919391749877821,y=−0.998759726357413,y'=0.00116651428951809
evalf15evalmapa→rhs−lhsa,convertIRRcons,D,convertt1000c2..−1,D
0.00231906479435229,2.32145880546108×10−6
The error in the constraints is now O(10^(-4)).
Method of Index Reduction with Projection
This approach is similar to the Method of Index Reduction, but at each step the solution is projected onto the constraints of the problem. The magnitude of the projection is monitored as part of the error control strategy.
Enforces constraints, removing the constraint drift problem
Requires a modified solver to add a projection step and error control.
Depending on the stability of system, the projection step may introduce errors in the solution.
The Maple rkf45_dae and rosenbrock_dae solvers operate in this mode by default.
dsn:=dsolveopdsys∪ICs,numeric,maxfun=0
t1000d:=dsn1000
t1000d≔t=1000.,λ=4.93952122803426,x=−0.00631755178476623,x'=0.0981561165511285,y=−0.999980043769625,y'=−0.000620117984924986
Now the constraints are satisfied well within the tolerance at any solution point.
evalf15evalmapa→rhs−lhsa,convertIRcons,D,convertt1000d2..−1,D
6.01946159584088×10−10,7.39595684608879×10−10,−6.84481804569259×10−10
Direct Solution Approaches
These methods work directly with the original DAE system, and do not require differentiation of the constraint. They are typically limited to handling DAEs of index 1, 2, or 3.
Does not require symbolic manipulation of the system
Does not indicate constraints on the initial data, and simply fails if these are not satisfied.
Limited to problems with a maximum index.
The core solver is usually more complex than for the other approaches.
Maple has the mebdfi solver (Modified Extended Backward Difference Formula Implicit). This solver has the ability to fully solve implicit DAE systems. For the pendulum example, we get:
dsn:=dsolveopdsys∪ICs,numeric,method=mebdfi,maxfun=0
dsn≔procx_mebdfi...end proc
t1000e:=dsn1000
t1000e≔t=1000.,λ=4.93950094847455379,x=−0.00630802673444276138,x'=0.0979430737998922551,y=−0.999980105346507164,y'=−0.000617845741768820308
evalf15evalmapa→rhs−lhsa,convertIRcons,D,convertt1000e2..−1,D
−2.290093×10−9,−5.921958654×10−9,−3.2170225×10−7
Important: The constraints are all satisfied within tolerance even though they are not explicitly used in the solution process.
Comparison of All Results to a High Accuracy Solution
The pendulum is an excellent example. We have a transformation that converts it to a single nonlinear ODE in theta, which we can solve to high accuracy using our gear method,. This allows us to compare our DAE solutions.
pode
θ''=−sinθπ2
tr0:=subsrt=1,tr
tr0≔x=sinθ,y=−cosθ
tr1:=∂∂ttr0
tr1≔x'=θ'cosθ,y'=θ'sinθ
evalconverttr0∪tr1,Dt=0|converttr0∪tr1,Dt=0,ICs
−1=−cosθ0,0=Dθ0sinθ0,0=sinθ0,110=Dθ0cosθ0
From the last two we get theta(0)=0, and so:
eval,theta0=0;
−1=−1,0=0,110=Dθ0
Now obtain a high accuracy solution to this problem using the 'gear' method.
Digits:=truncevalhfDigits:dsn:=dsolveⅆ2ⅆt2θt+π2sinθt=0,θ0=0,Dθ0=110,numeric,method=gear,abserr=1.10-13,relerr=1.10-13:tt≔time:t1000≔dsn1000;time−tt
t1000≔t=1000.,θ=−0.00629176884678665,θ'=0.0980270351645652
2.203
Convert to x, y coordinates, and obtain the lambda value:
evaltr0 ∪ tr1,t10002..−1
x'=0.0980250949044482,y'=−0.000616759376763182,x=−0.00629172733550273,y=−0.999980206887684
ans:=∪evalfevalc3,
ans≔x'=0.0980250949044482,y'=−0.000616759376763182,λ=4.93950917526204,x=−0.00629172733550273,y=−0.999980206887684
Compare with the basic index reduction solution, which is expected to be poor:
mapa→rhs−lhsa,evalt1000a2..−1,ans
−0.213730093168983,0.0208383735128951,−0.00964827590638218,−0.0448386849864990,0.00172009863781146
Compare with index reduction with removal, which is expected to be moderately poor:
mapa→rhs−lhsa,evalt1000b2..−1,ans
0.0189350436752778,−0.00608104544483369,0.00122027857612550,0.00178300876049147
Compare with index reduction with removal and additional error correction (extension methods with projection=false):
mapa→rhs−lhsa,evalt1000c2..−1,ans
0.00486593290462611,0.0189386325741037,−0.00608591991666610,0.00122048053027168,0.00178327366628127
Compare with index reduction with projection (extension method in default mode):
mapa→rhs−lhsa,evalt1000d2..−1,ans
0.0000120527722291541,−0.0000258244492635012,0.000131021646680274,1.63118059393064×10−7,−3.35860816180408×10−6
Compare with direct method (mebdfi method):
mapa→rhs−lhsa,evalt1000e2..−1,ans
−8.22678748502170×10−6,−0.0000162993989400338,−0.0000820211045558933,1.01541177133235×10−7,−1.08636500563820×10−6
In summary, we see that the Maple default rkf45_dae method and mebdfi method give the best agreement with the solution.
Other Problems
The Car Axis System
This multi-body system models the behavior of car suspension on a bumpy road. It consists of 4 second order ODEs and two index-3 DAEs.
restart;PDEtoolsdeclareprime=t,quiet:
Auxiliary variables:
yb:=sin10t10:xb≔1−yb2:ll≔xlt2+ylt2:lr≔xrt−xb2+yrt−yb2:
DAE system:
eqs:=ⅆ2ⅆt2xlt2000=12ll−1⋅xlt+lam1t⋅xb+2⋅lam2t⋅xlt−xrt,ⅆ2ⅆt2ylt2000=12ll−1⋅ylt+lam1t⋅yb+2⋅lam2t⋅ylt−yrt−12000,ⅆ2ⅆt2xrt2000=12lr−1⋅xrt−xb−2⋅lam2t⋅xlt−xrt,ⅆ2ⅆt2yrt2000=12lr−1⋅yrt−yb−2⋅lam2t⋅ylt−yrt−12000,0=xltxb+yltyb,0=xlt−xrt2+ylt−yrt2−1
eqs≔0=xlt−sin10t2+10010+yltsin10t10,0=xlt−xrt2+ylt−yrt2−1,xl''2000=12xlt2+ylt2−1xlt+lam1t−sin10t2+10010+2lam2txlt−xrt,xr''2000=12xrt−−sin10t2+100102+yrt−sin10t102−1xrt−−sin10t2+10010−2lam2txlt−xrt,yl''2000=12xlt2+ylt2−1ylt+lam1tsin10t10+2lam2tylt−yrt−12000,yr''2000=12xrt−−sin10t2+100102+yrt−sin10t102−1yrt−sin10t10−2lam2tylt−yrt−12000
Initial conditions and variables:
ini:=xl0=0,Dxl0=−12,yl0=12,Dyl0=0,xr0=1,Dxr0=−12,yr0=12,Dyr0=0,lam10=0,lam20=0;vrs:=xlt,ylt,xrt,yrt,lam1t,lam2t:
ini≔lam10=0,lam20=0,xl0=0,xr0=1,yl0=12,yr0=12,Dxl0=−12,Dxr0=−12,Dyl0=0,Dyr0=0
Exact solution (solution computed via gear with parameter removal and 1e-22 accuracy):
yans:=0.0493455784275240921227327851163,−0.0770583684035920842786873048316,0.496989460230008106764570177300,0.00744686659206841657568543179878,1.04174252488542611518584788701,0.0175568157535417366330123387591,0.373911027265365819359535268815,0.770341043779601063189353416447,−0.00473688659085332651484983657657,−0.00110468033125956583962543762007:
tt:=time:ds1:=dsolveopeqs,opini,vrs,numeric;time−tt
ds1≔procx_rkf45_dae...end proc
1.125
Solution with default tolerances:
tt:=time:dsol1:=ds13;time−tt;
dsol1≔t=3.,xlt=0.0493455779014815,xl'=−0.0770590607488527,ylt=0.496989454931791,yl'=0.00743988527830008,xrt=1.04174242180064,xr'=0.0175614035187829,yrt=0.373910194846429,yr'=0.770371407475209,lam1t=−0.00473687080377615,lam2t=−0.00110467301491489
0.110
Compare:
evalf16maprhs,dsol12..−1−yans
−5.2604258×10−10,−6.9234526065×10−7,−5.2982167×10−9,−6.981313768339×10−6,−1.03084791×10−7,4.58776524117×10−6,−8.324189364×10−7,0.0000303636956078,1.5787077181×10−8,7.316344672×10−9
Solution with tight tolerances:
tt:=time:ds2:=dsolveopeqs,opini,vrs,numeric,abserr=1.10-11,relerr=1.10-11,maxfun=0;time−tt;tt≔time:dsol2:=ds23;time−tt; evalf16maprhs,dsol22..−1−yans
ds2≔procx_rkf45_dae...end proc
0.047
dsol2≔t=3.,xlt=0.0493455784275599,xl'=−0.0770583684067605,ylt=0.496989460230369,yl'=0.00744686656072501,xrt=1.04174252488285,xr'=0.0175568158504510,yrt=0.373911027244694,yr'=0.770341044422276,lam1t=−0.00473688659029054,lam2t=−0.00110468033098021
0.750
3.582×10−14,−3.16846×10−12,3.606×10−13,−3.1343407×10−11,−2.575×10−12,9.690928×10−11,−2.06719×10−11,6.426745×10−10,5.62785×10−13,2.79361×10−13
Trajectory of x1,y1:
plotsodeplotds1,xlt,ylt,0..3,numpoints=500
Trajectory of x2,y2:
plotsodeplotds1,xrt,yrt,0..3,numpoints=500
The Double Pendulum
The double pendulum is a pendulum with an additional weight attached to the first pendulum bob by a string. It is a system with 4 second order ODEs and two index-3 DAEs. (x1(t),y1(t)) are the coordinates of the first pendulum bob, and (x2(t),y2(t)) are the coordinates of the second, which is attached to the first.
restart; PDEtoolsdeclareprime=t,quiet:
Set masses, lengths of strings, and g:
m1≔1:l1≔1: m2≔1: l2 ≔ 12: g≔9.8:
DAE system and initial conditions:
sys:=m1⋅ⅆ2ⅆt2x1t+2⋅λ1t⋅x1t+2⋅λ2t⋅x1t−x2t,m1⋅ⅆ2ⅆt2y1t+m1⋅g+2⋅λ1t⋅y1t+2⋅λ2t⋅y1t−y2t,m2⋅ⅆ2ⅆt2x2t−2⋅λ2t⋅x1t−x2t,m2⋅ⅆ2ⅆt2y2t+m2⋅g−2⋅λ2t⋅y1t−y2t,x1t2+y1t2−l12,x1t−x2t2+y1t−y2t2−l22;ics:=x10=0,Dx10=−52,y10=−l1,Dy10=0,x20=0,Dx20=3,y20=−l1−l2,Dy20=0:
sys≔x2''−2λ2tx1t−x2t,x1t−x2t2+y1t−y2t2−14,x1t2+y1t2−1,x1''+2λ1tx1t+2λ2tx1t−x2t,y2''+9.8−2λ2ty1t−y2t,y1''+9.8+2λ1ty1t+2λ2ty1t−y2t
Obtain the trajectory for the second mass:
dsn:=dsolvesys∪ics,numeric;plotsodeplotdsn,x2t,y2t,0..5,numpoints=1000
Construct an animation of the system, where tf is the final time and nfr is the number of frames:
tf≔5: nfr≔201: pl≔NULL: ps≔0.02:foritonfrdodsv:=dsn5i−1nfr−1:x1v,x2v,y1v,y2v:=seqrhsdsvj,j=4,6,8,10:pl:=pl,CURVES0,0,x1v,y1v,x2v,y2v,COLORRGB,1,0,0,POLYGONSx1v+ps,y1v+ps,x1v+ps,y1v−ps,x1v−ps,y1v−ps,x1v−ps,y1v+ps,x1v+ps,y1v+ps,x2v+ps,y2v+ps,x2v+ps,y2v−ps,x2v−ps,y2v−ps,x2v−ps,y2v+ps,x2v+ps,y2v+ps,COLORRGB,0,0,1:end do:
PLOTANIMATEpl,AXESSTYLENONE
Illustration of Constraint Drift
Consider the simple pendulum system and the worst rkf45_dae solution (differential=true, projection=false):
restart; PDEtoolsdeclarext,yt,λt,prime=t,quiet:dsys:=ⅆ2ⅆt2xt=−2⋅λt⋅xt,ⅆ2ⅆt2yt=−2⋅λt⋅yt−π2,xt2+yt2=1,y0=−1,Dx0=110,Dy0=0,x0=0,λ0=1200+1π22
dsys≔x''=−2λx,y''=−2λy−π2,x2+y2=1,y0=−1,Dx0=110,Dy0=0,x0=0,λ0=1200+π22
dsn:=dsolvedsys,numeric,differential=true,projection=false,maxfun=0
Since this is a pendulum, the trajectory (x(t),y(t)) should trace a portion of a circle repeatedly. Constraint drift causes a deviation this behavior.
plotsodeplotdsn,xt,yt,0..100,numpoints=2000,frames=100
The default solution mode works much better:
dsn:=dsolvedsys,numeric,maxfun=0
Return to Index for Example Worksheets
Download Help Document
