Poissons Equation Consider the linear system An p1 where An

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)

Solution

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.