Poisson\'s Equation Consider the linear system A_n = p_1 where An is an n x n ma
ID: 3927922 • Letter: P
Question
Poisson's Equation Consider the linear system A_n = p_1 where An is an n x n matrix with 2's on the main diagonal, 1's directly above and below the main diagonal and O's everywhere else. For instance, A_s is A_s [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/partial differential x^2 = p. This equation appears very often in physics. Construct the matrix A_50 in Matlab. Try reading the documentation for diag(). Make the vector p according to the formula rho-j = 2 (1 - cos (23 pi j/51)) sin (23 pi j/51). (a) Write down the matrix form of the Jacobi iteration = k + c. Concatenate the matrix M and the vector c and save the resulting 50 times 51 matrix as A1.dat. (b) Use Jacobi iteration to solve for given an initial guess of a column of ones. Continue to iterate the Jacobi method until every term in the vector is within 10^-4 of the previous iteration. I.e., norm(phi(:, k+l) - phi(:, k), Inf)Explanation / Answer
Answer:
See the code below to create matrix "An" and vector "rho":
-------------------------------------------------
n=50; %dimensions of matrix
%elements of main diagonal
main_diag_elems=2*ones(1,n);
%create diagonal matrix
An=diag(elems);
%set elements above main diagonal
An((n+1):(n+1):(n*n-1))=-1;
%set elements below main diagonal
An((2):(n+1):(n*n-1))=-1;
%calculation of vector rho
j=1:n;
rho=2*(1-cos((23*pi)/51))*sin(((23*pi).*j)/51);
------------------------------------------------
a) See the code below:
-------------------------------------------
%Jacobi form in matrix notation
%phi=M*phi'+c; where phi is a 1*n order vector
%and c is also a 1*n order vector
phi=ones(1,n); %initial value of phi
C=5; %value of constant. You can specify accoring to yous.
c=5*ones(1,n); %vector of constants
%assuming M to be An, as M not specified.
%You can set or create M as required.
M=An;
%Concatenating M to c
Mc=[M,c'];
dlmwrite("A1.dat",Mc);
-----------------------------------
Note: Look at and do accordingly as per inline comments.
b) See the code below:
----------------------------------------
%Jacobi iteration
k=1; %first iteration
phi_k=phi; %value of phi for a given iteration
while(1)
phi=M*phi_k'+c;
if(norm((phi-phi_k),Inf)<=1e-4)
break;
end
phi_k=phi;
k=k+1;
end
%final iteration as column vector
dlmwrite("A2.dat",phi');
%total no. of iterations
dlmwrite("A3.dat",k);
------------------------------------------------
c), d) These can also be performed in the same way as shown under a) and b).
Note: M needs to be specified clearly what it is and how to calculate.