Modeling Physics with Differential-Algebraic Equations Lecture 3 - - PowerPoint PPT Presentation

Modeling Physics with Differential-Algebraic Equations

Lecture 3

Numerical Integration of DAEs

COMASIC (M2) January 7th, 2019

Khalil Ghorbal khalil.ghorbal@inria.fr

• K. Ghorbal (INRIA)

1 COMASIC M2 1 / 22

Summary of lecture 2

Index Reduction

• Given a DAE F(x, ˙

x, t), we have seen how to perform a structural analysis to numerically compute ˙ x function of x. The structural nonsingularity ensures that, generically, one can perform the computation following the block order suggested by the BLT

• decomposition. Thus, one is able to compute the numerical values of

the derivatives given a consistent state of the system and carry on with a standard numerical integration.

• Note: the structural index is not always equal to the differentiation

index.

• K. Ghorbal (INRIA)

2 COMASIC M2 2 / 22

Summary of lecture 2

Index Reduction

• Given a DAE F(x, ˙

x, t), we have seen how to perform a structural analysis to numerically compute ˙ x function of x. The structural nonsingularity ensures that, generically, one can perform the computation following the block order suggested by the BLT

• decomposition. Thus, one is able to compute the numerical values of

the derivatives given a consistent state of the system and carry on with a standard numerical integration.

• Note: the structural index is not always equal to the differentiation

index.

• K. Ghorbal (INRIA)

2 COMASIC M2 2 / 22

Structural Analysis (Complement)

Consider the following DAE ˙ z − ˙ xy − x ˙ y + 2x + y − 3 = 0 z − xy = 0 x + y − 2 = 0

• Pantelides: structural index 1
• Differentiation index is 0 (simple linear system)
• (Hidden) Cancellation problems are undecidable in general
• K. Ghorbal (INRIA)

3 COMASIC M2 3 / 22

Structural Analysis (Complement)

Consider the following DAE ˙ z − ˙ xy − x ˙ y + 2x + y − 3 = 0 z − xy = 0 x + y − 2 = 0

• Pantelides: structural index 1
• Differentiation index is 0 (simple linear system)
• (Hidden) Cancellation problems are undecidable in general
• K. Ghorbal (INRIA)

3 COMASIC M2 3 / 22

Structural Analysis (Complement)

Consider the following DAE ˙ z − ˙ xy − x ˙ y + 2x + y − 3 = 0 z − xy = 0 x + y − 2 = 0

• Pantelides: structural index 1
• Differentiation index is 0 (simple linear system)
• (Hidden) Cancellation problems are undecidable in general
• K. Ghorbal (INRIA)

3 COMASIC M2 3 / 22

Structural Analysis (Complement)

Consider the following DAE ˙ z − ˙ xy − x ˙ y + 2x + y − 3 = 0 z − xy = 0 x + y − 2 = 0

• Pantelides: structural index 1
• Differentiation index is 0 (simple linear system)
• (Hidden) Cancellation problems are undecidable in general
• K. Ghorbal (INRIA)

3 COMASIC M2 3 / 22

Semi-Explicit DAE

Index reduction transforms a fully implicit DAE to a semi-explicit DAE ˙ x = f (x, y, t) = g(x, y, t)

Integration Schemes

• Backward Differentiation Formula (BDF)
• Orthogonal Collocation
• K. Ghorbal (INRIA)

4 COMASIC M2 4 / 22

Outline

1 BDF Method 2 Collocation (Overview) 3 Brief Introduction to Modelica

• K. Ghorbal (INRIA)

4 COMASIC M2 4 / 22

Backward Differentiation Formula (BDF)

Fixed Step Size

• Implicit: next value not explicitly given.
• Linear multistep: the next value is linearly related to the immediate

previous (backward) values (eventually more than one).

• For a fixed step size δ > 0, let tn = t0 + nδ, and xs the approximate
• f the exact x(ts).
• The general Backward Differentiation Formula:

˙ xn = linear combination of xn, xn−1, . . . , x0

• For a Cauchy problem ˙

