CS / Math 243



CS / Math 243 Homework set 4 (updated 10-20-08) Due October 23

1. (a) Using your results in problem 5 of homework set 3 implement a variable step size version of crank-nicholson in Matlab. Save your code as cranknic4var.m. To help you in doing this on my web site I have cranknic4.m ( a fixed step size crank-nicholson for example 4) and explicit4var.m ( a variable step size implementation of the simplest explicit method). You will want to combine features of both of these programs. Run your program with h = 0.02, tol= 0.01, and tf = 1.5. Follow the run with compare4 and turn in a print out of the resulting graph.

b) Try different values m with k = tf / m in cranknic4 (which uses a fixed step size). Determine the smallest value of m (within 10 is fine) so that the maximum error in cranknic4.m is less than less than the maximum error in your runs from part (a). This maximum error is the largest entry in the fourth graph (lower right corner) produced by compare4. Use tic and toc to determine which is fastest – the variable step size program in part (a) or the fixed step size program in part (b).

2. Run your cranknic4var with larger tol. If you make the tol large enough does the automatic step size selection fail. Describe why.

3. Stiff differential equations are equations that require very small time steps to insure stability when using explicit methods. Typically stiff equations have a portion of the solution that is quickly varying and another portion that is slowly varying. If an implicit method is used with a stiff equation the step size is only small in the region that is quickly varying. However for an explicit method stability issues force the step size to be small in regions where the solution is slowly varying. The explicit method then often takes an enormous number of steps. Try explicit4var with tol = 0.01 say. You will probably find that the method takes too long (or runs out of memory trying to save usol). If you are tired of waiting type control-C in Matlab to stop the program. Print out tcur and step and turn these values in. Describe the progress of the method.

4. Matlab has a variety of ODE solvers which are described on a separate handout. These can be used to solve parabolic (and hyperbolic) PDE’s. When a PDE is solved using a standard ode solver the solution technique is called the method of lines. I have included a file mol.m on my web site that does the method of lines for problem 4 (with h fixed at h = 0.02) You can choose one of Matlab’s ode solvers by uncommenting the desired method.. Type “help ode113.m” (or whatever method) for more information on each Matlab ode solver. The ode solvers need to define the right side of the (system of) ODE’s. I did this in program exp4prime.m (on my web site). In this case the right hand side is a finite difference approximation to u xx + g(t) where g(t) = sin(pi*x)*( L * sech(L*(t-1))^2 + pi^2*tanh(L*(t-1))) and u xx is approximated with our standard central difference scheme. Compare4 still works after running mol.m

a) With tol = .01 run the mol program with ode15s.m uncommented. Discuss the step size and error graphs. Also use tic and toc to get the time of the run.

b) Reduce the tol so that you get approximately the same maximum error (use compare4 to observe this maximum error) as the run in problem 1 (a) . Now use tic and toc to compare the run times. Who is faster – Matlab and its professionally written code or your cranknic4var? Discuss some possible reasons why this is true.

c) Try my mol program with Matlab’s ode113.m, which is an explicit method, using tol = .01. You will need to wait several minutes for the solution to finish. Use tic and toc to time the method. Which is faster – the implicit method or the explicit method? By how much? Briefly discuss the Matlab statistics output from ode15s.m and ode113.m. Try to find out more about this output by searching the Matlab documention (using help or searching for Matlab web site). Does this explain your times? Discuss.

d) Try one other method from Matlab’s set of methods for stiff equations. Choose tol for this method so that the method has about the same maximum error (as determined by looking at compare4) as in part 4b. Use tic and toc to determine if ode15s is faster or if the other method is faster.

5. Solving a nonlinear parabolic PDE: On the first homework assignment set we looked at the heat equation if thermal conductivity was a function of u (in fact this does happen in practice). The resulting differential equation was nonlinear:

[pic], where g(t,x) is a source term.

Consider this PDE with zero Dirichlet boundary conditions, k(u) = 2 + u and with g(t,x) chosen so that u(t,x) = sin(pi*x) * tanh( L* (t-1) ) is the true solution. This requires:

g(t,x) = (I used Maple to do this because I knew I would make a mistake by hand)

[pic]

To reduce the chance of typos one can also use Maple to generate Matlab code using:

> [pic]

> [pic]

> [pic]

This creates the following for g(t,x):

sin(pi * x) * L - sin(pi * x) * L * tanh(L * (t - 1)) ^ 2 + 0.2e1 * sin(pi * x) * pi ^ 2 * tanh(L * (t - 1)) + sin(pi * x) ^ 2 * pi ^ 2 * tanh(L * (t - 1)) ^ 2 - cos(pi * x) ^ 2 * pi ^ 2 * tanh(L * (t - 1)) ^ 2;

which I enhanced (changing a few ^’s to .^’s to allow calculation for a vector of x values):

sin(pi * x) * L - sin(pi * x) * L * tanh(L * (t - 1)) ^ 2 + 0.2e1 * sin(pi * x) * pi ^ 2 * tanh(L * (t - 1)) + sin(pi * x). ^ 2 * pi ^ 2 * tanh(L * (t - 1)) ^ 2 - cos(pi * x) .^ 2 * pi ^ 2 * tanh(L * (t - 1)) ^ 2;

a) Discretize the above PDE using the standard approximation for [pic]and the approximation [pic][pic]and use the discretization to create a system of ordinary differential equations (ODE’s) in time.

b) Find Jacobian matrix for the right hand side of the ODE system. Is the Jacobian matrix sparse? Describe the sparsity.

(c) Create a function similar to my ex4prime and mol (at ) except make them for the above PDE. Let h = .02, L = 50, tol = .01 and solve the above PDE using the method of lines over the time interval [0,1.5]. You may want to use my code above for g(t,x). Try ode45 and ode15s. Compare the run times and accuracy. (Note: compare4 should still work to get the accuracy).

(d) By default ode15s uses finite differences to approximate the Jacobian (which is needed for implicit solvers and nonlinear equations). Write a Matlab function FJac.m that will return the n-1 by n-1 matrix J = the exact Jacobian. The first line should be “function J = FJac(t,u)”. Change the parameters in the ode solvers to allow use of the true Jacobian rather than a finite difference approximation by setting

options = odeset('RelTol',tol,'AbsTol',tol,'Stats','on','Jacobian',@FJac);

Run ode15s and compare the run times with the result of the runs in part (c). If the run times are too quick for a meaningful time comparison reduce tol to a smaller value.

(e) Include information about the sparsity pattern of the Jacobian in ode15s and compare the run times again. You will need to use odeset to add JPattern information.

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download