Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations Dr. Abebe Geletu Ilmenau University of Technology Department of Simulation and Optimal Processes (SOP)

Winter Semester 2011/12

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Topics Newton’s Law: m¨ x = − Fl x m¨ y = mg Fl y Conservation of mechanical energy: x2 + y 2 = l2

(DAE)

x˙ 1 = x3 x˙ 2 = x4 F x1 ml F x˙ 4 = g x2 l 0 = x2 + y 2 − l 2 . x˙ 3 = −

1 2

Numerical Methods of Ordinary Differential Equations Initial Value Problems (IVPs)

Lecture 3 Introduction1 to Numerical Methods for Differential and Differential Algebraic Equations TU Ilmenau

Single Step Methods Multi-step Methods

Initial Value Problems General form of an initial value problem (IVP) is x(t) ˙ = f (t, x(t)), t0 ≤ t ≤ tf ,

(1)

x(t0 ) = x0

(2)

(given),

where x and f can be vector functions:     f1 (t, x(t)) x1 (t)  x2 (t)   f2 (t, x(t))      x(t) =  .  , f (t, x(t)) =  . .. .  .    . xn (t) fn (t, x(t)) If f is (explicitly) time independent; i.e. f (t, x(t)) = f (x(t)), then the system (1) - (2) is said to be autonomous. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

x(t0 ) = x0 is the initial condition (start time of the motion) The expression x(t) ˙ represents the first order derivative of x(t) d with respect to time; i.e. x(t) ˙ = dt x(t). Thus, the IVP (1) - (2) represents a system of first order differential equations. If f (t, x(t)) = A(t)x(t), then we have x(t) ˙ = A(t)x(t), x(t0 ) = x0 is a first order linear IVP. In the autonomous case we have A(t) = A. Example: Equation of motion of a simple pendulum: θ˙ = ω

g ω˙ = − sin(θ) L ⇒   θ˙ ω˙

 θ(0) = π4 ω , g − L sin(θ) ω(0) = 0

 =

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Commonly, numerical methods are developed for systems of first order differential equations. ⇒ any higher order differential equation should be written as a system of first order differential equations. NB: The Matlab ODE Toolbox works only with systems of first order differential equations. Strategy: Convert any n-th order initial value problem into a system of first-order IVP by introduction new variables.   dx d2 x dn−1 x dn x(t) , = f t, , . . . , dtn dt dt2 dtn−1 dx(t) dn−1 x(t) x(t0 ) = a1 , x(t ˙ 0 ) = a2 , = a3 , . . . , = an . dt t=t0 dtn−1 t=t0 Introduce new variables: 2 x1 = x, x2 = x, ˙ x3 = d dtx(t) 2 , . . . , xn =

dn−1 x(t) . dtn−1

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Then the n-th order IVP reduces to x˙ 1 = x2 ,

x1 (t0 ) = a1

x˙ 2 = x3 , .. .

x2 (t0 ) = a2

x˙ n = f (t, x1 , x2 , . . . , xn ), xn (t0 ) = an . Example: The equation of motion of a simple pendulum (neglecting air resistance) with mass m is a 2nd-order equation mLθ¨ + mg sin(θ(t)) = 0, (3) ˙ θ(t0 ) = θ0 , θ = ω0 . (4) ˙ Then Introduce ω = θ. θ˙ = ω,

θ(t0 ) = θ0 g ω˙ = − sin(θ), ω(t0 ) = ω0 . L

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Objective: To determine a continuous function x∗ (t) that satisfies the dynamic equation (1); i.e. x∗ (t) = f (t, x∗ (t)) and the initial condition (2); i.e. x∗ (t0 ) = x0 . !!! In general, the solution function x∗ (t) of the IVP is not easy to be analytically determined.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Basic idea of numerical solution of differential equations: Select a set of equally spaced time points t1 , t2 , . . . , tN from the interval [t0 , tf ] (tf freely chosen for IVPs). That is, t1 = t0 + h , t2 = t1 + h, and so on; so that h = tk+1 − tk , k = 0, 1, . . . , N − 1. Approximate the continuous solution x∗ (t) of (1) & (2) by values at the discrete time points t1 , t2 , . . . , tN ; i.e. x(t1 ) = x1 , x(t2 ) = x2 , . . . , x(tN ) = xN . determine the values x1 , x2 , . . . , xN through the IVP (1) - (2). Note: more discretization points from [t0 , tf ] implies better accuracy of the approximation; but demands more computation effort, depending on the algorithm used, it may be necessary to solve a system of linear or nonlinear equations. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Euler’s Method Motivation: for a single IVP x(t) ˙ = f (t, x(t)), x(t0 ) = x0 ; i.e. x(t) ∈ R, 0 ≤ t ≤ tf . Given the point (t0 , x0 ) the value x(t ˙ 0 ) = dx(t) is the slope of dt t=t0

the tangent line to the graph of x(t) at the point (t0 , x0 ).