x = f (x, t), x(t0) = x0, one obtains an implicit equation for xn: xn =

q

• n=1

an x1−n + δ bq f (xn, tn) where the an and bq depends only on the order q.

• E.g., q = 1 (BDF1): xn = xn−1 + δ f (xn, tn)

(a.k.a. Backward Euler)

• K. Ghorbal (INRIA)

5 COMASIC M2 5 / 22

Backward Differentiation Formula (BDF)

Fixed Step Size

• Implicit: next value not explicitly given.
• Linear multistep: the next value is linearly related to the immediate

previous (backward) values (eventually more than one).

• For a fixed step size δ > 0, let tn = t0 + nδ, and xs the approximate
• f the exact x(ts).
• The general Backward Differentiation Formula:

˙ xn = linear combination of xn, xn−1, . . . , x0

• For a Cauchy problem ˙

x = f (x, t), x(t0) = x0, one obtains an implicit equation for xn: xn =

q

• n=1

an x1−n + δ bq f (xn, tn) where the an and bq depends only on the order q.

• E.g., q = 1 (BDF1): xn = xn−1 + δ f (xn, tn)

(a.k.a. Backward Euler)

• K. Ghorbal (INRIA)

5 COMASIC M2 5 / 22

Backward Differentiation Formula (BDF)

Fixed Step Size

• Implicit: next value not explicitly given.
• Linear multistep: the next value is linearly related to the immediate

previous (backward) values (eventually more than one).

• For a fixed step size δ > 0, let tn = t0 + nδ, and xs the approximate
• f the exact x(ts).
• The general Backward Differentiation Formula:

˙ xn = linear combination of xn, xn−1, . . . , x0

• For a Cauchy problem ˙

x = f (x, t), x(t0) = x0, one obtains an implicit equation for xn: xn =

q

• n=1

an x1−n + δ bq f (xn, tn) where the an and bq depends only on the order q.

• E.g., q = 1 (BDF1): xn = xn−1 + δ f (xn, tn)

(a.k.a. Backward Euler)

• K. Ghorbal (INRIA)

5 COMASIC M2 5 / 22

Backward Differentiation Formula (BDF)

Fixed Step Size

• Implicit: next value not explicitly given.
• Linear multistep: the next value is linearly related to the immediate

previous (backward) values (eventually more than one).

• For a fixed step size δ > 0, let tn = t0 + nδ, and xs the approximate
• f the exact x(ts).
• The general Backward Differentiation Formula:

˙ xn = linear combination of xn, xn−1, . . . , x0

• For a Cauchy problem ˙

x = f (x, t), x(t0) = x0, one obtains an implicit equation for xn: xn =

q

• n=1

an x1−n + δ bq f (xn, tn) where the an and bq depends only on the order q.

• E.g., q = 1 (BDF1): xn = xn−1 + δ f (xn, tn)

(a.k.a. Backward Euler)

• K. Ghorbal (INRIA)

5 COMASIC M2 5 / 22

Backward Differentiation Formula (BDF)

Fixed Step Size

• Implicit: next value not explicitly given.
• Linear multistep: the next value is linearly related to the immediate

previous (backward) values (eventually more than one).

• For a fixed step size δ > 0, let tn = t0 + nδ, and xs the approximate
• f the exact x(ts).
• The general Backward Differentiation Formula:

˙ xn = linear combination of xn, xn−1, . . . , x0

• For a Cauchy problem ˙

x = f (x, t), x(t0) = x0, one obtains an implicit equation for xn: xn =

q

• n=1

an x1−n + δ bq f (xn, tn) where the an and bq depends only on the order q.

• E.g., q = 1 (BDF1): xn = xn−1 + δ f (xn, tn)

(a.k.a. Backward Euler)

• K. Ghorbal (INRIA)

5 COMASIC M2 5 / 22

Coefficients of BDFq

Taylor Series (Sundials Implementation)

