First Order Partial Differential Equations

��

�� First order partial differential equations in chemical engineering arise in models of processes that are dominated by convection rather than by dissipative processes. For example, consider a process that can be described by the 1-D convective diffusion equation with chemical reaction ∂c ∂t

∂c

+V

= ∂x

∂2 c ∂x2

+ f (c, k)

(1)

This process involves multiple time scales. There is a diffusive time scale: τD = L2 , a convective time scale: τC = L/ V, and a reaction time scale τR = 1 / k. Here k is the reaction rate constant, with units of inverse time. When diffusion and convection are important, the ratio of the respective time scales for those processes gives the Péclet number e = τD /τC = V L / . When reaction and convection are important, the corresponding dimensionless group is the convective Damköler number, a = τR /τC = V /L k . Thus in the limit of large Péclet numbers with finite convective Damköler number, say of order unity, the system is described by a first order PDE which when expressed in suitable dimensionless variables can be written as ∂c ∂c + Pe = f (c) (2) ∂t ∂x In these notes we discuss the solution of first order PDEs described by a

∂u ∂x

+b

∂u ∂x

=c

(3)

In particular we use the above equation to motivate the basic solution method we will used, called the method of characteristics.

�� Consider the following first order PDE a

∂u ∂x

+b

∂u ∂y

=c

(4)

where a, b, c are in general functions of x, y, and possibly also u(x,y) . Suppose we have found a solution to Eq. (4), then tuple {x, y, u(x, y)} defines a surface  in a Cartesian coordinate system, where the position vector r to any point on the surface is given by r (x, y) = x i + y j + u (x, y) k

(5)

From calculus we know that the unit normal to the surface is given by the following formula n=

∂r

∂r ⨯

∂x

∂r 

∂y

∂r ⨯

∂x

∂y

(6)

2 ��

PrimerMethodOfCharactristics.nb

A simple calculation shows that ∂r

∂r ⨯

∂x

∂y

= i+

∂u

k ⨯ j+

∂x

∂u ∂y

k =±

∂u ∂x

i+

∂u ∂y

j-k

(7)

Thus if u(x,y) is a solution to Eq. (4), then the vector p = (a i + b j + c k) must be orthogonal to the surface normal as  ∂∂xu i +

p·n ≡ p·

∂r ∂x

∂u ∂y

⨯

j - k

∂r ∂y

=a

∂u ∂x

+b

∂u ∂y

-c = 0

(8)

Thus we can conclude that the vector p must lie in the tangent plane of the solution surface  . So our task to finding a solution to Eq. (4) is equivalent to finding a surface  such that at every point on the surface the vector p lies in the tangent plane. Now consider a curve  that lies in the surface . Let its position vector r(x, y) be defined as in Eq. (5). It will be convenient to parameterize the curve in terms of arc length s along the curve, i.e., r (x (s), y (s)) = x (s) i + y (s) j + u[x (s), y (s)] k

(9)

Then the tangent vector local to the curve  is given by ⅆr

ⅆx =

ⅆs

ⅆs

i+

ⅆy ⅆs

j+

ⅆu ⅆs

k

(10)

Hence if our vector p = (a i + b j + c k) is going to be in the tangent plane it must also be every tangent to the curve . Thus the following parametric equations must hold ⅆx ⅆy ⅆu = a[x (s), y (s)], = b[x (s), y (s)], = c[x (s), y (s)] (11) ⅆs ⅆs ⅆs Hence by solving this system of ODEs, we are assured that p is tangent to the curve , which in turn lies in the solution surface . The curve  is called an integral curve (aka a solution curve) of the vector field p = (a i + b j + c k) . Integral curves are also called characteristic space curves of the PDE (Ref 4). The curve parameterized by {x(s), y(s)} along which x' (s) = a, y' (s) = b holds is called a characteristic ground curve for the PDE. Sometimes space curves and ground curves are used interchangeably and called characteristics. We will consider several example calculations.

�� Suppose we have the following PDE ∂u ∂t

+C

∂u ∂x

=0

(12)

where C is a constant. To solve this equation, we introduce the parameter s along the curve such that the origin of the curve satisfies x = x0 , t = t0 at s = 0

(13)

From Eq. (11) we know that a solution curve  satisfies the following ODEs ⅆx ⅆs

= C,

ⅆt ⅆs

= 1,

ⅆu ⅆs

