How to write 1D matlab code to solve poission's equation by multigrid method

1 view (last 30 days)
Hello Friends, I am developing a code to solve 1D Poisson's equation in matlab by multigrid method. I am getting the answer but not accurately. Please, help me to overcome with this difficulties. Here is my code for two grid method....
%%clear command and workspace window
clc;
clear; close all;
tic;
%%parameters
Nx = 256; % no of grid on x axis
Lx = pi/2 ; % length of the system along x direction
tol1 = 1.e-12; % tolerence
tol2 = 1.e-8; % tolerence
error1 = 1;
error2 = error1; % error
dx = Lx/(Nx);
x = (0:dx:Lx)' ; % length of x axis
%%boundary condition
u = ones(Nx+1,1);
u(1) = 1; % top
u(end) = 1; % bottom
f = -u;
uold = u;
%%exact solution
u_exact = sin(x)+cos(x);
p = 0;
residual = zeros(Nx+1,1);
nc2 = length(residual(1:2:end));
nc4 = length(residual(1:4:end));
cor_res2 = zeros(nc2,1);
cor_res4 = zeros(nc4,1);
ec2 = zeros(nc2,1);
ec4 = zeros(nc4,1);
%%iterate till GS converges
itrc = 6;
% while max(error2(:))> tol2
for p = 1:5
% p = p+1;
f = -u;
uold = u;
for k = 1 : itrc
uold = u;
for i = 2:Nx
u(i) = 0.5.*( uold(i+1) + u(i-1)- dx^2.*f(i) );
end
end
%%residual on fine grid
for k =1:itrc
for i = 2:Nx
residual(i) = f(i)-( u(i-1) - 2.*u(i) + u(i+1) )./dx^2;
end
end
%%coarse resdual calculation
cor_res2(1) = 0.5.*residual(1)+0.25.*residual(2);
for i = 2:nc2-1
cor_res2(i) = 0.25.*residual(2.*i-1)+ 0.5.*residual(2*i)+ 0.25.*residual(2*i+1);
end
cor_res2(end) = 0.25.*residual(end-1)+0.5.*residual(end);
%%coarse grid error calculation
z = 0;
while max(error1(:))> tol1
z = z+1;
ecold = ec2;
for i = 2:nc2-1
ec2(i) = 0.5.*( ec2(i-1) + ec2(i+1) + dx^2 .* cor_res2(i) );
error1 = abs((ecold-ec2)./ec2);
end
if (error1<tol1)
break
end
z = z+1;
end
nc4 = length(1:2:cor_res2);
cor_res4 = zeros(nc4,1);
cor_res4(1) = 0.5.*cor_res2(1)+0.25.*cor_res2(2);
for i = 2:nc4-1
cor_res4(i) = 0.25.*cor_res2(2.*i-1)+ 0.5.*cor_res2(2*i)+ 0.25.*cor_res2(2*i+1);
end
cor_res4(end) = 0.25.*cor_res2(end-1)+0.5.*cor_res2(end);
%%coarse grid error calculation
q = 0;
while max(error1(:))> tol1
q = q+1;
e4cold = ec4;
for i = 2:nc4-1
ec4(i) = 0.5.*( ec4(i-1) + ec4(i+1) + dx^2 .* cor_res4(i) );
error2 = abs((e4cold-ec4)./ec4);
end
if (error1<tol1)
break
end
q = q+1;
end
err_true = u_exact-u;
%%fine grid error by interpolation
err4 = interp1( 1:4:Nx+1, ec4, 1:Nx+1,'linear' ); % error on fine grid
err2 = interp1( 1:2:Nx+1, ec2, 1:Nx+1,'linear' ); % error on fine grid
u = u - err2'-err4'; % correction
for k = 1:itrc
uold2 = u;
for i = 2:Nx
u(i) = 0.5.*( uold2(i+1) + u(i-1)- dx^2.*f(i) );
error2 = abs((uold2-u)./u);
end
end
if (error2<tol2)
break
end
end
plot(x,u,x,u_exact,x,error2);
% plot(residual);
% plot(err);
title({'Two point GS';['Number of iterations: ',num2str(p)]});
toc;

Answers (0)

Categories

Find more on Polar Plots in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!