Ixn = xn (Identity) Nxn = xn+1 (Forward Shift) N −1xn = xn−1 (Backward Shift) Dxn = ˙ xn (Differential) ∆ = I − N (Backward Operator) Observe that N = (I − ∆)−1 (Operator Algebra) Nxn = xn+1 = xn + δDxn + δ2 2! D2xn + δ3 3! D3xn + · · · = (I + δD + δ2 2! D2 + δ3 3! D3xn + · · · )xn = eδDxn Thus: N = eδD

• K. Ghorbal (INRIA)

6 COMASIC M2 6 / 22

δD = ln(N) = ln

• (I − ∆)−1

= ∆ + 1 2∆2 + 1 3∆3 + · · · δ ˙ xn = δDxn = ∆xn + 1 2∆2xn + 1 3∆3xn + · · · BDFq: truncate at order q, for instance for q = 2 δ ˙ xn = ∆xn + 1 2∆2xn = (xn − xn−1) + 1 2(xn − 2xn−1 + xn−2) Thus δ ˙ xn = 3 2xn − 2xn−1 + 1 2xn−2 xn = 4 3xn−1 − 1 3xn−2 + 2 3δf (xn, tn)

• K. Ghorbal (INRIA)

7 COMASIC M2 7 / 22

BDF for Semi-Explicit DAE of index 1

˙ x = f (x, y, t) 0 = g(x, y)

Numerical integration using the BDF2 Scheme

Suppose xn−1 and xn−2 are known, xn and yn are computed by numerically solving the following system (e.g. Newton’s methods) at each iteration: xn = 4 3xn−1 − 1 3xn−2 + 2 3δf (xn, yn, tn) 0 = g(xn, yn)

• BDF converges if m ≤ 6: xi − x(ti) = yi − yi(ti) = o(δm)
• BDF requires a consistent initial condition
• K. Ghorbal (INRIA)

8 COMASIC M2 8 / 22

BDF for Semi-Explicit DAE of index 1

˙ x = f (x, y, t) 0 = g(x, y)

Numerical integration using the BDF2 Scheme

Suppose xn−1 and xn−2 are known, xn and yn are computed by numerically solving the following system (e.g. Newton’s methods) at each iteration: xn = 4 3xn−1 − 1 3xn−2 + 2 3δf (xn, yn, tn) 0 = g(xn, yn)

• BDF converges if m ≤ 6: xi − x(ti) = yi − yi(ti) = o(δm)
• BDF requires a consistent initial condition
• K. Ghorbal (INRIA)

8 COMASIC M2 8 / 22

Newton’s Methods: Generalization

Could be used to numerically approximate the roots of F : Rk → Rk (Fi continuously differentiable functions):

• Multiply by the inverse of the Jacobian

xn+1 = xn − J−1

F (xn)F(xn)

• Or solve the system of linear equations

JF(xn)(xn+1 − xn) = −F(xn) Under some assumptions, the method converges quadratically towards a root of F.

• K. Ghorbal (INRIA)

9 COMASIC M2 9 / 22

Outline

1 BDF Method 2 Collocation (Overview) 3 Brief Introduction to Modelica

• K. Ghorbal (INRIA)

9 COMASIC M2 9 / 22

Collocation

• Collocate means approximate/study an unknown function by means
• f other“simpler”functions
• For instance using polynomial or trigonometric functions

Weierstrass Theorem (1885)

If x(t) is a continuous function on [a, b], then for any given ǫ > 0, there exists a polynomial p(t) such that max

t∈[a,b]|x(t) − p(t)| < ǫ

The theorem has a constructive proof using Bernstein polynomials.

Integration by Collocation

Unlike BDF, Collocation could be used to construct at once n points (t1, x1), . . . , (tn, xn) such that xi approximates the function x(t) at ti for all i. Useful for Consistent Initialization.

• K. Ghorbal (INRIA)

10 COMASIC M2 10 / 22

Collocation

• Collocate means approximate/study an unknown function by means
• f other“simpler”functions
• For instance using polynomial or trigonometric functions

Weierstrass Theorem (1885)

If x(t) is a continuous function on [a, b], then for any given ǫ > 0, there exists a polynomial p(t) such that max

