Nx = input('Number of nodes in x-direction');Ny = input('Number of nodes in y-direction');dx = 1/(Nx - 1);dx2 = dx*dx;dy = 1/(Ny - 1);dy2 = dy*dy;x = 0:dx:1;y = 0:dy:1;for i = 1:Nx for j = 1:Ny phi_an(i,j) = 500*exp(-50*((1-x(i))^2 + y(j)^2)) + 100*x(i)*(1-y(j)); endendSphi = zeros(Nx,Ny);for i = 1:Nx for j = 1:Ny Sphi(i,j) = 50000*exp(-50*((1-x(i))^2 + y(j)^2))*(100*((1-x(i))^2 + y(j)^2)-2); endend% Defining solution vector with intial guess
phi = zeros(Nx,Ny);phi(1,1:Ny) = 500*exp(-50*(1+y.^2)); % Left BC
phi(Nx,1:Ny) = 100*(1-y) + 500*exp(-50*y.^2); % Right BC
phi(1:Nx,1) = 100*x + 500*exp(-50*((1 - x).^2)); % Bottom BC
phi(1:Nx,Ny) = 500*exp(-50*((1-x).^2 + 1)); % Top BC
% Setting up the link coefficients
ae = 1/dx2;aw = 1/dx2;an = 1/dy2;as = 1/dy2;ao = -(2/dx2 + 2/dy2);toc%% Jacobi Method
tic% Matrix to store new values
phin = zeros(Nx-2,Ny-2);%Setting up the iterations
mat = zeros(1,1000000);R = zeros(Nx-2,Ny-2);for n = 1:1000000 R2 = 0; for i = 2:Nx - 1 for j = 2:Ny - 1 phin(i-1,j-1) = (Sphi(i,j) - ae*phi(i-1,j) - aw*phi(i+1,j) - an*phi(i,j+1) - as*phi(i,j-1))/ao; end end for r = 2:Nx-1 for s = 2:Ny-1 R(r-1,s-1) = Sphi(r,s) - ae*phi(r-1,s) - aw*phi(r+1,s) - an*phi(r,s+1) - as*phi(r,s-1) - ao*phi(r,s); R2 = R2 + R(r-1,s-1)*R(r-1,s-1); end end R2 = sqrt(R2); mat(n) = R2; if R2<10^-6 break end phi(2:Nx-1,2:Ny-1) = phin(1:Nx-2,1:Ny-2); end[X,Y] = meshgrid(x,y);contour(X,Y,phi);
The code to the above system of equations is the one above the image.
Best Answer