I'm just trying to solve a LINE SOR (PDE EQUATION), but I keep getting this error ??? Attempted to access d(50); index out of bounds because numel(d)=49.
Error in ==> thomas at 7 d(i) = (d(i) - a(i)*d(i-1))/(b(i) - a(i)*c(i-1));
Error in ==> laplacsorlinep1 at 96 ukp1(j,2:Nx-1) = thomas(a,b,c,d)
%---------------------------PROGRAM------------------------------------------------%
%
if true
end-------------------------------------------------------------------
clear
clc
% define variables %-------------------------------------------------------------------% dx = 0.1; %define dx B = 5; %define number of grides B H = 5; %define number of grides Y a = 0:dx:H; %define grides in x direction b = 0:dx:B; %define grides in y direction data=zeros(length(a),length(b)); %boundary condition %-------------------------------------------------------------------% for j = 1 for i = 1 : size(data,1) data(i,j) = 0; %boundary condition for west side end end for j = size(data,2) for i = 1 : size(data,1) if i < (1.8/dx) data(i,j) = 0; %boundary condition for east side elseif i > (2/dx) data(i,j) = 100; else i >= (1.8/dx) & i <= (2/dx) data(i,j) = 50; end end end for i = 1 for j = 2 : size(data,2) data(i,j) = 0; %boundary condition for top side end end for i = size(data,1) for j = 1 : size(data,2)-1 if j < (1/dx) data(i,j) = 0; %data(i,j-1) + 1; %boundary condition for bottom side elseif j > (1.2/dx) data(i,j) = 100; else j >= (1/dx) & j <= (1.2/dx) data(i,j) = 50; end end end
%-------------------------------------------------------------------% % define variables %data=xlsread('data.xlsx'); Nx = size(data,2); % no nodes in the x direction Ny = size(data,1); % no nodes in the y direction tol = 1d-6; % tolerance err = 1; % error k=0; % interation counter omega = 1; %relaxation parameter %-------------------------------------------------------------------% % define u array u = data; ukp1 = u; %-------------------------------------------------------------------% %calculate relaxation parameter for SOR method while err> tol D = 4*eye(Nx-2); L = -omega*diag(ones(Nx-3,1),-1); U = -omega*diag(ones(Nx-3,1),1); T = inv(D)*(L + U); newomega = 2/(1 + sqrt(1 - max(abs(eig(T)))^2)); err = abs(newomega - omega); omega = newomega; end err = 1; %calculate lower, main and upper diagonal vectors a(2:Nx-2) = -omega; b(1:Nx-2) = 4; c(1:Nx-3) = -omega;
%output Omega fprintf('omega = %1.4f\n',omega)
% line SOR method while err > tol % update iteration counter k = k + 1; fprintf(' %6i |',k); for j = 2 : Ny - 1 %calculate RBS vector for i = 2 : Nx - 1 d(i-1) = 4*(1 - omega)*u(j,i) + omega*ukp1(j-1,i) + omega*u(j+1,i); end d(1) = d(1) + omega*ukp1(j,1); d(Nx-2) = d(Nx-2) + omega*ukp1(j,Nx); %solve tridiagonal system using thomas algoritm ukp1(j,2:Nx-1) = thomas(a,b,c,d);
for i = 2 : Nx - 1
fprintf('%10.6f &',ukp1(j,i))
end
err = sqrt(sum(sum((ukp1 - u).^2)));
fprintf(' %8.6f \n',err)
u = ukp1;
end
end
imagesc(ukp1);figure(gcf);
title('2-D Laplace''s equation')
xlabel('B')
ylabel('H')
%--------------------------FUNCTION---------------------------------------------------------% function x = thomas(a,b,c,d) N = length(b); c(1) = c(1)/b(1); d(1) = d(1)/b(1); for i = 2 : N - 1 c(i) = c(i)/(b(i) - a(i)*c(i-1)); d(i) = (d(i) - a(i)*d(i-1))/(b(i) - a(i)*c(i-1)); end d(N) = (d(N) - a(N)*d(N-1))/(b(N) - a(N)*c(N-1)); x(N) = d(N); for i = N - 1 : -1 : 1 x(i) = d(i) - c(i)*x(i+1); end
