MATLAB for differential equation
Joe Bob
Mon lab: 4:30-6:50
Lab 3
Exercise 1
(a) Create function M-file for banded LU factorization
function [L,U] = luband(A,p) % LUBAND Banded LU factorization % Adaptation to LUFACT % Input: % A diagonally dominant square matrix % Output: % L,U unit lower triangular and upper triangular such that LU=A
n = length(A); L = eye(n); % ones on diagonal
% Gaussian Elimination for j = 1:n-1 a = min(j+p,n); for i = j+1:a L(i,j) = A(i,j)/A(j,j); % Row multiplier b = min(j+p-1,n); A(i,j:b) = A(i,j:b) - L(i,j)*A(j,j:b); end end
U = triu(A);
end (b) Invoke function in command window
>> A = [1 5 3 -1; 2 4 9 9; 1 1 -1 -3; 4 3 10 3] % declare matrix A
A =
1 5 3 -1
2 4 9 9
1 1 -1 -3
4 3 10 3
>> luband(A, 4) % call luband function from command window
ans =
1.0000 0 0 0
2.0000 1.0000 0 0
1.0000 0.6667 1.0000 0
4.0000 2.8333 1.7500 1.0000
Exercise 2
(a) Create script-file that runs luband function (Lab3_ex2.m)
A = 10*eye(6)+6*diag(ones(5,1),-1)+6*diag(ones(5,1),1); [Lband, Uband] = luband(A,2); [L, U] = lu(A);
disp('The original matrix A is ') disp(A)
disp('Using luband we have') disp('L = ') disp(Lband) disp('U = ') disp(Uband) disp('A = ') disp(Lband*Uband)
disp('Using the built in lu we have') disp('L = ') disp(L) disp('U = ') disp(U) disp('A = ') disp(L*U)
(b) Output in command window
>> Lab3_ex2
The original matrix A is
10 6 0 0 0 0
6 10 6 0 0 0
0 6 10 6 0 0
0 0 6 10 6 0
0 0 0 6 10 6
0 0 0 0 6 10
Using luband we have
L =
1.0000 0 0 0 0 0
0.6000 1.0000 0 0 0 0
0 0.9375 1.0000 0 0 0
0 0 1.3714 1.0000 0 0
0 0 0 3.3871 1.0000 0
0 0 0 0 -0.5813 1.0000
U =
10.0000 6.0000 0 0 0 0
0 6.4000 6.0000 0 0 0
0 0 4.3750 6.0000 0 0
0 0 0 1.7714 6.0000 0
0 0 0 0 -10.3226 6.0000
0 0 0 0 0 13.4875
A =
10 6 0 0 0 0
6 10 6 0 0 0
0 6 10 6 0 0
0 0 6 10 6 0
0 0 0 6 10 6
0 0 0 0 6 10
Using the built in lu we have
L =
1.0000 0 0 0 0 0
0.6000 1.0000 0 0 0 0
0 0.9375 0.7292 -0.2153 -0.3704 1.0000
0 0 1.0000 0 0 0
0 0 0 1.0000 0 0
0 0 0 0 1.0000 0
U =
10.0000 6.0000 0 0 0 0
0 6.4000 6.0000 0 0 0
0 0 6.0000 10.0000 6.0000 0
0 0 0 6.0000 10.0000 6.0000
0 0 0 0 6.0000 10.0000
0 0 0 0 0 4.9954
A =
10.0000 6.0000 0 0 0 0
6.0000 10.0000 6.0000 0 0 0
0 6.0000 10.0000 6.0000 0 0
0 0 6.0000 10.0000 6.0000 0
0 0 0 6.0000 10.0000 6.0000
0 0 0 0 6.0000 10.0000
Exercise 3
(a)
(b)
(c)
.
.
.
.
Exercise 4
.
etc.