=0

(14)

Integrating the above equations gives x = C s + K1 ,

t = s + K2 , u = K3

(15)

We can evaluate the constants of integration K1 and K2 from the initial conditions Eq. (13) imposed on


��

3

the solution curve. Thus we have x = C s + x0 ,

t = s + t0 , u = K3

(16)

We can eliminate the parameter s to get x = C (t - t0 ) + x0 ⟹ x0 = x - C (t - t0 )

(17)

From the last equation in Eq. (14) we note that u(x, y) is a constant along the solution curve  ( since ⅆu/ⅆs = 0 ) This means that u(x, y) must be some function of the "initial condition" specified for that curve . Thus we can write u = F (x0 , t0 )

(18)

where F is an arbitrary function that depends on the solution curve.

�� The initial-value or Cauchy problem is define as follows a

∂u ∂x

+b

∂u ∂t

=c

(19)

subject to the initial condition IC :

u (x, 0) = f (x)

A more general IC can be constructed by defining a curve Γ parameterized by s such that       Γ: x = p (s), t = q (s), u = u (s), s1 < s < s2

(20)

(21)

The initial data curve Γ is subject to the constraint that ⅆp ⅆq (22)  - b  ≠ 0 along Γ ⅆs ⅆs   This ensures that {p (s), q (s)} is not a characteristic ground curve of PDE Eq. (4). When then refer to Γ as the initial data curve for the PDE Eq. (19). Thus our problem statement becomes ∂u ∂u +b =c a (23) ∂x ∂t a

subject to the initial condition IC :

   u (x, t) = u p (s), q (s) = f (s)

(24)

Thus the solution to our initial value problem is a surface with the following parametric representation:    x = x s, s, t = t s, s, u = u s, s (25) Here s is the parameter that establishes the curve Γ that defines the initial data for the problem. Along the solution curve s varies and is assigned the value s = 0 on Γ. Thus the surface u(x, t) is found by solving the following system of ODEs ⅆx   = a, x 0, s = p (s) ⅆs ⅆt ⅆs

= b,

  t 0, s = q (s)

(26)

4 ��


ⅆu ⅆs

= c,

  u 0, s = f (s)

Once we have solved this system of equations we have the solution in the form   u s, s = u s (x, t), s (x, t)

(27)

In the neighborhood of Γ, we have x = xs, s, and t = ts, s. We can then ask whether we can solve these two equations to find s = s(x, t) and s = s(x, t). This can be answered using the implicit function theorem, namely det

∂x ∂s ∂t ∂s

∂x  ∂s ∂t  ∂s

≠0

(28)

Expanding Eq. (28) we get ∂x ∂t ∂x ∂t ∂t ∂x = a  -b  ≠ 0 -   ∂s ∂s ∂s ∂s ∂s ∂s

(29)

This is nothing more than the constraint that the initial data curve must satisfy. If this constraint is not satisfied then the solution to the initial value problem may not exist or have multiple solutions. Suppose our initial data is given by (30) u (x, 0) = Cos (x), x ≥ 0 Thus the parameter s represents the positive x-axis at t = 0 and we denote values along the x-axis by x0

and our system of equations becomes ⅆx = a, x (0, x0 ) = x0 ⅆs ⅆt ⅆs ⅆu ⅆs

= b,

t (0, x0 ) = 0

= c,

u (0, x0 ) = f (x0 )

(31)

Suppose our initial data is given by (32) u (0, t) = g (t), t ≥ 0 Thus the parameter s represents the positive t-axis and we denote values along the t-axis by t0 and our system of equations becomes ⅆx = a, x (0, t0 ) = 0 ⅆs

ⅆt ⅆs ⅆu ⅆs

= b,

t (0, t0 ) = t0

= c,

u (0, t0 ) = g (t0 )

In the next examples we will show how these ideas can be applied to different example problems.

(33)


��

5

�� Let us define the following initial value problem (IVP) ∂u ∂t

+C

∂u ∂x

= 0, -∞ < x < ∞, t > 0 (34)

u (x, 0) = f (x) As in the previous example, the solution is found by solving the characteristic equations ⅆx ⅆs

ⅆt

= C,

ⅆs

ⅆu

= 1,

ⅆs

=0

(35)

We saw in the previous example that u(x, t) is constant along characteristic curves , and that the value of u is determined by the value of x0 given by the initial data. That is u (x0 , 0) = f (x0 ), where - ∞ < x0 < ∞