Figure: Euler’s Method Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Since the function x(t) is not yet known, the derivative (slope) can be approximated by x(t0 + h) − x(t0 ) x(t ˙ 0) ≈ h Thus, using x(t ˙ 0 ) = f (t0 , x0 ), we have x(t0 + h) − x(t0 ) = f (t0 , x(t0 )) h This implies x(t0 + h) = x(t0 ) + hf (t0 , x(t0 )) ⇒ x1 = x(t1 ) = x(t0 + h) = x(t0 ) + hf (t0 , x(t0 )) From this we obtain the value of x1 = x(t1 ). By the same procedure x2 can be also determined. In general, once xk is known, xk+1 can be found from xk+1 = xk + hf (tk , xk ), k = 0, 1, . . . , m − 1. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Euler’s Algorithm Select N . (t −t ) Set h = f N 0 ; x0 = x(t0 ) FOR k = 0 : (n − 1) xk+1 = xk + hf (tk , xk ); tk+1 = tk + h; ENDFOR. Matlab: function [t,x]=myEuler(@f,x0,t0,tf,N) h=(tf-tn)/N; t=[t0:h:tf]; %discrete time points between, and including, x(1)=x0; for k=2:N x(k)=x(k-1) + h*feval(f,t(k-1),x(k-1)); end; Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Simple Pendulum with Euler’s Method function time=MyPenEuler(N,tf) t=cputime; m=1; L=9.81; g=9.81; beta=0; t0=0; theta0=(1/10)*pi; omega0=0; h=(tf-t0)/N; tt=[t0:h:tf]; ttheta=zeros(size(tt)); oomega=zeros(size(tt)); ttheta(1)=theta0; oomega(1)=omega0; Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Simple Pendulum with Euler’s Method...

for k=1:N oomega(k+1)=oomega(k)+h*... (-beta/m/L*oomega(k)-g/L*sin(ttheta(k))); ttheta(k+1)=ttheta(k)+h*oomega(k); end time=cputime-t; plot(tt,oomega,tt,ttheta) title(’Simple pendulum using Euler’’s method’)

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Advantages/Disadvantages of Euler’s Method Advantages: Euler’s Method is simple and direct Can be used for nonlinear IVPs Disadvantages: it is less accurate and numerically unstable. Approximation error is proportional to the step size h. Hence, good approximation is obtained with a very small h. This requires a large number of time discretization leading to a larger computation time. Usually applicable to explicit differential equations.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Improvements on Euler’s Method (A) The Mid-Point Method Let k1 = hf (ti , xi ). Note that: k1 = hf (ti , xi ) = xi+1 − xi . But to compute k1 , it is not necessary to know the value of xi+1 .

Figure: Use also the mid point Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Improvements on Euler’s Method... • The point (ti + 21 ti , xi + 12 k1 ) is the mid-point of the segment that connects the points (ti , xi ) and (ti+1 , xi+1 ). • Now, instead of the single slope at (ti , xi ), we use the slopes at (ti , xi ) and at (ti + 12 ti , xi + 12 k1 ) to approximate xi+1 . ⇒ This leads to the Mid-Point Runge-Kutta Algorithm: The mid-point Runge-Kutta Method k1 = hf (ti , xi ) 1 1 k2 = hf (ti + h, xi + k1 ) 2 2 xi+1 = xi + k2 .

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Improvements on Euler’s Method... (B) Methods that use Intermediate Points Question: Can we improve the mid-point RK method by freely choosing any point from the segment connecting (ti , xi ) and (ti+1 , xi+1 )?

Figure: Choose freely any intermediate point Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

• A weighted sum of f (ti , xi ) and f (ti + ch, xi + mk1 ) provideds more slope information than f (ti , xi ) alone. Scheme: (RK2)

x(t + h) = x(t) + hfe(t, x)

where fe(t, x) = [w1 f (t, x) + w2 f (t + ch, x + mk1 )] . How to determine w1 , w2 , c and h? We use: (i) Taylor approximation of f (t + ch, x + mk1 ): f (t + ch, x + mk1 ) = f (t, x) + chft (t, x) + mk1 fx (t, x)x0

(∗)

Put this into (*) to obtain x(t+h) = x(t)+h [w1 f (t, x) + w2 {f (t, x) + chft (t, x) + mk1 fx (t, x)})] ⇒ (with k1 = hf (t, x))

x(t+h) = x(t)+h (w1 + w2 ) f (t, x)+w2 ch2 ft (t, x)+mw2 h2 fx (t, x)f (t, x) Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Improvements on Euler’s Method... (ii) 2nd- order Taylor approximation of x(t + h): x(t + h) = x(t) + hx0 (t) +

h2 00 x (t) 2

(∗∗)

h2 00 x (t) (∗∗) 2 Recall x0 (t) = x(t) ˙ = f (t, x) ⇒ x00 (t) = ft (t, x) + fx (t, x)f (t, x). Putting these into (**) we obtain x(t + h) = x(t) + hx0 (t) +

x(t + h) = x(t) + hf (x, t) +

h2 h2 ft (t, x) + fx (t, x)f (t, x)) 2 2

Comparing (*) and (**) we find

(∗∗)

w1 + w2 = 1, cw2 = 12 , mw2 = 12 .

A system of three equations with four unknowns ⇒ This system has infinitely many solutions. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Improvements on Euler’s Method... • By fixing one of the unknowns w1 , w2 , c and m, we can solve the system of equations for the other three. ⇒ Hence, we can construct several RK2 methods. Example: • If we fix w1 = 0, then w2 = 1, c = 12 and m = 21 . Hence, (RK2) reduces to the mid-point RK2 method. w1 0

w2 1

c

m

1 2 1 3

1 2 2 3

1 2

1 2

1

1

3 4

3 4

Method Modified Euler Method Heun’s Method Ralston’s Method

