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
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
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
Note that the exact solution of (3) is
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
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?