(36)

But we showed in Example 1 that x0 = x - C t. Thus the solution is u (x, t) = f (x - C t)

(37)

Note that at t = 0 the solution satisfies the IC.

�� Consider the following PDE that describes the function u(t, x): 2xt

∂u

∂u +

∂x

∂t

= u,

-∞ < x < ∞, t > 0

(38)

with the initial condition imposed on the x-axis u (0, x) = x

(39)

The characteristic equations are ⅆx ⅆs

= 2 x t,

ⅆt ⅆs

= 1,

ⅆu ⅆs

=u

(40)

subject to the initial conditions at s = 0 x (0, x0 ) = x0 , t (0, x0 ) = 0, u (0, x0 ) = x0

(41)

In this example the function u(t, x) is not constant along characteristic curves . The solution strategy is to first solve for t then use the result to eliminate t from the equation for x as follows ⅆt = 1 ⟹ t = s as t = 0 at s = 0. (42) ⅆs Then we have ⅆx ⅆs

2

= 2 x t = 2 x s ⟹ x = A ⅇs , with x = x0 at s = 0

(43)

Hence the parametric equation for x is x = x 0 ⅇt

2

Finally we solve for u along the characteristic to get

(44)

6 ��


ⅆu ⅆs

⟹ u = A ⅇs

=u

(45)

Evaluating at t = 0 (s = 0) we have u (0, x0 ) = x0 ⟹ A = x0

(46)

so that the general solution is 2

u (x, t) = x0 ⅇt = x ⅇ-t ⅇt = x ⅇt-t

2

(47)

Let us use Mathematica to check that this is indeed a solution 2

D[u[x, t], t] + 2 x t D[u[x, t], x] ⩵ u[x, t] /. u → #1 ⅇ#2-#2 & // Simplify True

�� In this example we consider an initial boundary value problem (IBVP) u2

∂u

∂u +

∂x

∂t

= 0,

x > 0, t > 0

(48)

subject to the following conditions u (0, x) =

x, x>0

(49)

u (t, 0) = 0, t > 0 This is called an initial-boundary value problem. The characteristic equations for the parametric functions xs, s, ts, s, us, s: ⅆx ⅆs

= u2 ,

ⅆt ⅆs

= 1,

ⅆu ⅆs

=0

(50)

subject to the initial conditions at s = 0 along the positive x-axis: x (0, x0 ) = x0 , t (0, x0 ) = 0, u (0, x0 ) =

x0

(51)

The last equation shows that along characteristics u is a constant that depends on the initial condition u (s, x0 ) = F (x0 ) =

x0

(52)

Since u is a constant we can integrate to find x ⅆx ⅆs

ⅆt

= u ⟹ x = u 2 s + x0 ,

ⅆs

=0 ⟹s=t

(53)

Thus x = u 2 t + x0

(54)

Substituting for x0 in Eq. (52) gives u (s, x0 ) =

x0 =

x - u2 t ≡ u (t, x)

(55)

This is an implicit equation for u. We can square both sides to get u2 = x - u2 t and then solving for u we get x u2 = 1+t

(56)

(57)


��

7

Note that this function automatically satisfies the BC at x = 0 Let us check that this is a solution u[x, t]2 D[u[x, t], x] + D[u[x, t], t] ⩵ 0 /. u →

#1

1/2

1 + #2

& // Simplify

True

�� Consider the following nonlinear first order PDE 1 ∂u x ∂x

∂u

+ u2

=u

∂y

(58)

with the initial data u = x on y = 0,

-∞ < x < ∞

(59)

Thus our initial data lies on the "curve" Γ along the x-axis. The parametric equations we need to solve for xs, s, ts, s, us, s ⅆx = ⅆs

1 x

ⅆy

,

ⅆs

ⅆu

= u2 ,

ⅆs

=u

(60)

subject to the initial conditions at s = 0 along the x-axis: x (0, x0 ) = x0 , y (0, x0 ) = 0, u (0, x0 ) = x0

(61)

Solving for x gives x ⅆx = ⅆs ⟹

x2 2

= s + K1

(62)

Applying the initial condition that x = x0 at s = 0 gives x2 2

= s+

x0

(63)

2

Solving the parametric equation for u and applying the initial condition at s=0 gives ⅆu ⅆs

