Help writing code in matlab Consider the linear system An ph

Help writing code in matlab

Consider the linear system A_n phi = rho, where A_n is an n times n matrix with 2\'s on the main diagonal, -1\'s directly above and below the main diagonal and 0\'s everywhere else. For instance, A_5 is A_5 = [2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2]. This is a discretized version of Poisson\'s equation: partial differential^2 phi/partial differential x^2 = rho. This equation appears very often in physics. Construct the matrix A_50 in Matlab. Try reading the documentation for diag(). Make the vector rho according to the formula rho_j = 2 (1 - cos (23 pi/51)) sin (23 pi j/51). Write down the matrix form of the Jacobi iteration phi_k + 1 = M_phi_k + c. Concatenate the matrix M and the vector c and save the resulting 50 times 51 matrix as A1.dat. Use Jacobi iteration to solve for phi given an initial guess of a column of ones. Continue to iterate the Jacobi method until every term in the vector phi is within 10^-4 of the previous iteration. I.e., norm(phi(:, k+1) - phi(:.k), Inf)

Solution

MATLAB CODE FOR QUESTION a) AND b)

clear all; %clearing the workspace
n = 50; % size of matrix
A50 = diag(2*ones(1,n))... %Main diagonal with entries 2
-diag(ones(1,n-1),1)... %Upper diagonal with entries -1
-diag(ones(1,n-1),-1); % Lower diagonal with entries -1
Pj = 2*(1-cos(23*pi/51))*sin(23*pi*(1:50)/51); % Creating the vector Pj
% Create variable to store unknowns
Phi = zeros(n,1); % column vector
% Decomposing the A matrix in to D,L and U
D = diag(diag(A50)); % D is the diagonal matrix of A50
L = tril(A50,-1); % Lower diagonal matrix of A50
U = triu(A50,1); % Upper diagonal matrix of A50
M = D\\(-L-U); % Calculating M in the given equation
c = D\\Pj\'; % Calculating c in the given equation
A1 = [M c]; % Creating the 50x51 matrix to save
save A1.dat A1 % saving the 50x51 matrix as A1.dat
ext = 1; % loop exiting flag/
while ext % loop will continue until ext becomes 0
PhiOld = Phi; % save current values to calculate error later
Phi = (M*Phi+c); % Jacobi iteration
if( norm(Phi - PhiOld) <= 10^-4) % Check for a given accuracy
save A2.dat Phi; % If yes save the value of Phi in A2.dat
save A3.dat ext; % Save the number of iterations
ext =0; % Exit the loop by making ext = 0
else % If Not increment the iteration count
ext = ext + 1; % incrementing the iteration index
end
end

MATLAB CODE FOR c) AND d)

clear all; %clearing the workspace
n = 50; % size of matrix
A50 = diag(2*ones(1,n))... %Main diagonal with entries 2
-diag(ones(1,n-1),1)... %Upper diagonal with entries -1
-diag(ones(1,n-1),-1); % Lower diagonal with entries -1
Pj = 2*(1-cos(23*pi/51))*sin(23*pi*(1:n)/51); % Creating the vector Pj
% Create variable to store unknowns
Phi = zeros(n,1); % column vector
% Decomposing the A matrix in to D,L and U
L = tril(A50); % Lower diagonal matrix of A50
U = triu(A50,1); % Upper diagonal matrix of A50
M = L\\(-U); % Calculating M in the given equation
c = L\\Pj\'; % Calculating c in the given equation
A4 = [M c]; % Creating the 50x51 mtrix to save
save A4.dat A4 % saving the 50x51 matrix as A4.dat
ext = 1; % loop exiting flag/
while ext % loop will continue until ext becomes 0
PhiOld = Phi; % save current values to calculate error later
Phi = (M*Phi+c); % Gauss-Seidel iteration
if( norm(Phi - PhiOld) <= 10^-4) % Check for a given accuracy
save A5.dat Phi; % If yes save the value of Phi in A5.dat
save A6.dat ext; % Save the number of iterations in A6.dat
ext =0; % Exit the loop by making ext = 0
else % If Not increment the iteration count
ext = ext + 1; % incrementing the iteration index
end
end

OUTPUT

The programs have no output, instead of print the result the program saving result as dat files. If you want to check the result run the program and load the dat files A1 - A6 and use disp(A1), disp(A2) etc to check the result