Scientificc computing exam 2021 PDF

Preview

Summary

Sample SolutionChair of Scientific Computing in Computer Science Department of Informatics Technical University of MunichEsolutionPlace student sticker hereNote: - During the attendance check a sticker containing a unique code will be put on this exam. - This code contains a unique number that assoc...

Description

Chair of Scientific Computing in Computer Science Department of Informatics Technical University of Munich

Note:

Esolution

• During the attendance check a sticker containing a unique code will be put on this exam. • This code contains a unique number that associates this exam with your registration number. • This number is printed both next to the code and to the signature field in the attendance check list.

Place student sticker here

Date: Time:

Friday 26th February, 2021 08:00 – 09:30

tio

IN2005 / Endterm Michael Bader

lu

Exam: Examiner:

n

Scientific Computing 1

Working instructions

So

• This exam consists of 18 pages with a total of 5 problems. Please make sure now that you received a complete copy of the exam. • The total amount of achievable credits in this exam is 50 credits. • Detaching pages from the exam is prohibited. • Allowed resources:

– your own handwritten notes

e

– all printed or electronic scripts, textbooks or lecture material

Sa m pl

• The working time for the exam is 90 minutes plus 30 minutes submission time. • Subproblems marked by * can be solved without results of previous subproblems. • The small scale of boxes printed next to each subproblem indicates the maximum number of credits awarded for this subproblem. Do not write into these boxes! • 18 credits will be sufficient to pass the exam.

• Answers are only accepted if the solution approach is documented. Give a reason for each answer unless explicitly stated otherwise in the respective subproblem. • Do not write with red or green colors nor use pencils. • During the exam, no questions concerning the content of the exam questions will be answered. • If you think that a question or exercise text contains an error or ambiguity, then choose a correction or variant that allows you to complete the solution – state how you have corrected or interpreted the exercise.

General Disclaimer for all solution examples: the given solution text is not intended as a “perfect answer” to the exercise; it only outlines the main aspect of the solution in compact form. It is assumed that you look up further details in the lecture material where necessary.

Left room from

to

/

– Page 1 / 18 –

Early submission at

Problem 1

Hatching of Pigeons (8 credits)

A pair of pigeons decided to settle on a balcony and lay their eggs there. On the next day, the pigeons discovered that the balcony is visited by people and got stressed about it. We model the pigeons’ behavior by Markov chains. Suppose there are three states in our model: 1. The pigeons abandon their nest and do not return. 2. The pigeons are worried about the people and consider to abandon the nest. 3. The pigeons are not afraid of people and do not consider to abandon their nest.

(1)

n

The transition matrix of the corresponding Markov chain is given by:   1 5/16 0 P = 0 1/2 1/4  . 0 3/16 3/4

tio

The transition matrix allows us to find the probabilities of the statesfn+1 after n + 1 days when we know the probabilities after n days:    1 5/16 0 f1 (n) fn+1 = Pfn = 0 1/2 1/4  f2 (n) , f3 (n) 0 3/16 3/4

0

lu

where f1 (n), f2 (n), f3 (n) are the probabilities of the corresponding three states after n days.

a) * Sketch the diagram of the described Markov chain corresponding to transition matrix (1). Use Figure 1.1 to make your sketch.

So

1

Sa m pl

e

2

Figure 1.1: Sketch the diagram of the described Markov chain. Two empty figures are provided in case you are not satisfied with your first attempt and want to produce a second. In case you do, use the second empty figure and make sure you cross out your worse attempt.

See Figure 1.1 left.

– Page 2 / 18 –

b) * Compute the probability distribution fn in the long term limit n → ∞. Does it depend on the initial distribution? Explain your answer with reasons/arguments.

0 1

Suggested variants of answers:

2

1. The first state is absorbing. As soon as the system appears at this states, it stays in it forever. In the long term limit, the system eventually will get to the first state independent on the initial state. Therefore,  ⊤ we conclude that f (n → ∞) = 1 0 0 , and it does not depend on the initial distribution.

n

2. We assume that in the long term limit the probability distribution converges to a certain stationary solution. Therefore, we solve the following equation to find the stationary solution   1 5/16 0 0 1/2 1/4  fs = fs . 0 3/16 3/4

tio

Furthermore, the sum of fs components must be equal one. Then the stationary solution is fs =  ⊤ 1 0 0 . This solution corresponds to λ1 = 1 eigenvalue, all other eigenvalues as shown in c) are smaller than one. Therefore, we conclude that contributions of other eigenvectors decays in time and the long term solution does not depend on initial condition.