= u ⟹ u = K2 ⅇs = x0 ⅇs

(64)

The parametric equation for y using the IC y = 0 at s = 0 on  gives ⅆy ⅆs

= x20 ⅇ2 s ⟹y =

x20 2

ⅇ2 s + K 3 =

x20 2

e2 s - 1

(65)

We can eliminate the arc length parameter s by using Eq. (63) and Eq. (64) to solve for s in terms of x and u in terms of s s= Thus we can write

x2 2

-

x20 2

,

Log (u) = Log (x0 ) + s

(66)

8 ��


x2 = 2 Log (u) + x20 - 2 Log (x0 ) y=

x20

u2

2

x20

(67)

-1

Consider next the characteristic curve  that passes through x0 = 1, y = 0. The parametric equations for this characteristic are x2 = 2 Log (u) + 1, y =

(u2 - 1)

(68)

2

Shown below is a plot of the characteristic that passes through x0 = 1 and y=0 ParametricPlot

2 Log[u] + 1 ,

u2 - 1

, {u, 1, 5}, PlotStyle → Thick, 2 Frame → True, FrameLabel → {Style["x", 16], Style["y", 16]}, RotateLabel → False, AspectRatio → 1, FrameTicksStyle -> Directive[Black, 16], FrameStyle -> Directive[Black, Thickness[0.0025]]

12 10 8 y

6 4 2 0 1.0

1.2

1.4

1.6

1.8

2.0

x We can use Mathematica's Eliminate function to eliminate x0 from the two parametric equations and thereby obtain to a relationship between x, y and u Eliminatex2 == 2 Log[u] + x20 - 2 Log[x0 ], y == u2 - 2 y ⩵ ⅇu

2

-x2 -2 y

x20

u2

2

x20

- 1 , x0 

u2

Note that we are unable to write u as an explicit function of x and y. Let us check that our implicit function is indeed a solution of the PDE. First, we define the following equation to represent the implicit function:


eqn = - 2 y + u[x, y]2 ⩵ ⅇ-x - 2 y + u[x, y]2 ⩵ ⅇ-x

2

2

-2 y+u[x,y]2

-2 y+u[x,y]2

��

9

u[x, y]2

u[x, y]2

Then we compute the derivative ∂u/∂x and ∂u/∂y: Dx = SolveD[eqn, x], u(1,0) [x, y] 2

u(1,0) [x, y] →

ⅇu[x,y] x u[x, y] - ⅇx

2

+2 y

2

2

+ ⅇu[x,y] + ⅇu[x,y] u[x, y]2



Dy = SolveD[eqn, y], u(0,1) [x, y] u(0,1) [x, y] →

- ⅇx u[x, y] - ⅇx

2

2

2

+2 y

+ ⅇu[x,y] u[x, y]2

+2 y

+ ⅇu[x,y] + ⅇu[x,y] u[x, y]2 

2

2



We substitute the values of the derivatives into the PDE and find that our implicit function of u(x,y) is indeed a solution. 1 x

u(1,0) [x, y] + u[x, y]2 u(0,1) [x, y] ⩵ u[x, y] /. Dx /. Dy // Simplify

{{True}}

�� In this example we consider the following problem ∂u

∂u +

∂t

∂x

=0

0 < x < ∞, t > 0

(69)

subject to the conditions BC : IC :

u (0, t) = t, u (x, 0) =

t>0 (70)

x2 0 ≤ x ≤ 2 x 2≤x

The characteristics for this problem are ⅆt ⅆs

= 1,

ⅆy ⅆs

=1,

ⅆu ⅆs

=0

(71)

Solving these equations gives x = s + K1 , t = s + K2 , u = K3

(72)

For each BC-IC we have a set of initial conditions for Eq. (71). The BC in Eq. (70) represents the t-axis at x=0. Thus we take s = t0 along the t-axis, and the initial conditions for xs, s, ts, s, us, s at s = 0 are x (0, t0 ) = 0, t (0, t0 ) = t0 , u (0, t0 ) = t0

(73)

Apply these ICs to Eq. (72) gives K1 = 0, K2 = t0 , K3 = t0 This means that the characteristic curves emanating from the t-axis at x = 0 are

(74)

10 ��


x = s, t = s + t0 , u = t0

(75)

Eliminating t0 gives u (x, t) = t - x

(76)