Heun’s Method: k1 = f (ti , xi )   xi+1 = xi + h 12 f (ti , xi ) + 12 f (ti + h, xi + k1 ) Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Runge-Kutta Methods The above schemes are known as second-order Runge-Kutta methods, since x(t + h) has been approximated by a second degree Taylor polynomial. I For an accurate approximation use more intermediate points between (ti , xi ) and (ti + h, xi + k1 ) and include more slope information. I Requires higher degree Taylor approximations of x(t + h).

Figure: Use several intermediate points Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Runge-Kutta Methods... General Scheme: [RK(m)]

xi+1 = xi + hfe(ti , xi );

where fe(ti , xi ) = w1 f (ti , xi ) + w2 f (ti + α1 h, xi + β1 k1 ) + . . . + wm f (ti + αm h, xi + βm km ).

• The values k1 , k2 , . . . , km are found at the intermediate points. • w1 + w2 + . . . + wm = 1. • The relation between ω1 , . . . , ωm , α1 , . . . , αm , β1 , . . . , βm is determined by comparing the m-th order Talyor approximation of x(t + h) and the scheme RK(m)

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

I Hence, the resulting method is called an m-th order Runge-Kutta method. I There are third-order Runge-Kutta (RK3), fourth-order Runge-Kutta methods (RK4), etc. The classical 4-th order Runge-Kutta Method has the form: (RK4s) k1 = f (ti , xi ) 1 1 k2 = hf (ti + h + xi + k1 ) 2 2 1 1 k3 = hf (ti + h + xi + k2 ) 2 2 k4 = hf (ti + h, xi + k3 ) 1 1 1 1 xi+1 = xi + k1 + k2 + k3 + k4 6 3 3 6

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

The classical 4th order Runge-Kutta Method function myRK4(N,tf) m=5; %Mass L=9.81; %length of the string g=9.81; %gravitaional acceleration beta=0; % t0=0; theta0=(1/10)*pi; omega0=0; h=(tf-t0)/N tt=[t0:h:tf]; ttheta=zeros(size(tt)); oomega=zeros(size(tt)); ttheta(1)=theta0; oomega(1)=omega0; Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