Sa m pl

e

So

lu

3. Matrix P is column-stochastic, therefore, it has at least one eigenvalueλ1 = 1 and the absolute values of all other eigenvalues are less than one. Indeed, fromdet(P − λI) = 0, we find that λ1 = 1 , λ2 = 7/8 , and λ3 = 3/8, where |λ2 | < 1 and |λ3 | < 1 . As a result, the long term limit solution does not depend on the initial state and will converge to the eigenvector corresponding toλ1 = 1. Therefore, the corresponding  ⊤ eigenvector f1 = 1 0 0 is exactly the long term limit we are looking for.

– Page 3 / 18 –

0 1 2 3 4

c) * On average, the squabs (baby pigeons) appear after 18 days. Find the probability that parent pigeons will not abandon the nest during this period, if at the first day they are at state 3. Hint: the eigenvalues and eigenvectors of the transition matrix are:  ⊤ λ1 = 1 fe1 = 1 0 0  ⊤ (2) λ2 = 7/8 fe2 = −5 2 3  ⊤ λ3 = 3/8 fe3 = −1 2 −1 −8 Furthermore, you can use the approximations λ218 = 0.1 and λ18 . 3 = 2 · 10

f0 =

3 X

n

We can write down the initial probabilities vector as a linear combination of eigenvectors fj : aj fej .

Then after 18 iterations fn+1 = Pfn we obtain: f18 =

3 X

aj λj18fej .

j=1

lu

We find the coefficients a1 , a2 , and a3 from the linear system

tio

j=1

a1 fe1 + a2 fe2 + a3 fe3 = f0 .

⊤ For our initial distribution vector f0 = 0 0 1 we find that a1 = 1, a2 = 1/4, and a3 = −1/4. Now, we find the probability that in 18 days the nest is abandoned:

So



18 18 f18 (1) = a1 λ18 1 fe1 (1) + a2 λ2 fe2 (2) + a3 λ3 fe3

= 1 · 118 · 1 −

1 · 4

 18 7 8

Sa m pl

e

We conclude that the squab will appear with probability 1 − f18 (1) ≈ 0.12.

– Page 4 / 18 –

·5+

1 · 4

 18 3 8

· 1 ≈ 0.88.

Problem 2

ODEs: Analysis and Numerical Methods (13 credits)

Consider an electrical circuit containing a resistor, an inductor, and a capacitor, as in Figure 2.1. Such circuits are called RLC series circuits. They can be modelled by a second-order, constant-coefficient ODE: L

d2 q(t ) dt 2

+R

dq(t ) 1 + q(t ) = E(t ), C dt

(3)

tio

n

where L is the inductance, R is the resistance, C is the capacitance, E(t ) is the variable electric potential of the electric source, and q(t ) is the unknown charge on the capacitor. As a reminder, the electric current passing through t) such a circuit is given by I(t ) = dq( . dt

lu

Figure 2.1: RLC serial circuit

dq(t ) = dt



So

a) * Transform equation (3) in an equivalent system of first-order ODEs, under the assumption thatL 6= 0. Write the system of ODEs in matrix-vector form. 0

1

1 − LC

− LR 

q(t ) +

q(t ) I(t )



Sa m pl

e

q(t ) =



– Page 5 / 18 –



0

E(t ) L



0 1 2

0 1

b) Physical constraints on real-life RLC series circuits are, of course, thatR > 0, L > 0 , and C > 0. Show that this is already a sufficient condition for the given ODE to have a stable critical point in the homogeneous case, i.e. E(t ) = 0.

2

Method 2: Computing the generic eigenvalues   λ det (λI − A )=0 ˙ ⇒  1

LC

λ1,2

R ± =− 2L

−1 λ + RL s 

  R 1 2  ˙  = λ + L λ + LC =0

R 2L

2

−

So

lu

  ⇒ ℜ(λ1,2 ) < 0 checking both cases: λ1,2 ∈ R ∧ λ1,2 ∈ C \ R

1 LC

n

Method 1: Using principal minors λ1 + λ2 = tr(A ) = − RL < 0 1 >0 λ1 · λ2 = det (A ) = LC   ⇒ ℜ(λ1,2 ) < 0 ∀λ1,2 ∈ C

tio

3