along all characteristics passing through x = 0, t ≥ 0 Consider next the initial conditions for characteristics emanating from the x-axis. The initial conditions for xs, s, ts, s, us, s at s = 0 are x (0, x0 ) = x0 , t (0, x0 ) = 0, u (0, x0 ) = x20 for 0 ≤ x0 < 2 u (0, x0 ) = x0

(77)

for 2 ≤ x0

Now for 0 ≤ x ≤ 2, and t=0, we have at s=0 x0 = K1 , K2 = 0, u = K3 = x20

(78)

Thus x = s + x0 , t = s, u = x20

(79)

This gives after eliminating x0 and s the solution on all characteristics passing through 0 ≤ x ≤ 2, t = 0 u (x, t) = (x - t)2

(80)

And finally, for x>2, we have x = s + x0 , t = s, u = x0

(81)

This gives after eliminating x0 and s the solution on all characteristics passing through x > 2, t = 0 u (x, t) = (x - t) Note that all characteristics have slope ⅆ x/ⅆ t = 1, and hence do not intersect. Here is a plot of the characteristics and the various zones. Also show is the line at t0 = 1.25

(82)


plt1 = Plot{x, x - 2}, {x, 0, 5}, PlotRange → {0, 4}, AspectRatio → 1, PlotStyle → {{Red, Thick}, {Blue, Thick}, {Magenta, Thick}}, Frame → True, FrameLabel → {Style["x", 16], Style["t", 16]}, RotateLabel → False, Epilog → TextStyleForm"u(0,x)=x2 ", FontSize -> 16, {1, .15}, Text[StyleForm["u(0,x)=x", FontSize → 16], {3, .15}], TextStyleForm"u(t,x)=(x-t)2 ", FontSize -> 16, {2.5, 1.5}, Text[StyleForm["u(t,x)=(x-t)", FontSize -> 16], {4, 1}], Text[StyleForm["u(t,x)=(t-x)", FontSize -> 16], {1, 2.5}], Dashing[{0.02, 0.02}], Line[{{0, 1.25}, {5, 1.25}}], AbsolutePointSize[8], Black, Point[{1.25, 1.25}], Point[{3.25, 1.25}], FrameTicksStyle -> Directive[Black, 16], FrameStyle -> Directive[Black, Thickness[0.0025]]

4

3 u(t,x)=(t-x) t 2 u(t,x)=(x-t)2 1

0

u(t,x)=(x-t)

u(0,x)=x2 0

1

u(0,x)=x 2

3

4

5

x It is a simple matter to find the solution at any time by solving the following set of equations Remove[u]

��

11

12 ��


u[x_, t_] := (t - x) /; t ≥ x; u[x_, t_] := (x - t) /; t < x && x > t + 2; u[x_, t_] := (x - t)2 /; t < x ≤ t + 2; plt1 = Plot{x, x - 2}, {x, 0, 5}, PlotRange → {0, 4}, AspectRatio → 1, PlotStyle → {{Red, Thick}, {Blue, Thick}, {Magenta, Thick}}, Frame → True, FrameLabel → {Style["x", 16], Style["t", 16]}, RotateLabel → False, Epilog → TextStyleForm"u(0,x)=x2 ", FontSize -> 16, {1, .15}, Text[StyleForm["u(0,x)=x", FontSize → 16], {3, .15}], TextStyleForm"u(t,x)=(x-t)2 ", FontSize -> 16, {2.5, 1.5}, Text[StyleForm["u(t,x)=(x-t)", FontSize -> 16], {4, 1}], Text[StyleForm["u(t,x)=(t-x)", FontSize -> 16], {1, 2.5}], Dashing[{0.02, 0.02}], Line[{{0, 1.25}, {5, 1.25}}], AbsolutePointSize[8], Black, Point[{1.25, 1.25}], Point[{3.25, 1.25}], FrameTicksStyle -> Directive[Black, 16], FrameStyle -> Directive[Black, Thickness[0.0025]]; plt2 = Plot[Evaluate[u[x, 1.25]], {x, 0, 5}, PlotStyle → {Gray, Thick}, Frame → True, FrameLabel → {"x", "u(x,t)"}, RotateLabel → False]; Show[ plt1, plt2]

4

3 u(t,x)=(t-x) t 2 u(t,x)=(x-t)2 1

0

u(t,x)=(x-t)

u(0,x)=x2 0

