python计算: Project 2: Shock

1 Introduction

Project 2: Shock

Figure 1: A sonic shock

When a plane travels at a speed higher than that of sound, i.e. at supersonic speed, a shock wave forms around it. This shock consists of a very fast variation in space of physical quantities like pressure. It is what gives rise to the known “sonic-boom” that can be heard when a supersonic plane flies over you. Shocks are intrinsically related to nonlinear waves and indeed waves breaking on the ocean are produced by a similar mechanism. The sound of an explosion is also an effect of a generated shock wave.

In this project we will compute viscous shock waves using a famous model problem for turbulence known as the “viscous Burger’s equation”. Since Burger’s equation is set in one space dimension we can apply the tools that we have seen in the course to compute approximate solutions

1

2 The problem

In this project the aim is to write a program that computes an approximate solution to the boundary value problem

−ε d2 u (x) + u(x) du (x) = 0 in (−L, L) (1) dx2 dx

u(−L) = a (2) u(L) = b (3)

where ε ∈ R+. This equation is known as the stationary form of the viscous Burger’s equation that was introduced as a model problem for turbulence. It is well known that its solution can exhibit shocks, that will converge to discontinuities in the solution as ε → 0.

Verify that for the case L = ∞ the exact solution is given by
u(x) = −tanh􏰳 x 􏰴. (4)

2ε
We can use this solution to impose the boundary conditions at −L and L.

2.1 Numerical discretisation

Following what we did in class, we approximate u : (0, 1) 􏰲→ R as a vector U ∈ RN so that for h = 1/N and i ∈ {0,…,N}, Ui ≈ u(xi) where xi = −L+ 2L∗i×h. Observe that the solution is known for i = 0 and i = N since they correspond to the boundary values U0 = u(−L) = a = −tanh(−L/(2ε)) and UN = u(L) = b = − tanh(L/(2ε)). Since these two values are known through (2)-(3).

We recall the finite difference approximations
d2u(xi)≈(D2U)i =(Ui+1−2Ui+Ui−1)/h2 i=1,…,N−1 (5)

dx2

and

du(xi)≈(DU)i =(Ui+1−Ui−1)/(2h), i=1,…,N−1. (6) dx

Using these approximations of the derivatives we may write an approximate (finite dimensional) problem:

−ε(D2U)i+Ui(DU)i=F, i=1,…,N−1 (7) U0 = u(−L) (8) UN = u(L). (9)

2

Observe that the system (7) corresponds to N − 1 linear equations and can therefore be written on the form

AU + C (U ; U ) = F (10) • U ∈ RN+1 is the approximate solution and F ∈ RN+1 is a data term

where
dependingonaandbsuchthatF0 =aandFN =b.

• A ∈ RN+1×N+1 is the system matrix depending on the parameter ε. The first line of A consists of 1 followed by N zeroes. Lines i = 1,…,N−1 are defined by the first term in the first line of (7) together with (5). Finally line N consists of first N zeroes and then 1.
• C : RN+1 􏰲→ RN+1 is the nonlinear function defined by C0 = CN = 0, Ci(V ; U) = Vi(DU)i, i = 1, . . . , N − 1.

2.2 Solving the equations

The system (10) represents a nonlinear system of equations. This problem can be solved by fixed point iteration. Observe that the grid spacing h must be taken smaller than ε, otherwise the grid can not resolve the solution and the iterative algorithm may fail to converge.

1. Fix a tolerance Tol.
2. SetU0 =0,κ=1.
3. Solve the linear system

   (A + C(Uκ−1))Uκ = F (11) where C(Uκ−1) is the matrix where the first and last rows are zero and

   thelinesi=1,…,N−1aredefinedby(C(Uκ−1)Uκ)i =Uκ−1(Uκ − i i+1

Uκ )/(2h). The linear system (11) can be solved using the iterative i−1

solver of homework 3, exercise 4, or using the Crout algorithm for the solution of a trilinear system (see below).

4. Compute the increment ∆ U κ =

1 􏰤N 􏰥2

N − 1 􏰞 | U iκ − U κ − 1 | 2 . i

i=0

5. If ∆Uκ >Tol, increase κ by one and return to point 3 above to continue iterations. Otherwise the iteration has converged and we leave the fixed point iteration loop.

3

3 The project

The aim of the present project is to use the elements above to write a program that computes approximations of the differential equation (1)-(3) using the discrete formulation (7)-(9). Below follows a point by point outline of how the program can be constructed and what it should deliver.

3.1 Program structure

1. Open the file shock.bvp for read.
2. Read the equation data ε, L. from the file.
3. Read the parameters of the computation, N, Nlev and Tol.
4. An example of the file (shock.bvp), showing how the data should be organised, reads (here ε = 0.01, L = 1, N = 100, Nlev = 4 and T ol = 0.01):

        0.01
        1
        100
        4
        0.01
   

5. Construct the nonlinear system (7)-(9) for the given data.
6. Solve the nonlinear system, using the fix point procedure described above. For the linear solver either use the iterative solver that you programmed in week 3 or a Crout factorization. A Python code for the Crout factorization is available for download on the moodle page.

3.2 Computation and output

1. Compute the approximation U ∈ R2kN−1 with h = 1/(2k N) for k = 0,…, Nlev-1.
2. for each approximate solution compute the approximate error ei = u(xi) − Ui for all nodes xi on the level. Then compute

emax =max|ei|, i

􏰤 􏰥1 2

el2 = h􏰞|ei|2 i

4

.

4

3.

Report the results in a table on the form: h el2 eoc(el2) emax eoc(emax)

followed by the values of the different quantities for the different mesh- sizes h. Here “eoc(e)” stands for ‘experimental order of convergence and is defined by

eoc(ek) = log(ek−1/ek) log(2)

where ek is the error on level k. Clearly the eoc can not be computed for the first computation, k = 0.

If the error behaves as e = O(hα), eoc is an approximation of α. You would expect to get eoc(el2)≈ 2. If you do not observe this a common problem is that the tolerance in the iterative solver(s) is too large.

Plot the exact solution given by (4) and the approximate solution in the interval (−L, L).

Finally observe that the boundary conditions (8) – (9) are artificial since we use that we know the exact solution. Repeat the computa- tions for one of the above cases using the realistic boundary conditions u(−L) = 1, u(L) = −1. Is (4) still a good approximation of the exact solution?

What to submit for the assessment

4. 5.

You should submit the following items in a zipped folder:

• A .pdf document with details on the members of the group, names and student numbers. Here you may also give some information on the project, difficulties encountered etc.
• A well commented code that produces the expected output as outlined in the discussion above.
• A document containing the output of the program run using the input file (available on the moodle page): shock.bvp.
• Include both text (i.e. the tables discussed in the previous section) and graphical output (graphics can be saved using the buttons on the bottom of the Figure window, or using the commands exemplified in the program matplotexample.py available on the moodle page).

  5