Moving forward, we will consider a particular case of our RLC series circuit, in the form of the initial value problem       0 0 0 1 , b= , q(0) = . (4) ˙ ) = Aq(t ) + b, A = q(t −18 −6 180 9

1

0=A ˙ qˆ + b ⇒ qˆ = −A−1 b

Sa m pl

2

c) * Find all critical points of the ODE in (4).

e

0

A

1 = 18

ˆq =

 

−1

10 0

  −6 −1 18

0

=

 1 −3

– Page 6 / 18 –

1

1 − 18

0



d) Apply the method of Heun (second order Runge-Kutta) to the ODE in (4) using a time stepτ . Express the result in the form q(n+1) = M(τ ) q(n) + d(τ ).

0 1 2

Method of Heun:

3 (n)

=q

+ τ f (t

q(n+1) = q(n) +

For our ODE: q(n+1) = q(n) + (n+1)

q

A2 =

τ2 2

108

,q ) (n)

f (t



(n)

+ f (t (n+1) , qII )





2

(n)

, Ab =

−9τ 2 + 1 54τ 2 − 18τ





+

τ2 2



180 −1080

⇒ 



Ab + τ b

−3τ 2 + τ q(n) + 9τ 2 − 6τ + 1



e Sa m pl – Page 7 / 18 –

5



6



Aq(n) + b + A q(n) + τ Aq(n) + b

A + τA + I q

18

(n) , qI )



+b

90τ 2 −540τ 2 + 180τ

So

=



2

τ

2

  −18 −6

(n+1)

q

=



τ

(n )

n

=q

4 (n )

tio

(n)



lu

(n) qI (n) qII

Problem 3

Numerical Methods for ODEs (8 credits)

Consider the direction fields representing different numerical schemes for a 1D differential equationdp /dt = f (t , p), as given on page 9.

1

The different numerical schemes are demonstrated in Figure 3.1. These schemes are used to approximate the solution at times tn = 0.5 · n . Fill in the empty fields (formula, number of the figure, implicit or explicit, one- or multi-step) in the table below.

2

Hints:

0

3

1. The approximate values are marked by green dots.

4

2. The bold arrows correspond to the necessary function f (t , p) evaluations to get the approximate value of p at the following time steps.

n

5

3. The green dashed lines are used to connect the bold arrows to the place where they are applied, in case they are computed at different points.

7

4. Multi-step methods may require several points as the initial conditions.

tio

6

8

Ralston: pn+1 = pn + τ

2

BDF-2: pn+1 =

3

4 p 3 n



1

f (tn , pn ) + 34 f tn + 32 τ , pn + 32τ f (tn , pn ) 4

− 13 pn−1 + 32 τ f (tn , pn )

Explicit Euler:

4

e

pn+1 = pn + τ f (tn , pn )

Second order modified Euler:



Sa m pl

pn+1 = pn + τ f tn + 12 τ , pn + 21 τ f (tn , pn )

5



So

1

Figure Nr.

Second order Adams-Bashforth:

Implicit or Explicit

One or Multi-step

8

explicit

one-step

7

explicit

multi-step

6

explicit

one-step

1

explicit

one-step

5

explicit

multi-step

2

implicit

one-step

4

explicit

one-step

3

implicit

one-step

lu

Nr. Method – Formula

pn+1 = pn + 32 τ f (tn , pn ) − 21τ f (tn − τ , pn−1 )

6

Crank-Nicolson (Trapezoidal rule):

pn+1 = pn + 12 τ (f (tn , pn ) + f (tn+1 , pn+1 ))

7

Heun:

pn+1 = pn + τ

8

1

f (tn , pn ) + 12 f (tn + τ , pn + τ f (tn , pn )) 2

Implicit Euler:



pn+1 = pn + τ f (tn + τ , pn+1 )

– Page 8 / 18 –

n to S e Sa Figure 3.1: Direction field for the system of ODEs.

– Page 9 / 18 –

Problem 4

1D Flow Equation (12 credits)

For a 1D fluid that flows with very small velocity u(x, t ), we can derive the following 1D PDE: ∂u ∂2u − C 2 = g(t ) ∂x ∂t

(5)

where C ∈ R+ denotes the diffusion constant of the fluid, g(t ) is the pressure gradient −∂ p /∂ x . 0 1

n

2

a) * Develop a numerical scheme to solve equation (5): use a first-order discretisation for the temporal derivative, a time-implicit, spatial second-order finite difference discretization for the spatial second-order derivative and a time-explicit evaluation of the function g(t ). You may further assume an equidistant discretization of the unit interval in space,x ∈ [0, 1], using N + 1 grid points and a resulting mesh size h := 1/N . The time step shall be denoted by τ, the discrete velocity values foru(t , x) by uin := u(n τ , ih), i = 0, ..., N , and the discrete representation of the function g(t ) by g n := g(n τ ).

