1. Solve a Large, Transient System Using LU Decomposition.
In the attached fileTransientHeatFlux_n100.matare an A matrix and b vector right-hand-side representing a model of a plate held initially at 100C, with edges at 0C [note that the model was developed using the finite difference method, but this is not important for your solution]. Solution of the matrix system is equivalent to taking a step in time; that is, the b vector may be seen as b(0), and the solution of A x = b is really solving A b(1) = b(0).
Use the MATLAB function lu() to perform LU decomposition of the A matrix, then perform 50 solves (i.e. 50 time steps), computing the b vector at each new time from the previous b vector. Use the MATLAB backslash function (\) to perform matrix solves with the L and U matrices provided from lu(). Overwrite the old b vector (don't save them).
At each timestep, visualize the current solution (b) using the following commands:
axis manual; axis([0 1 0 1 0 100]);
drawnow;
bb = reshape(b,N,N);
surf(xx,yy,bb);
Note that N, xx, and yy are also provided in the .mat file above. Given that the u vector in the .mat file is u(0), print out the surface shape (into your submitted PDF) at the 10th, 20th, and 30th timestep.
Comment 03/10: You've been given an A matrix and b vector that describe a physical system. The xx and yy variables only describe the grid - they are not used in solution. Further, when you solve the system, the solution (the "x" vector in lectures) is really b at the next time step. So to advance even further, you use that vector you just solved for as the new b vector, and solve again.
Comment 03/11: Units of the plot are °Celsius (height) and meters (width and length).
2. Solve a Large System Using Iterative Methods.
For the following two sub-problems, use the built-in simple MATLAB matrix algebra rules, i.e. if you're multiplying a matrix L by a vector xi, it's just L*xi.
- Code up the Jacobi Method as presented in lecture and use it to solve the large matrix of problem 3 from HW4.Return the iteration countto solve this problem and a plot of the solution.
- Modify the Jacobi method to use the most current information (i.e. implement the Gauss-Seidel method) and solve problem 3 of HW4 once more.Again include the iteration countand solution plot.
Regarding convergence, I suggested the norm(x1-x0)/norm(x0) criterion today. Note further that you should use a lower tolerance with this, like 1e-4 (in class I said 1e-2 would be sufficient, but I mis-spoke).
EXTRA CREDIT
For 20 points extra credit, in problem 1, use tic() and toc() - RTM if you don't know these - to
- Time the LU Decomposition;
- Time the 50 solves using the L and U matrices by substitution;
- Time 50 solves using direct solution (i.e. Gauss elimination each time).
- Report the overall times, the per-solve time (i.e. the 50 solve times divided by 50), and comment.
Be sure to reset the initial vector (b) prior to the second solve, and, be ready for a long wait for the second set of solves (using Gauss elimination each time). This is why we prefer LU Decomposition.