Math 450 Homework 3 - Spring 2009

Click here for Math450 homework page

HOMEWORK #3 due Wednesday March 4

Text references are to
   Introduction to partial differential equations with MATLAB
   Jeffrey Cooper (1998) (Birkhauser)

NOTE: In the following, ^ means superscript, _ (underscore) means subscript, and Sum(i=1,9) means the sum for i=1 to 9.

1.  Consider the system

u_t + c u_x = 0     x > 0,  t > 0         (1)
u(x,0) = f(x) (x > 0),   u(0,t) = g(t) (t > 0)

where c > 0, the functions f(x) and g(t) are continuous at zero, and f(0)=g(0)=1.

    (i) Find the characteristics of equation (1) parametrized by the location x_0 where the characteristic intercepts the X axis. Note that x_0 can be positive, negative, or zero.

    (ii) Show that x_0 > 0 if and only if x > ct, in which case u(x,t)=f(p_1(x,t)) where p_1(x,t) > 0 (and find p_1(x,t)), and x_0 < 0 if and only if x < ct, in which case u(x,t)=g(p_2(x,t)) where p_2(x,t) > 0, and find p_2(x,t). (Note that p_1(x,t) represents a point on the X-axis while p_2(x,t) represents a point on the T axis.)


2. (i) Do Section 3.3 Problem 4 parts (a,b) (p92).
    (ii) Do Section 3.3 Problem 10 (p95). (The two problems are related.)


3.  Do Section 3.4 Problem 10 (p100). (This derives an exact solution for the heat equation for x>0 with a Robin-type boundary condition at x=0.)


4.  Consider an example of Burger's equation with viscosity

u_t + u u_x = (1/4) u_xx     t > 0,   all x         (2)
u(x,0) = f(x)     all x

Consider the traveling-wave solution u(x,t)=f(x-ct) derived in Section 3.5 with lim(x to -infinity) f(x)=1 and lim(x to +infinity) f(x)=0.

    (i) Find the values of u(t/2,t) for t=1,2,3,4,5,6,7. Where do these values seem to be in reference to the moderated shock wave that should be coming in from the left?

    (ii) Plot the values of u(2,t) for 0 ≤ t ≤ 8. These show the position of the shock wave at the point x=2 for this range of t. In particular, what are the values when t=1,4,7? (For t=1,4,7, you can use the approximation exp(3)=20.)


5.  (Essentially Section 3.7 Project 1 (p110)). Write a MATLAB program that implements the fully-explicit Euler update scheme (3.33) in the text for the system

u_t = k u_xx     0 < x < 10,  t > 0         (3)
u(x,0) = f(x) = sin(pi*x/10)     0 < x < 10
u(0,t) = u(10,t) = 0

Note that the exact solution of (3) is

u(x,t) = sin(pi*x/10) exp(-k(pi/10)^2 t)         (4)

Consider the update scheme (3.33) in the text in terms of the parameter rr=2s=2k(delt)/(delx)^2 instead of in terms of s in the text. Then the stability condition for the scheme is rr le 1 instead of s le 1/2.

    (i) Set k=1 for simplicity and write a MATLAB program to generate plots of the estimated curves u(x,t) for t=0,5,10,15,20.
    (Hints: (a) It should be easier to modify the program cl.m from Chapter 2 (on the Math450 Web site) than the Crank-Nicolson program heat3.m from Chapter 4. Delete all of the initial `input` statements (for input from the user) except for the input values of delx then delt. Then set n1=n2=n3=n4=5/delt, so that t0,t1,t2,t3,t4 = 0,5,10,15,20. Have the program compute t1,t2,t3,t4 as in cl.m and then display n1,n2,n3,n4 and 0,t1,t2,t3,t4 to make sure that you and your program are on the same track. Make the initial values u=u(x,0)=f(x) for f(x) in equation (3) above. Write a MATLAB program file f.m returning y=f(x) in (3). Alternatively, the line

f = inline('sin(pi*x/10)');

within your program (before it is used) has the same effect as writing a program file f.m.
    (b) Replace x = -1:delx:6 in cl.m by x=0:delx:10, since the range of X values is now [0,10] instead of [-1,6].
    (c) In updates from u=u(x,t) to u=u(x,t+delt), substitute the fully-explicit Euler scheme (3.33) in the text for the ``upwind'' Euler scheme in cl.m in all four loops. Replace the statement j=[2,jmax] by j=[2,jmax-1] before the updates u(j)=..., since you don't want your program to reach outside the interval [0,10] on either the left or the right. Since u(0) and u(10) are not updated, they remain at 0. As in cl.m, save the estimated values of u(x,t) for t=0,5,10,15,20 in the arrays snap0,...,snap4, which cl.m then displays.)

    (iii) For delx=0.50, run your program with delt=0.125 so that rr=1, and again with delt=0.250 so that rr=2. How does the graphed output compare with the exact solution in equation (4) above? Which value of delt appears to give output that is closer to the exact solution (4)?
    (Note: MATLAB allows you to enter 1/2 for 0.50 and 1/8 for 0.125 in response to `input` statements.)

    (iv) After the last of the four loops updating u(j), have your program calculate and display the sum of the absolute values of the differences between u(j) and the exact solution in (4) above for t=t4. What is the sum of the absolute errors in both cases?
    (Hint: After the last loop, set ut4=u*exp(-(pi/10)^2*t4) and display SumAbsErr = sum(abs(u-ut4)). If the numerical scheme converges, this value should be small.)

    (v) Do steps (iii-iv) again with delx=0.25 and delt=1/32 (so that rr=1) and delt=1/16 (so that rr=2). What is the sum of the absolute errors in both cases? Do you get similar results? Are the errors smaller?

  • Top of this page