1

u(0,x)=x 2

3

4

5

x Note the black dots indicate when dependence of the solution on a particular family of characteristics changes. We can explain the behavior in the above plots by examining the characteristics that pass through the dashed line. As we increase x at t = 1.25 we move along the dashed line. The value of u is determined


��

13

by the characteristic that passes through (x* , 1.25). For x* < 1.25 we are in a zone that has characteristics emanating from the t-axis; thus u[x,t] decreases. When we reach x* = 1.25 we are at the boundary between characteristics that emanate from the t-axis versus those that emanate from the x axis (see green point). Any further increase in x* means that the value of u(x* ,t) is given by u(x, t) = (x - 1.25)2 . Thus u(x,t) increases quadratically with increasing x until x* = 1.25 + 2. Then the value of u(x,t) is given by u(x, t) = (x* - 1.25) and u(x,t) increases linearly with x. Here is a plot of these equations for t = {0, 0.5, 1.0, 1.5, .…, 4} Plot[Evaluate[Table[u[x, j], {j, 0, 4, .5}]], {x, 0, 8}, PlotStyle → {{Gray, Thick}}, Frame → True, FrameLabel → {Style["x", 16], Style["u(x,t)", 16]}, RotateLabel → False, FrameTicksStyle -> Directive[Black, 16], FrameStyle -> Directive[Black, Thickness[0.0025]]]

8 6 u(x,t) 4 2 0

0

2

4

6

8

x

�� In this example we consider the following problem xu

∂u ∂x

+t u

∂u = -xt ∂t

(83)

where the initial data is specified on the curve Γ in the t-x plane IC :

u (x, t) = 5

on Γ

(84)

where the curve Γ has the following parametric form:    1 p (s), q (s) = s,   s Here is a plot of the initial data curve in the x - t plane

(85)

14 ��


��

ParametricPlots,

1

, {s, 0.1, 4}, Frame → True, s FrameLabel → {Style["x", 16], Style["t", 16]}, PlotStyle → Red, Epilog → Text[Style["u=5", 16], {1, 1.5}] 2.5

t

2.0

u=5

1.5

��

1.0

0.5

0

1

2

3

4

x The characteristic equations for this system are given by ⅆx ⅆs

= x u,

ⅆt ⅆs

= t u,

ⅆu ⅆs

= -x t

(86)

subject to the initial conditions that on the curve Γ u=5

(87)

The initial conditions for xs, s , ts, s , us, s then follow from the parametric form for Γ: 1     x 0, s = s, t 0, s =  , u 0, s = 5 s Note that the RHS of each characteristic equation depends on two variables. Thus they cannot be integrated directly. However, if we combine the first two equations we can effectively eliminate u 1 ⅆx 1 ⅆt = x ⅆs t ⅆs

(88)

(89)

which can be integrated to give  ln (x) - ln (t) = K1 (s)

(90)

We can use the initial conditions to determine K1 . At s = 0, x = s, y = 1/s which gives  ln (s) - ln

1  s

 = K1 (s)

Thus we can write Eq. (90) as x 2 ln   = ln (s ) ⟹ t

(91)

2 x=s t

Let us plot these characteristics in the t-x plane

(92)


��

��

15

initPt[s_] := {s, 1 / s} Block{$DisplayFunction = Identity}, x plt1a = TablePlot , {x, initPt[s][[1]], 5}, PlotStyle → Blue, s2 PlotRange → {0, 20}, {s, 0.2, 2, .3}; 1 plt1b = ParametricPlots, , {s, 0.05, 5}, Frame → True, s FrameLabel → {Style["x", 16], Style["t", 16]}, PlotStyle → Red; Show[plt1b, plt1a, PlotRange → {0, 15}, AspectRatio -> 1] 15

t

10

��

5

0 0

1

2

3

4

5

x Now substituting this result into the characteristic equation for u, we get ⅆu ⅆs

2 = -s t2

(93)

And combining this equation with the characteristic equation for t we get 2 2 ⅆu s t ⅆu 2 == -s t ⟹u ⅆt ⅆt tu

(94)

Integrating we get 2  u2 + s t2 = K3 (s)

(95)

Applying the initial condition u=5 on Γ, gives 2 25 + s

1 2 s

= K3 ⟹ K3 = 26

(96)

2 Finally noting that s = x/t, we get the solution