t∈[a,b]|x(t) − p(t)| < ǫ

The theorem has a constructive proof using Bernstein polynomials.

Integration by Collocation

Unlike BDF, Collocation could be used to construct at once n points (t1, x1), . . . , (tn, xn) such that xi approximates the function x(t) at ti for all i. Useful for Consistent Initialization.

• K. Ghorbal (INRIA)

10 COMASIC M2 10 / 22

Collocation

• Collocate means approximate/study an unknown function by means
• f other“simpler”functions
• For instance using polynomial or trigonometric functions

Weierstrass Theorem (1885)

If x(t) is a continuous function on [a, b], then for any given ǫ > 0, there exists a polynomial p(t) such that max

t∈[a,b]|x(t) − p(t)| < ǫ

The theorem has a constructive proof using Bernstein polynomials.

Integration by Collocation

Unlike BDF, Collocation could be used to construct at once n points (t1, x1), . . . , (tn, xn) such that xi approximates the function x(t) at ti for all i. Useful for Consistent Initialization.

• K. Ghorbal (INRIA)

10 COMASIC M2 10 / 22

Collocation: main idea

If p(t) is a polynomial that approximates x(t), solution of our semi-explicit DAE, such that p(ti) = x(ti) = xi for some known ti, where i = 1, . . . , n. Then, for each i, we can compute xi, yi by solving ˙ p(ti) = f (p(ti), yi, ti) 0 = g(p(ti), yi)

• The collocation methods attempts to construct such p as well as

finding appropriate ti.

• Observe that the polynomial p may be dependent on xi, the

unknowns we want to determine.

• Again, Newton’s methods could be used to solve iteratively for p(ti)

and hence xi.

• K. Ghorbal (INRIA)

11 COMASIC M2 11 / 22

Polynomial Interpolation

Suppose we have n points (t1, x1), . . . , (tn, xn), then we construct the polynomial p(t) such that: p(ti) = xi, i = 1, . . . , n

• There exists a (unique) polynomial p of degree n interpolating the n

points p(t) =

n

• i=1

xiLi(t)

• Li are the Lagrange polynomials

Li(t) =

n

• k=0,k=i

t − tk ti − tk , i = 1, . . . , n satisfying in particular Li(tj) = δij (Kronecker delta)

• K. Ghorbal (INRIA)

12 COMASIC M2 12 / 22

Determining the instant (nodes) ti

We want to determine the instants ti ∈ [a, b] such that there exists positive weights wi which make the Gauss quadrature integral exact for the polynomial p, that is I[p] = b

a

p(t)dt =

n

• i=1

wip(ti) = Qn[p] Without loss of generality, we can suppose that a = −1 and b = 1 since b

a

f (t)dt = b − a 2 1

−1

f a + b 2 + t b − a 2

• dt .

Gauss Quadrature Fundamental Theorem (w(x) = 1)

If the ti are the zeros of the Legendre polynomial ℓn, then there exist n weights wi which make the Gauss quadrature integral exact for all polynomials of degree 2n − 1 or less (in particular for our polynomial p of degree n).

• K. Ghorbal (INRIA)

13 COMASIC M2 13 / 22

Legendre Polynomials

Rodrigues’ Formula

ℓn(t) = 1 2nn! dn dtn (t2 − 1)n

• ℓn has degree n
• Legendre polynomials are orthogonal

< ℓi, ℓj >= 1

−1

ℓi(t) ℓj(t)dt = 2 2n + 1δij

• ℓ0(t) = 1, ℓ1(t) = t, ℓ3(t) = 1

2(5t3 − 3t)

• There are effective methods to compute the zeros of ℓn and wi

(Golub-Welsch algorithm)

• Thus we have the ti which depend only on n and not on our
• riginal DAE system !.
• K. Ghorbal (INRIA)

14 COMASIC M2 14 / 22

Example (1/2)

˙ x = x + y 0 = x2 + y2 − 1

• Suppose [a, b] = [0, 0.1], n = 3
• The ti are derived from τi the roots of ℓ3(t): ti = a+b

2

+ τi b−a