uin+1 − uni

n+1 n+1 uin+1 + ui+1 −1 − 2u i = gn h2

Sa m pl

e

So

lu

τ

−C

tio

The first-order time-implicit scheme is given by the implicit Euler method, the second-order spatial discretization corresponds to the default three-point stencil for the 1D Laplacian:

– Page 10 / 18 –

b) Reformulate the discretised equation from a) for one particular time stepn → n + 1 in the form of a linear system of equations N X Aij ujn+1 = bin , i = 0, ...N, (6)

0

with matrix (Aij ) ∈ RN+1×N+1 . The term bin should contain all contributions from the previous time stepn . Assume homogeneous Dirichlet conditions at the boundaries, i.e. u0n+1 = 0, uNn+1 = 0. Make sure to give an exact definition of the matrix entries Aij and right hand side entries bi .n

3

Splitting the equation from above into information from time step n and n + 1 results in   1 2C C n+1 C n+1 1 − 2 ui − = g n + uni + 2 uin+1 − 2 ui+1 1 + h τ τ h h

      

h2

1

i = j and (i = 0 ∨ i = N)

0

otherwise.

lu

τ

Aij :=

tio

for all (inner and boundery) points i = 0, ..., N . The matrix A can thus be defined as  |j − i | = 1 and i 6= 0, N −hC2        1 + 2C if i = j and i 6= 0, N

Sa m pl

e

So

The right hand side bi n evolves at bi n := g n + τ1 uni for i = 1, ..., N − 1 and b0n = 0, bNn = 0.

– Page 11 / 18 –

n

j=0

1 2

4

0 1 2

c) Now consider the case that the pressure gradientg(t ) is equal to 0 . The homogeneous equation is then of the ∂u ∂2u form: − C ∂ x 2 = 0. ∂t As in the von Neumann stability analysis we assume that the numerical solution is of type: (n )

ui

3 4

= akn sin(π kih )

(7)

Derive an explicit formula for the coefficientak such that ui(n) solves the time-implicit scheme developed in a) and b). Hint: You may use the equality sin(A + B) + sin(A − B) = 2 sin(A ) cos(B).

n

Inserting the definition of ui(n) = a nk sin(π kih ) into the discrete update equation from b, and putting gn = 0 yields:   1 2C C n+1 C 1 sin( π k (i − 1)h) + − 2 ak akn+1 sin(π k (i )h) − 2 akn+1 sin(π k (i + 1)h) = aknsin(π kih ) + τ τ h2 h h

1

ak =

2C τ h2

.

(cos(π kh) − 1)

Sa m pl

e

So

lu

1−

tio

Dividing by akn sin(π kih ) results in an expression for ak :

– Page 12 / 18 –

d) Under which condition will the solution u decay? State whether and which constraints are imposed on the choice of τ and h .

0 1

The solution decays if it holds |ak | < 1. Using the monotony of cos(π kh ) , it is enough to consider the extreme ! cases cos(π kh) = ±1. We obtain !

ak

cos(π kh) =1

=

!

ak

cos(π kh)=−1

=

1 1 1+

4τ h2

Sa m pl

e

So

lu

tio

n

On the interval cos(π kh) ∈ (−1, 1), |ak | < 1 is consequently fulfilled for all choices ofτ and h. Our time-implicit algorithm is hence unconditionally stable.

– Page 13 / 18 –

2

Problem 5

Cruzeix-Raviart Elements for FEM (9 credits)

The 2-dimensional heat equation is given as: ∂ u = ∆u, ∂t where ∆u = 0

∂2 u ∂x2

+

(8)

∂2 u ∂y2

a) * Derive the weak form of Equation (8) and simplify it as much as possible.

1

n

multiplying with testfunction integrating

tio

∂ u = ∆u ∂t ∂ φi u = φi ∆u ∂t Z Z ∂ φi ∆u dΩ , ∀i φi u dΩ = ∂t ΩZ ZΩ ∂ φi u dΩ = − ∇φi · ∇u dΩ , ∀i ∂t Ω Ω

integration by parts; boundary=0

Sa m pl

e

So

lu

2

– Page 14 / 18 –

b) Using the Ritz-Galerkin approach, derive the linear system of equations. Mark the integrals corresponding to the mass and stiffness matrix a...