The Method of Characteristics Applied to Quasi-Linear PDEs

The method of characteristics applied to quasi-linear PDEs 18.303 Linear Partial Differential Equations Fall 2006 1 Motivation [Oct 26, 2005] Most of the methods discussed in this course: separation of variables, Fourier Series, Green’s functions (later) can only be applied to linear PDEs. However, the method of characteristics can be applied to a form of nonlinear PDE. 1.1 Traffic flow Ref: Myint-U & Debnath §12.6 Consider the idealized flow of traffic along a one-lane highway. Let ρ (x, t) be the traffic density at (x, t). The total number of cars in x1 ≤ x ≤ x2 at time t is Z x2 N (t) = ρ (x, t) dx (1) x1 Assume the number of cars is conserved, i.e. no exits. Then the rate of change of the number of cars in x1 ≤ x ≤ x2 is given by dN dt = rate in at x1 − rate out at x2 = ρ (x1 , t) V (x1 , t) − ρ (x2 , t) V (x2 , t) Z x2 ∂ = − (ρV ) dx x1 ∂x where V (x, t) is the velocity of the cars at (x, t). Combining (1) and (2) gives  Z x2  ∂ρ ∂ + (ρV ) dx = 0 ∂t ∂x x1 1 (2) and since x1 , x2 are arbitrary, the integrand must be zero at all x, ∂ρ ∂ + (ρV ) = 0 ∂t ∂x (3) We assume, for simplicity, that velocity V depends on density ρ, via   ρ V (ρ) = c 1 − ρmax where c = max velocity, ρ = ρmax indicates a traffic jam (V = 0 since everyone is stopped), ρ = 0 indicates open road and cars travel at c, the speed limit (yeah right). The PDE (3) becomes   ∂ρ 2ρ ∂ρ +c 1− =0 (4) ∂t ρmax ∂x We introduce the following normalized variables u= ρ ρmax , t̃ = ct into the PDE (4) to obtain (dropping tildes), ut + (1 − 2u) ux = 0 (5) The PDE (5) is called quasi-linear because it is linear in the derivatives of u. It is NOT linear in u (x, t), though, and this will lead to interesting outcomes. 2 General first-order quasi-linear PDEs Ref: Guenther & Lee §2.1, Myint-U & Debnath §12.1, 12.2 The general form of quasi-linear PDEs is A ∂u ∂u +B =C ∂x ∂t (6) where A, B, C are functions of u, x, t. The initial condition u (x, 0) is specified at t = 0, u (x, 0) = f (x) (7) We will convert the PDE to a sequence of ODEs, drastically simplifying its solu­ tion. This general technique is known as the method of characteristics and is useful for finding analytic and numerical solutions. To solve the PDE (6), we note that (A, B, C) · (ux , ut , −1) = 0. 2 (8) Recall from vector calculus that the normal to the surface f (x, y, z) = 0 is ∇f . To make the analogy here, t replaces y, f (x, t, z) = u (x, t) − z and ∇f = (ut , ux , −1). Thus, a plot of z = u (x, t) gives the surface f (x, t, z) = 0. The vector (ux , ut , −1) is the normal to the solution surface z = u (x, t). From (8), the vector (A, B, C) is the tangent to this solution surface. The IC u (x, 0) = f (x) is a curve in the u − x plane. For any point on the initial curve, we follow the vector (A, B, C) to generate a curve on the solution surface, called a characteristic curve of the PDE. Once we find all the characteristic curves, we have a complete description of the solution u (x, t). 2.1 Method of characteristics We represent the characteristic curves parametrically, x = x (r; s) , t = t (r; s) , u = u (r; s) , where s labels where we start on the initial curve (i.e. the initial value of x at t = 0). The parameter r tells us how far along the characteristic curve. Thus (x, t, u) are now thought of as trajectories parametrized by r and s. The semi-colon indicates that s is a parameter to label different characteristic curves, while r governs the evolution of the solution along a particular characteristic. From the PDE (8), at each point (x, t), a particular tangent vector to the solution surface z = u (x, t) is (A (x, t, u) , B (x, t, u) , C (x, t, u)) . Given any curve (x (r; s) , t (r; s) , u (r; s)) parametrized by r (s acts as a label only), the tangent vector is   ∂x ∂t ∂u . , , ∂r ∂r ∂r For a general curve on the surface z = u (x, t), the tangent vector (A, B, C) will be different than the tangent vecto (xr , tr , ur ). However, we choose our curves (x (r; s) , t (r; s) , u (r; s)) so that they have tangents equal to (A, B, C), ∂x = A, ∂r ∂t = B, ∂r ∂u =C ∂r (9) where (A, B, C) depend on (x, t, u), in general. We have written partial derivatives to denote differentiation with respect to r, since x, t, u are functions of both r and s. However, since only derivatives in r are present in (9), these equations are ODEs! This has greatly simplified our solution method: we have reduced the solution of a PDE to solving a sequence of ODEs. 3 2 f(x) 1.5 1 0.5 0 −3 −2 −1 0 x 1 2 3 Figure 1: Plot of f (x). The ODEs (9) in conjunction with some initial conditions specified at r = 0. We are free to choose the value of r at t = 0; for simplicity we take r = 0 at t = 0. Thus t (0; s) = 0. Since x changes with r, we choose s to denote the initial value of x (r; s) along the x-axis (when t = 0) in the space-time domain. Thus the initial values (at r = 0) are x (0; s) = s, t (0; s) = 0, u (0; s) = f (s) . (10) 3 Example problem [Oct 28, 2005] Consider the following quasi-linear PDE, ∂u ∂u + (1 + cu) = 0, ∂t ∂x u (x, 0) = f (x) where c = ±1 and the initial condition f (x) is f (x) = (      1, |x| > 1 =  2 − |x| , |x| ≤ 1    1, x < −1 2 + x, −1 ≤ x ≤ 0 2 − x, 0 < x ≤ 1 1, x>1 The function f (x) is sketched in Figure 1. To find the parametric solution, we can write the PDE as   ∂u ∂u (1, 1 + cu, 0) · , , −1 = 0 ∂t ∂x Thus the parametric solution is defined by the ODEs dt = 1, dr dx = 1 + cu, dr 4 du =0 dr with initial conditions at r = 0, t = 0, x = s, u = u (x, 0) = u (s, 0) = f (s) . Integrating the ODEs and imposing the ICs gives t (r; s) = r u (r; s) = f (s) x (r; s) = (1 + cf (s)) r + s = (1 + cf (s)) t + s 3.1 Validity of solution and break-down (shock formation) To find the time ts and position xs when and where a shock first forms, we find the Jacobian: ! xr xs ∂ (x, t) = det J = ∂ (r, s) tr ts = ∂x ∂t ∂x ∂t − = 0 − (cf ′ (s) r + 1) = − (cf ′ (s) t + 1) ∂r ∂s ∂s ∂r Shocks occur (the solution breaks down) where J = 0, i.e. where t=− The first shock occurs at 1 (s) cf ′  1 ts = min − ′ cf (s)  In this course, we will not consider what happens after the shock. You can find more about this in §12.9 of Myint-U & Debnath. We now take cases for c = ±1. For c = 1, since min f ′ (s) = −1, we have ts = − 1 =1 min f ′ (s) Any of the characteristics where f ′ (s) = min f ′ (s) = −1 can be used to find the location of the shock at ts = 1. For e.g., with s = 1/2, the location of the shock at ts = 1 is       1 1 1 1 1+ = 1+ 2− 1 + = 3. xs = 1 + f 2 2 2 2 Any other value of s where f ′ (s) = −1 will give the same xs . 5 For c = −1, since max f ′ (s) = 1, we have ts = 1 =1 max f ′ (s) Any of the characteristics where f ′ (s) = max f ′ (s) = 1 can be used to find the location of the shock at ts = 1. For e.g., with s = −1/2, the location of the shock at ts = 1 is       1 1 1 1 xs = 1 − f − 1− = 1− 2− 1 − = −1. 2 2 2 2 Any other value of s where f ′ (s) = 1 will give the same xs . 3.2 Solution Method (plotting u(x,t)) Since r = t, we can rewrite the solution as being parametrized by time t and the marker s of the initial value of x: x (t; s) = (1 + cf (s)) t + s, u (; s) = f (s) We have written u (; s) to make clear that u depends only on the parameter s. In other words, u is constant along characteristics! To solve for the density u at a fixed time t = t0 , we (1) choose values for s, (2) compute x (t0 ; s), u (; s) at these s values and (3) plot u (; s) vs. x (t0 ; s). Since f (s) is piecewise linear in s (i.e. composed of lines), x is therefore piecewise linear in s, and hence at any given time, u = f (s) is piecewise linear in x. Thus, to find the solution, we just need to follow the positions of the intersections of the lines in f (s) (labeled by s = −1, 0, 1) in time. We then plot the positions of these intersections along with their corresponding u value in the u vs. x plane and connect the dots to obtain a plot of u (x, t). Note that for c = 1, the s = −1, 0, 1 characteristics are given by s = −1 : x = (1 + cf (−1)) t − 1 = 2t − 1 s = 0 : x = (1 + cf (0)) t + 0 = 3t s = 1 : x = (1 + cf (1)) t + 1 = 2t + 1 These are plotted in Figure 2. The following tables are useful as a plotting aid: 1 t= 2 s = −1 0 u= 1 2 3 x= 0 2 6 1 1 2 t = ts = 1 s = −1 0 1 u= 1 2 1 x= 1 3 3 A plot of u (x, 1/2) is made by plotting the three points (x, u) from the table for t = 1/2 and connecting the dots (see middle plot in Figure 3). Similarly, u (x, ts ) = u (x, 1) is plotted in the last plot of Figure 3. Repeating the above steps for c = −1, the s = −1, 0, 1 characteristics are given by s = −1 : x = (1 − f (−1)) t − 1 = −1 s = 0 : x = (1 − f (0)) t + 0 = −t s = 1 : x = (1 − f (1)) t + 1 = 1 These are plotted in Figure 4. We then construct the tables: 1 t= 2 t = ts = 1 s = −1 0 u= 1 2 x = -1 − 12 1 1 1 s = −1 0 1 u= 1 2 1 x = −1 −1 1 As before, plots of u (x, 1/2) and u (x, 1) are made by plotting the three points (x, u) from the tables and connecting the dots. See middle and bottom plots in Figure 5. Note that for c = 1 the wave front steepened, while for c = −1 the wave tail steepened. This is easy to understand by noting how the speed changes relative to the height u of the wave. When c = 1, the local wave speed 1 + u is larger for higher parts of the wave. Hence the crest catches up with the trough ahead of it, and the shock forms on the front of the wave. When c = −1, the local wave speed 1 − u is larger for higher parts of the wave; hence the tail catches up with the crest, and the shock forms on the back of the wave. 4 Solution to traffic flow problem [Oct 31, 2005] The traffic flow PDE (5) is ut + (1 − 2u) ux = 0 7 (11) t 1 0.5 0 −3 −2 −1 0 x 1 2 3 Figure 2: Plot of characteristics for c = 1. u(x,0) 2 1 u(x,0.5) 0 −3 2 −1 0 1 2 3 4 −2 −1 0 1 2 3 4 −2 −1 0 1 2 3 4 1 0 −3 2 u(x,1) −2 1 0 −3 x Figure 3: Plot of u(x, t0 ) with c = 1 for t0 = 0, 0.5 and 1. 8 t 1 0.5 0 −3 −2 −1 0 x 1 2 3 Figure 4: Plot of characteristics with c = −1. u(x,0) 2 1 u(x,0.5) 0 −4 2 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 1 0 −4 2 u(x,1) −3 1 0 −4 x Figure 5: Plot of u(x, t0 ) with c = −1 for t0 = 0, 0.5 and 1. 9 and has form (6) with (A, B, C) = (1 − 2u, 1, 0). The characteristic curves satisfy (9) and (10) ∂x = 1 − 2u, ∂r ∂t = 1, ∂r ∂u = 0, ∂r x (0) = s, t (0) = 0, u (0) = f (s) . Integrating gives the parametric equations t = r + c1 , u = c2 , x = (1 − 2u) r + c3 = (1 − 2c2 ) r + c3 Imposing the ICs gives c1 = 0, c2 = f (s), c3 = s, so that t = r, u = f (s) , x = (1 − 2f (s)) r + s = (1 − 2f (s)) t + s (12) We can now write x (t; s) = (1 − 2f (s)) t + s, u (; s) = f (s) Again, the traffic density u is constant along characteristics. Note that this would change if, for example, there was a source/sink term in the traffic flow equation (11), i.e. ut + (1 − 2u) ux = h(x, t, u) where h(x, t, u) models the traffic loss / gain to exists and on-ramps at various posi­ tions. 4.1 Example : Light traffic heading into heavier traffic Consider light traffic heading into heavy traffic, and model the initial density as   α, x≤0 
3 (13) u (x, 0) = f (x) = − α x + α, 0 ≤ x ≤ 1 4   3 , x≥1 4 where 0 ≤ α ≤ 3/4. The lightness of traffic is parametrized by α. We consider the case of light traffic α = 1/6 and moderate traffic α = 1/3. From (12), the characteristics are [DRAW]   (1 − 2α) t + s, s≤0  x= (1 − 2α − 2 (3/4 − α) s) t + s, 0 ≤ s ≤ 1   −t/2 + s, s≥1 10 For α = 1/6, we have x=    2 t + s, 3
2 − 76 s t + 3 − 12 t + s, s≤0 s, 0 ≤ s ≤ 1 s≥1    1 t + s, 3
1 5 − s t+ 3 6 1 − 2 t + s, s≤0 s, 0 ≤ s ≤ 1 s≥1   For α = 1/3, we have x=   Again, for fixed times t = t0 , plotting the solution amounts to choosing an appropriate range of values for s, in this case −2 ≤ s ≤ 2 would suffice, and then plotting the resulting points u (t0 , s) versus x (t0 , s) in the xu-plane. The transformation (r, s) → (x, t) is non-invertible if the determinant of the Ja­ cobian matrix is zero, ! ! xr xs 1 − 2f (s) −2f ′ (s) r + 1 ∂ (x, t) = det = det = 2f ′ (s) r − 1 = 0. ∂ (r, s) tr ts 1 0 (14) Solving for r and noting that t = r gives the time when the determinant becomes zero, 1 . (15) t=r= ′ 2f (s) Since times in this problem are positive t > 0, then shocks occur if f ′ (s) > 0 for some s. The first such time where shocks occur is tshock = 1 . 2 max {f ′ (s)} (16) In the example above, the time when a shock first occurs is given by substituting (13) into (16), 1 1
. tshock = = 2 max {f ′ (s)} 2 34 − α Thus, lighter traffic (smaller α) leads to shocks sooner! The position of the shock at tshock is given by 1 −α xshock = (1 − 2α) tshock = 23 . −α 4 11