%The classic 4th order Runge Kutta method for k=1:N k1=h*(-beta/m/L*oomega(k)-g/L*sin(ttheta(k))); m1=h*oomega(k); k2=h*(-beta/m/L*(oomega(k)+k1/2)-g/L*(sin(ttheta(k)+m1/2 m2=h*(oomega(k)+m1/2); k3=h*(-beta/m/L*(oomega(k)+k2/2)-g/L*(sin(ttheta(k)+m2/2 m3=h*(oomega(k)+k2/2); k4=h*(-beta/m/L*(oomega(k)+k3/2)-g/L*(sin(ttheta(k)+m3/2 m4=h*(oomega(k)+k3); oomega(k+1)=oomega(k)+1/6*(k1+2*k2+2*k3+k4); ttheta(k+1)=ttheta(k)+1/6*(m1+2*m2+2*m3+m4); end % Plot the solutions. plot(tt,oomega,tt,ttheta) title(’Simple pendulum solved with 4th order Runge Kutta’) Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

In general the 4-th order Runge-Kutta Methods take the form: (RK4) k1 = f (ti , xi ) k2 = hf (ti + c2 h, xi + a21 k1 ) k3 = hf (ti + c3 h, xi + a31 k1 + a32 k2 ) k4 = hf (ti + c4 h, xi + a41 k1 + a42 k2 + a43 k3 ) xi+1 = xi + ω1 k1 + ω2 k2 + ω3 k3 + ω4 k4 where the unknown parameters are organized as 0 c2 c3 c4

a21 a31 a41 ω1

a32 a42 ω2

a43 ω3

ω4

known as Butcher’s table. As a result it is possible to construct various 4th order Runge-Kutta methods. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Question: When do RK methods fail? Example: x˙ 1 = −100x1

(5)

x˙2 = −0.1x2 .

(6)

The state x1 (t) = x1 (0)e−100t and x2 (t) = x2 (0) = e−0.1t . x2 (t) decays faster than x1 (t). So it is not reasonable to compute such systems with the same step length h.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Advantages and Disadvantages of Simple RK Methods Advantages they are easy to implement the are stable Disadvantages they require relatively large computer time error estimation are not easy to be done the above simple RK methods do not work well for stiff differential equations (eg. linear differential equations with widely spread eigenvalues) In particular, they are not good for systems of differential equations with a mix of fast and slow state dynamics. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Adaptive Runge-Kutta Methods The choice of a good step size h is crucial. In particular, this is important for stiff differential equations: • too large h may lead to inaccurate approximation of the solution; • too small h forces the algorithm use large computation resources. ⇒ The selection of the step length h is a compromise between numerical accuracy and speed of computation. I Adaptive adjustment of the step size h reduces (local) truncation error for approximating x(ti+1 ) by xi+1 : kx(ti+1 ) − xi+1 k Recall the Taylor approximation: h3 ... h4 .... h2 x (ti ) + x (ti ) + . . . x ¨(ti ) + 2 3! 4! = xi + C1 h + C2 h2 + C3 h3 + C4 h4 + . . .

x(ti + h) = x(ti ) + hx(t ˙ i) +

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Adaptive Runge-Kutta Methods... • Adaptive strategies adjust the step-length from iteration to iteration by testing the accuracy iterates. • Usually this test of accuracy is done by comparing lower and higher order RK algorithms. ⇒ Hence, adaptive RK methods combine lower and higher order RK methods to adjust the step-size h.

(RK3) k1 = f (ti , xi ) 1 1 k2 = hf (ti + h + xi + k1 ) 2 2 k3 = hf (ti + h, xi − k1 + 2k2 ) h RK3 xi+1 = xi + (k1 + 4k2 + k3 ) . 6 Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

(RK4s) k1 = f (ti , xi ) 1 1 k2 = hf (ti + h + xi + k1 ) 2 2 1 1 k3 = hf (ti + h + xi + k2 ) 2 2 k4 = hf (ti + h, xi + k3 ) h RK4 xi+1 = xi + (k1 + 2k2 + 2k3 + k4 ) 6 RK4 RK4 • Both xi+1 and xi+1 are approximations for x(ti+1 ). Thus the local truncation errors (using the Taylor approximation of x(t + h) with ti+1 = ti + h): 3 from (RK3) is x(ti+1 ) − xRK3 i+1 (ti+1 ) = C3 h + . . . RK4 4 from (RK4) x(ti+1 ) − xi+1 = C4 h + . . . ⇒ (subtracting the second from the first) RK3 3 xRK4 i+1 − xi+1 ≈ C3 h Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

I Hence, the truncation error can be estimated in two useful ways: (i) on the hand, the local truncation error is can be approximated by RK3 3 ErrT runc = |xRK4 i+1 − xi+1 | ≈ kC1 kh

(7)

In the next iteration step we would like to choose a new step-size hnew so that the truncation error below a give tolerance: ErrT runcnew ≈ kC1 kh3new < ErrT ol.

(8)

I The new discretization step-size hnew should be chosen in such a way that the resulting truncation ErrTruncnew lies below a given tolerance level ErrTol. I The step-size should be adjusted relative to the truncation error committed in the previous step and the given error tolerance for the next step (taking ratios of terms in equations (7) and (8)) kC1 kh3new ErrT ol < . kC1 kh3 ErrT runc Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Adaptive Runge-Kutta Methods... ⇒  hnew < h

ErrT ol ErrT runc

1/3

In particular the new h can be computed adaptively by  hnew = γh

ErrT ol ErrT runc

1/3

for a fixed γ such that 0 < γ < 1. (ii) on the other hand, comparing the algorithms (RK3) and (EK4s)s, RK3 ErrT runc = |xRK4 i+1 − xi+1 | =

h [−2k2 + k3 + k4 ] 6

These lead to an adaptive 3-4 Runge-Kutta algorithm. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Adaptive Runge-Kutta Methods...an RK34 • Set t0,x0, initial step size h, ErrorT ol, tf and γ, while t < tf do • Compute k1 = f (t, x) k2 = hf (ti + 12 h + xi + 12 k1 ) k3 = hf (ti + h, xi + 21 k2 ) k4 = hf (ti + h, xi + k3 ) • Set ErrT runc ← h6 k2k2 − k3 − k4 k. if ErrT runc ≤ ErrorT ol then • Set xi+1 ← xi + h6 (k1 + 2k2 + 2k3 + k4 ) • Set t ← t + h end if
ErrorT ol 1/3 • New step-size: h ← γh ErrT runc end while Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Common adaptive methods are the Runge-Kutta-Fehlberg(RKF) formulas. • RKF23 (a combination of 2nd and 3rd order RKs), RKF45, etc. Matlab ode functions: I ode23 is an adaptive Runge-Kutta methods and a combination of a 2nd and 3rd order Runge-Kutta methods. I ode45 is an adaptive RK known as the Dormand-Prince method. Example: Mass-springer dumper system: Mx ¨ + C x˙ + kx = F with mass M = 5kg, dumper value C = 10N s/m, spring constant k = 1000N/m, F = sin(20t)N with forcing frequency ω = 20. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Mass-springer dumper system...Matlab To solve 5¨ x + 10x˙ + 100x = sin(20t), with x(0) = 0, x(0) ˙ = 0. ⇒ x˙ 1 = x2 , x˙ 2 = −2x2 − 20x1 + 0.2 ∗ sin(20t). function springDumper(N,tf) t0=0; h=(tf-t0)/N; x10=0; x20=0; tt=t0:h:tf; [t,x]= ode23(@diffEq,tt,[x10,x20]); subplot(2,1,1), plot(t,x(:,1)); title(’Result from Ode23’) [t,x]= ode45(@diffEq,tt,[x10,x20]); subplot(2,1,2), plot(t,x(:,1),’r’); title(’Result from Ode45’) end Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Mass-springer dumper system...Matlab...

function dxdt=diffEq(t,x) dx1dt = x(2); dx2dt = -2*x(2)-20*x(1)+0.2*sin(20*t); dxdt=[dx1dt;dx2dt]; end

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Adaptive Runge-Kutta Methods...Advantages/Disadvantages Advantages: • provide better precision in the numerical approximation state trajectories; • step-sizes can be automatically adapted to the dynamics of the states; ⇒ process dynamics with stiff-differential equation models can be efficiently solved Disadvantages: • such algorithms may take enormous computation time; • simulation of dynamic systems with several state equations can be computationally demanding; Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multi-step and Implicit Methods Euler, Runge-Kutta and adaptive RK methods are single-step and explicit methods. • Multi-step methods: use a set of earlier iterates up to the current iterate xi to compute the next iterate xi+1 . ⇒ Earlier iterates xi , xi−1 , . . . , xi−m are used to compute xi+1 . For example: a two-step method uses: xi and xi−1 to compute xi+1 a three-step method uses: xi , xi−1 and xi−2 to compute xi+1 • Implicit methods: x˙ i+1 = f (ti+1 , xi+1 ), where x˙ i+1 = x(t)| ˙ t=ti+1 . Hence, xi+1 appears implcitily in f (ti+1 , xi+1 ). Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multi-step and Implicit Methods... I BDF is a well known efficient multi-step method: BDF = Backward Differentiation Formula I While IRK is a well known implicit method: IRK = Implicit Runge-Kutta Method I Both BDF and IRK are widely used for the direct discretization of continuous dynamic optimization problems. I Particularly, when the IVP involves a stiff-differential equations, BDF and IRK are methods of choice. Stiff differential equations: A differential equation x˙ = f (t, x) is stiff if the (Jacobi) matrix ∂f ∂x has eigenvalues with widely spread negative real parts. ⇒ some states change fast while others are slower. Examples: mass-spring systems with large spring constants, (chemical) process engineering models, etc. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multi-step and Implicit Methods...BDF I Adaptive RK methods can solve some stiff problems, but they are in general not efficient for a system with several states. (A) BDF: Given xi , xi−1 , . . . , xi+1−m to compute xi+1 . Find an m-th degree polynomial P that interpolates the m + 1 points (ti+1 , xi+1 ), (ti , xi ), (ti−1 , xi−1 ), . . . , (ti+1−m , xi+1−m ), A good choice for the interpolating polynomials P is m X P (t) = xi+1−j Lj (t) j=0

with the Lagrange polynomial  m  Y t − ti+1−l , j = 0, 1, . . . , m. Lj (t) = ti+1−j − ti+1−l l=0 l 6= j Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

BDF... Observe that P (ti+1−j ) = xi+1−j , j = 0, 1, . . . , m. In particular P (ti+1 ) = xi+1 . Thus, replace x˙ i+1 by P˙ (ti+1 ) to obtain P˙ (ti+1 ) = f (ti+1 , xi+1 ).

(∗)

But P˙ (ti+1 ) =

m X

xi+1−j L˙ j (ti+1 ) = xi+1 L˙ 0 (ti+1 ) +

j=0

m X

xi+1−j L˙ j (t)

j=1

Putting this into (*) we get: xi+1 = −

m X j=1

xi+1−j

L˙ j (ti+1 ) 1 + f (ti+1 , xi+1 ) L˙ 0 (ti+1 ) L˙ 0 (ti+1 )

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Define aj =

L˙ j (ti+1 ) 1 , bm = . ˙ ˙ L0 (ti+1 ) hL0 (ti+1 )

An m-step BDF Algorithm (BDF m): Given a1 , . . . , am ,bm xi+1 = −

m X

aj xi+1−j + bm hf (ti+1 , xi+1 )

j=1

Example:(a) A BDF 2 algorithm   4 1 2 xi+1 = xi − xi−1 + hf (ti+1 , xi+1 ) 3 3 3 (b) A BDF 3 algorithm   4 9 2 6 xi+1 = xi − xi−1 + xi−2 + f (ti+1 , xi+1 ) 3 11 11 11 Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

BDF... • If f is a nonlinear function, all BDF algorithms require the solution of nonlinear equaitons. =⇒ BDF is used in conjuction with a Newton method. I Matlab function: ode15s uses a BDF algorithm. Advantages: • are suitable for the solution of stiff differential equations • can be easily generalized to solve DAEs • BDF 2 up to BDF 6 are known to be stable • widely used in circuit simulation Disadvantages: • may require the solution of nonlinear equations • faces difficulties if the step-size h is constant (i.e. requires adaptation of the step size). Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multi-step and Implicit Methods...IRK (B) IRK: I If xi is an approximation of x(ti ), how to determine the approximation xi+1 of x(ti+1 )? Idea: Approximate x(ti+1 ) by xi+1 = P (ti+1 ), where P is a polynomial satisfying (the requirements): P (ti ) = xi P˙ (ti + τl h) = f (ti + τl h, P (ti + τl h)), l = 1, . . . , m, where τ1 , τ2 , . . . , τm ∈ (0, 1) and h = ti+1 − ti the steps size. • Note that since P is a polynomial, P˙ is also a polynomial.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

(9) (10)

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

IRK.. Now define P˙ as P˙ (ti + τ h) =

m X

f (ti + τl h, P (ti + τl h))Ll (τ ),

l=1

with  m  Y τ − τj Ll (τ ) = τl − τj j=1 l 6= j

(11)

• This special polynomial satisfies equation (10). • Integrating equation (9) on both sides Z τj Z τj m X P˙ (ti + τ h) = f (ti + τl h, P (ti + τl h)) Ll (τ ), j = 1, . . . , m. 0

l=1

0

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

IRK.. It follows that (i) P (ti +τj h) = xi +h

m X

Z

τj

f (ti +τl h, P (ti +τl h))

Ll (τ ), j = 1, . . . , m. 0

l=1

and for τj=1 (ii) P (ti + h) = xi + h

m X

Z f (ti + τl h, P (ti + τl h))

1

Ll (τ )dτ. 0

l=1

Now, define Z τj ajl = Ll (τ ) 0 Z 1 bl = Ll (τ ) 0 Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

(12) (13)

IRK.. An m-stage IRK Algorithm: Given the points τ1 , τ2 , . . . , τm ∈ [0, 1): kl = xi + h

m X

alj f (ti + τj h, kj ), l = 1, . . . , m

(14)

j=1

xi+1 = xi + h

m X

bl f (ti + τl h, kl )

(15)

l=1

• If τ1 , . . . , τm are known the vector b = (b1 , . . . , bm ) and the matrix A = (alj ) are easy to determine from equations (14) and (15). Butcher Tabeleau: τ

A b>

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

The difference in IRK methods lies: • on the selection of the parameters τ1 , . . . , τm • on special adjustments to reduce the size of the matrix A. • IRK methods are also known as collocation methods. Some examples: √ (a) Gauss-Legendre τ1 = 1/2 − 15/10, τ2 = 1/2 and √ τ3 = 1/2 + 15/10 (are roots of the 3rd degree Legendre Polynomila): √ √ √ 15 15 5 2 5 1/2 − 15/10 − − 36 9 15 36 √ √30 15 5 15 2 5 1/2 36 24 9√ 36 − 24 √ √ 5 5 1/2 + 15/10 36 + 3015 92 + 1515 36 5 18

4 9

(b) Radau IIA (c) Radau 5 (d) Lobatto IIIA (e) Diagonally implicit RK methods Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

5 18

IRK.. I Matlab function: ode23tb is an implementation of an IRK algorithm. Example: A continually stirred homoiothermal reactor

A continuously stirred reactor (CSTR) with feed A and product B. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

• • • • • •

Concentration of A in the reactor: CA The feed concentration: CAf Mixture temperature in the CSTR: T Bulk volume of the mixture in the CSTR V Density and specific heat constants: ρ, Cp k0 reaction constant

Constants: k0 = 3.493x107h−1 , E = 11843kcal/kmol, (−∆H) = 5960kcal/kmol, ρCp = 500kcal/m3/K U A = 150kcal/h/K, R = 1.987kcal/kmol/K, V = 1m3 , q = 1m3 /h, CAi = 10kmol/m3 , Ti = 298K, Tj = 298K

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

dCA dt dT dt

= =

q (CAi − CA ) k0 exp(−E/RT )CA V q (−∆H)k0 exp(−E/RT )CA UA (Ti − T ) + + (Tj − T ) V ρCp ρCp V

Objective: Study the dynamics of the CSTR. xo2 =[6,399]; df=@(t,x)cstr(x); [t2,x2]=ode15s(df,[0,100],xo2,[]); subplot(3,1,1); [t2,x2]=ode23tb(df,[0,100],xo2,[]); subplot(3,1,2); plotyy(t2,x2(:,1),t2,x2(:,2)); subplot(3,1,3); [t2,x2]=ode45(df,[0,100],xo2,[]); plotyy(t2,x2(:,1),t2,x2(:,2)); Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

IRK...Advantages/Disadvantages Advantages: • can be used to solve stiff differential equations and DAEs • widely used for the discretization of optimal control problems Disadvantages: • to guarantee stability of higher-order IRK some of the τl ’s need to be chosen outside [0, 1] (J. C. Butcher) • may require the solution of nonlinear equations implying intensive computations • adaptive selection of the step-size h can be resource intensive • multi-step IRK is still a research topic ⇒ The method of ”orthogonal collocation on finite elements” is a generalizations of IRK methods. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Boundary Value Problems Note that: The solution x∗ (t) of an IVP depends on the initial values t0 and x0 . If we change one or both of these values, the solution x∗ (t) also changes. This dependence is indicated by writing x∗ (t; t0 , x0 ). In general, a boundary value problems has the form (BV P )

x˙ = f (t, x), x(t0 ) = A, x(tf ) = B.

• The function f (t, x) is assumed continuous and differentiable. • The solution x(t) is required to satisfy the two condones x(t0 ) = A and x(tf ) = B known as boundary conditions. Objective: to compute a function x(t) that satisfies the BVP.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Boundary Value Problems (BVPs) Difficulties with BVPS: A BVP (in contrast to IVPs) can have several solutions (i.e., there can be several trajectories (paths) to arrive at x(tf ) = B starting from x(t0 ) = A.)

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Single-Shooting Strategy: BVPs are solved by solving a sequence of IVPs Hence, any one of: Euler, Runge-Kuta, BDF, IRK, etc. can be used in the solution strategy of BVPs.

Single-Shooting Idea of single-shooting: from the initial position x(t0 ) shoot with an angle to hit the target at x(tf ).

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Single-Shooting... • Question: What is the shooting angle? Note that: the shooting (aiming) angle is related with the slope x(t ˙ 0 ) of the trajectory x(t) (since tan α = x(t ˙ 0 )). • Nevertheless, x(t ˙ 0 ) is not given (unknown). That is, x(t ˙ 0 ) =? • So we make a guess and take x(t ˙ 0 ) = u. (But, with this guess (trial shooting angle) we may not hit the target). • Now, drop the end-time condition x(tf ) = B and consider the initial value problem x˙ = f (t, x), x(t0 ) = A, x(t ˙ 0 ) = u.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

• Use the transformation x1 = x, x2 = x˙ 1 ⇒ x˙ 2 = x ¨1 = x ¨ to write the problem as a system of first-order IVPs x˙ 1 = x2 , x1 (t0 ) = A

(16)

x˙ 2 = ft (t, x1 ) + fx (t, x1 )x2 , x2 (t0 ) = u.

(17)

• Solve this IVP  using any  one of the methods for IVPs to obtain a ∗ (t) x 1 solution x∗ (t) = . x∗2 (t) • If x∗1 (t) satisfies the end-point boundary condition at t = tf ; i.e. if x∗1 (tf ) = B then we are done: x∗1 (t) solves the (BVP). Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Otherwise, x∗1 (tf ) 6= B. • Our earlier guess is wrong and we have to adjust our shooting (aiming) angle, i.e. take another guess for u. • Question: How to make a good guess for u? • We known that the solution x∗ (t) depends on the initial values of the IVP, so we write x∗ (t; u). • The solution of the system of equations F (u) = x∗1 (tf , u) − B = 0 provides a good guess for u. In general, this can be a system of nonlinear equations. • We use the Newton’s method to solve such system of nonlinear equations. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

A Newton Algorithm Initialize: Set u0 and tolerance ε while kF (u(k) )k > ε do • Compute the vector F (u(k) ) and the Jacobian matrix F 0 (u(k) ). • Solve the system of linear equations F 0 (u(k) )d = −F (u(k) ) to obtain dk as a solution. • Set u(k+1) = u(k) + dk end while • To solve the system of nonlinear equations we need the values of F (u(k) ) and F 0 (u(k) ). • The value of F (u(k) ) is found as F (uk ) = x∗ (tf , uk ) − B. • It remains to compute the derivative F 0 (u(k) ). Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

We us the forward-difference numerical approximation for F 0 (uk ): F 0 (u(k) ) =

F (u(k) + ∆) − F (u(k) ) ∆

for a given ∆. Note that F (u(k) + ∆) = x∗1 (tf , u(k) + ∆) − B. • Since x∗1 (tf , u(k) + ∆) is unknown, we need to re-solve the problem (16) & (17) using x2 (t0 ) = u(k) + ∆ x˙ 1 = x2 , x1 (t0 ) = A

(18) (k)

x˙ 2 = ft (t, x1 ) + fx (t, x1 )x2 , x2 (t0 ) = u

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

+ ∆.

(19)

Single Shooting ... A Single-Shooting Algorithm Select u0 , ∆ > 0 and tolerance ε > 0 while kF (u(k) )k > ε do • Compute inital value problems (16) & (17) with x2 (tf ) = u(k) and (18) & (19) with x2 (tf ) = u(k) + ∆ • Compute the vector F (u(k) ) and the Jacobian matrix F 0 (u(k) ): F (u(k) + ∆) − F (u(k) ) ∆ • Solve the system of linear equations F 0 (u(k) ) =

F 0 (u(k) )d = −F (u(k) ) to obtain dk as a solution. • Set k ← k + 1 • Set u(k) ← u(k) + dk end while Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Single Shooting ... Advantages/Disadvantages

Advantages: • Simple to implement • effective when the interval [t0 , tf ] is short. Disadvantages: • for some problems the algorithm can become un- stable • requires good initial guesses for the IVPs • the solution strategy is highly influenced by the length of the shooting (distance) interval [t0 , tf ]. Hence, large numerber of iterations may be needed.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multiple-Shooting Method • Consider again the problem x˙ = f (t, x), x(t0 ) = A, x(tf ) = B.

(20)

• In the multiple-shooting method, the interval [t0 , tf ] is divided into subintervals [t0 , t1 ], [t1 , t2 ], . . . , [tm−1 , tm ] by using time-points t0 = t1 , t2 , . . . , tm−1 , tm = tf . • On each subinterval tj ≤ t ≤ tj+1 , j = 0, 1, . . . , m − 1, solve the initial value problems x˙ = f (t, x),

x(tj ) = uj , tj ≤ t ≤ tj+1 , j = 0, 1, . . . , m − 1, (21)

to get a solution x(t; tj , uj ) .

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multiple-Shooting Method...

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multiple-Shooting Method... • Test if all of the following equations hold true at the end points of each sub-interval x(t1 ; t0 , u0 ) = u1

(22)

x(t2 ; t1 , u1 ) = u2 .. .

(23)

x(tm−1 ; tm−2 , um−2 ) = um−1 x(tm ; tm−1 , um−1 ) = um = B. • If yes, we are done. The solution obtained satisfies the (BVP). • Otherwise, make a new and a better guess of the u’s.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

(24) (25) (26)

• To find a better guess for the u0 s, solve the system of (possibly nonlinear) equations: F1 (u) = x(t1 ; t0 , A) − u1 = 0

(27)

F2 (u) = x(t2 ; t1 , u1 ) − u2 = 0 .. .

(28)

Fm−1 (u) = x(tm−1 ; tm−2 , um−2 ) − um−1 = 0 Fm (u) = x(tm ; tm−1 , um−1 ) − B = 0. • Put all the guess values into a variable vector u as   u1  u2    u= .  .  .  um−1 • Note that, each of the uj ’s can be themselves vectors. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

(29) (30) (31)

• Hence, solve the system of equations   F1 (u)  F2 (u)    F (u) =  .  = 0 .  .  Fm (u) using the Newton’s Algorithm to find a better guess for u. • For a given u(k) , the vector F (u(k) ) can be computed using equations (27)-(31). Thus, we have   (k) x(t1 ; t0 , A) − u1   (k) F1 (u )   (k) (k) x(t2 ; t1 , u1 ) − u2   F2 (u(k) )       . (k)   .. F (u ) =  = ..     . (k) (k)   x(tm−1 ; tm−2 , um−2 ) − um−1  Fm (u(k) ) (k) x(tm ; tm−1 , um−1 ) − B Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

• The elements of the Jacobian matrix F 0 (u(k) ) can be computed using finite differences: ∂Fj Fj (u(k) + ∆el ) − F (u(k) ) = ∂ul ∆ where el is the l-th unit vector in RK and K = length(u1 ) + . . . + length(um ). • The value of Fj (u(k) + ∆el ) is obtained by solve the IVP x˙ = f (t, x),

(k)

x(tj ) = uj + ∆elj , tj ≤ t ≤ tj+1 , j = 0, . . . , m − 1.

• Once the values F (u(k) ) and F 0 (u(k) ) are available, solve the system of linear equaitons F 0 (u(k) )d = −F (u(k) ) to obtain a solution d(k) , and update the guess vector as u(k+1) = u(k) + d(k) . Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

• In fact, the Jacobian F 0 (u(k) ) has  −I X2 −I   X3  0 (k) F (u ) =    

the special structure  −I .. .

..

. Xm−1

    ,   −I  Xm

where (the following are sensitivities w.r.t. the initial values) ∂x(t2 ; t1 , u1 ) X2 (t2 ; t1 , u1 ) = ∂u (k)

(32)

u=u

.. .

(33) ∂x(tm−1 ; tm−2 , um−2 ) Xm−1 (tm−1 ; tm−2 , um−2 ) = (34) ∂u (k) u=u ∂x(tm ; tm−1 , um−1 ) Xm (tm ; tm−1 , um−1 ) = (35) ∂u (k) Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations u=u TU Ilmenau

A Multiple-Shooting Algorithm • Select shooting-points t0 < t1 < . . . , tm = tf , initial guess (0) (0) (0) (u1 , u2 , . . . , um ), ∆ > 0 and tolerance ε > 0 while kF (u(k) )k > ε do (k) • Solve the IVPs in equation (21) twice to obtain x(tj+1 ; tj , uj ) (k)

and X(tj+1 ; tj , uj ), j = 0, 1, . . . , m − 1 • Construct the vector the vector F (u(k) ) and the Jacobian matrix F 0 (u(k) ). • Solve the system of linear equations F 0 (u(k) )d = −F (u(k) ) to obtain dk as a solution. • Set k ← k + 1 • Set u(k) ← u(k) + dk end while Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multiple-Shooting ... Remarks Remark (Parallel Computation) I The most time-consuming step in the multiple-shooting method is the computation of the values of the vector F (u) and Jacobian matrix F 0 (u) for each Newton iteration. I The special structure of the Jacobian F 0 (u) (also of F (u)) indicates that: • each j-th block row of the Jacobian can be obtained through the solution of the IVP only on the j-th shooting interval. ⇒ the computation of the Jacobian can be done by solving smaller IVPs in parallel, correspond to the block rows of the Jacobian.

Choosing the shooting points I The optimal choice of the shooting points t1 , t2 , . . . , tm−1 can be problem dependent. Hence, they can be considered as variables to be determined through some computational scheme. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

Multiple-Shooting ... Advantages/Disadvantages Advantages: • it is efficient when the time-interval [t0 , tf ] too long; • it is easy to parallelize and, hence, applicable to problems of large dimensions; • accuracy of computations can be improved by considering additional intermediate points on the interval [tj , tj+1 ]; • it is applicable to solve stiff differential equations and DAEs Disadvantage: • For large-number of sub-intervals it can be computationally intensive; • The Newton-step may require the solution of a large set of linear equations requiring special linear-algebra strategies, • The choice of the shooting time points t1 , t2 , . . . , tm−1 has an impact on the quality of the solution. Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

References and Toolboxes References: • U. M. Ascherm, L. R. Petzold: Computer methods for ordinary differential equations and differential algebraic equations. SISM 1998. • L. F. Shampine and M. W. Reichelt: The Matlab ODE Suite, SIAM J. Sci. Comput., V. 18, pp. 1-22, 1997. • R. Ashino, M. Nagase: Behind and beyond the Matlab ODE Toolbox. Comput. Math. Appl. V. 40, pp. 491 – 512, 2000. • L. F. Shampine, J. Kierzenak:. Solving boundary value problems for ordinary differential equations in Matlab bvp4c. • W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery: Numerical recipes in C++ (3r ed.). Cambridge University Press, 2007.

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau

References and Toolboxes... Toolboxes: The Matlab ODE-Toolbox found under the MATrix LABoratory (MATLAB) software - The Math Works. SUite of Nonlinear and DIfferential/ALgebraic equation Solvers (Sundials) : https://computation.llnl.gov/casc/sundials/main.html Modelica : https://modelica.org/ Scialb: http://www.scilab.org/ EMSO (Environment for Modelling, Simulation and Optimisation): http://www.vrtech.com.br/rps/emso.html

Lecture 3 Introduction to Numerical Methods for Differential and Differential Algebraic Equations

TU Ilmenau