2

• τi ∈
• −
• 3

5, 0,

• 3

5

• wi =

b

a Li(t)dt (Lagrange Polynomials)

• wi ∈ { 1

36, 2 45, 1 36}

• p(t) = x1L1(t) + x2L2(t) + x3L3(t)
• ˙

p(t) = x1 ˙ L1(t) + x2 ˙ L2(t) + x3 ˙ L3(t) (˙ Li(tj) = δij)

• K. Ghorbal (INRIA)

15 COMASIC M2 15 / 22

Example (2/2)

System to Solve

˙ p(t1) = x1 + y1 ˙ p(t2) = x2 + y2 ˙ p(t3) = x3 + y3 = x2

1 + y2 1 − 1

= x2

2 + y2 2 − 1

= x2

3 + y2 3 − 1

One solution for [0, 0.1]: x1 = x2 = x3 =

1 √ 2

y1 = y2 = y3 = −1

√ 2

Note that in this particular simple example, the system has a singularity at (x, y) = (1, 0), so not all initial conditions and intervals [a, b] work. In general, such problems occur and one wants to find an initial condition that has no singularities at least as long as t ∈ [a, b].

• K. Ghorbal (INRIA)

16 COMASIC M2 16 / 22

Outline

1 BDF Method 2 Collocation (Overview) 3 Brief Introduction to Modelica

• K. Ghorbal (INRIA)

16 COMASIC M2 16 / 22

What’s Modelica ?

Specification Version 3.3 2014

Modelica is a freely available, object-oriented language for modeling of large, complex, and heterogeneous systems. It is suited for multi-domain modeling, for example, mechatronic models in robotics, automotive and aerospace applications involving mechanical, electrical, hydraulic control and state machine subsystems, process oriented applications and generation and distribution of electric power. Models in Modelica are mathematically described by differential, algebraic and discrete equations. Modelica is designed such that available, specialized algorithms can be utilized to enable efficient handling of large models having more than one hundred thousand equations. Modelica is suited and used for hardware-in-the-loop simulations and for embedded control systems.

• K. Ghorbal (INRIA)

17 COMASIC M2 17 / 22

Modelica

Modelica is a programming language not a tool.

• Domain Specific Language: modeling of physical systems
• Equation based (declarative, acausal, unlike State Flow Languages)
• Object Oriented (object: data + equations)
• Designed (and used) to model and simulate complex physical systems
• Textual and graphical editors

Brief History

• First design group meeting Fall 1996
• Language specification: v1. (1997) – v3.3 revision 1 (2014)
• Modelica Association: open, non profit organization established in

2000

• Modelica conference held annually since 2000
• K. Ghorbal (INRIA)

18 COMASIC M2 18 / 22

Modelica Tools

(Nonexhaustive lists)

Commercial Environments

• Dymola (Dassault Systemes)
• LMS Imagine.Lab Amesim (Siemens)
• SystemModeler (Wolfram)
• MapleSim (Maplesoft)

Open/Free Environments

• OpenModelica (OSMC)
• JModelica.org
• Modelicac (Scilab)
• K. Ghorbal (INRIA)

19 COMASIC M2 19 / 22

OpenModelica Environment

• K. Ghorbal (INRIA)

20 COMASIC M2 20 / 22

Crash Course (DrModelica)

• Hello World Example
• DAE Example
• Inheritance
• Connectors, Connections and Equations Generation
• K. Ghorbal (INRIA)

21 COMASIC M2 21 / 22

References

• Uri M. Ascher and Linda R. Petzold, Computer Methods for Ordinary

Differential Equations and Differential-Algebraic Equations, 1998.

• K. E. Brenan, S. L. Campbell and L. R. Petzold, Numerical Solution
• f Initial-Value Problems in Differential-Algebraic Equations, 1996.
• G. H. Golub and J. H. Welsch, Calculation of Gauss Quadrature

Rules, 1969.

• computation.llnl.gov/projects/sundials/ida
• www.openmodelica.org
• K. Ghorbal (INRIA)

22 COMASIC M2 22 / 22