u2 + x t = 26 ⟹ u (x, t) =

26 - x t for xt < 26

(97)

16 ��


�� In this section we illustrate how to solve these hyperbolic equations using a finite difference scheme. Again we consider the following problem ∂u ∂t

+a

∂u ∂x

=0

(98)

For our the finite difference method we define a grid with xi = i h and tj = j k. In the Lax-Wendroff method we use a Taylor series expansion to write u[i, j + 1] ≡ u[xi , tj + k] = u[i, j] + k

∂u ∂t

i,j

+

1 2

k2

∂2 u ∂t2

+ …

(99)

i,j

We then used Eq. (98) to eliminate the time derivatives in Eq. (99) as it follows from Eq. (98) that the spatial derivative and time derivative are related by ∂ ∂ = -a (100) ∂t ∂x Thus Eq. (99) can be expressed as u[i, j + 1] = u[i, j] - ka

∂u ∂x

i,j

+

1 2

k2 a2

∂2 u ∂x2

+ …

(101)

i,j

Finally we use central difference approximation ∂u ∂x

i,j

=

u[i + 1, j] - u[i - 1, j] 2h (102)

∂2 u ∂x2

= i,j

u[i + 1, j] - 2 u[i, j] + u[i - 1, j] h2

to replace the spatial derivatives in Eq. (101) . The finite difference formula is then u[i, j + 1] = 1 1 (a P) (1 + a p) u[i - 1, j] + 1 - a2 p2  u[i, j] - a p (1 - a p) u[i + 1, j] 2 2

(103)

where p = k/h. The method can be shown to be stable for 0 < a p ≤ 1. In the next section we apply the Lax-Wendroff method apply to Example Problem 6

�� In this section we illustrate how to solve these hyperbolic equations using a finite difference scheme. To apply our finite difference method we need to define a grid. Here we set k = 0.005 and h = 0.01 ��

k = 0.005; h = 0.01; p = k  h; a = 1; a p;

Next we need to define the IC and BC by specifying rules for the nodal values at t = 0 and x=0 . Note that the initial data has discontinuities.


��

��

17

Remove[u, plt]; u[i_ /; 0 < i h ≤ 2, 0] := i h2 u[i_ /; 2 < i h, 0] := i h u[0, j_] := j k

Next we define a function that implements the Lax-Wendroff method defined by Eq. (103). If you inspect this function, you will observe that it implements a recursive formula to find the nodal values at a give i,j node. To avoid a lengthy computation we use the trick of saving previous values of computed nodal values by using the construct u[_, j_] := u[i, j] =. ... ��

u[i_ /; i > 0, j_ /; j > 0] := u[i, j] = 1 1 (a p) (1 + a p) u[i - 1, j - 1] + 1 - a2 p2  u[i, j - 1] - a p (1 - a p) u[i + 1, j - 1] 2 2

To generate a solution we use Table and let i vary from 1-350 and j from 1-200. The results are then plotted using ListPlot and stored in the variable plt[j]. Because we are using a recursive formula, the method is limited to how large i and j can be. ��

Table[plt[j] = ListPlot[Table[{i h, u[i, j]}, {i, 1, 350}], PlotStyle → RGBColor[0, 0, 1], Joined → True , Frame → True, PlotRange → {0, 4.5}, FrameLabel → {Style["x", 16], Style["u(x,t)", 16]}], {j, 1, 200}];

Here are some representative solutions ��

Show[{plt[1], plt[30], plt[150], plt[200]}] 4

��

u(x,t)

3

2

1

0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

x The results show the propagation of the wave with the discontinuities generated by the initial conditions. If we compare these results with the exact results we see that the Lax-Wendroff method is capable of representing the discontinuities but with some smoothing due to the discretizing method.

�� The method of characteristics is covered in most textbooks on PDEs, though the detail of presentation might vary a lot. The following references were helpful in preparing these notes: ◼ H - K. Rhee, R. Aris & N. R. Amundson, First - Order Partial Differential Equations : Volume 1, Prentice Hall, 1986

18 ��


◼ C. Pozrikidis, Numerical Computation in Science and Engineering, Oxford University Press, 1998 ◼ C. Pozrikidis, Introduction to Theoretical and Computational Fluid Dynamics, Oxford University Press, 1997 ◼ D. A. Anderson, J.C. Tannehill and R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer, Hemisphere Publishing Company, 1984