could you help me in this code

1 view (last 30 days)
nana
nana on 4 May 2023
Commented: Walter Roberson on 4 May 2023
% Define grid parameters
h = 1/4; % Grid spacing
x = -1:h:1; % Grid points in x direction
y = -1:h:1; % Grid points in y direction
n = length(x); % Number of grid points
% Initialize solution matrix
U = zeros(n,n);
% Set boundary conditions
U(1,:) = 0; % u(x,y)=0 along y=1
U(n,:) = 0; % u(x,y)=0 along y=-1
U(:,1) = 0.5*U(:,2); % u(x,y)=(1/2)u, x=-1,-1<y<1
U(:,n) = -0.5*U(:,n-1); % u(x,y)=-(1/2)u, x=1,-1<y<1
% Construct matrix A
B = [-6 2 0 0 0; 1 -6 1 0 0; 0 1 -6 1 0; 0 0 1 -6 1; 0 0 0 2 -6.25];
I = eye(5);
A = kron(eye(n-2), B) + kron(diag(ones(n-3,1), -1), I) + kron(diag(ones(n-3,1), 1), I);
% Construct right-hand side vector rhs from boundary conditions
rhs = zeros(n^2, 1);
rhs(1:n) = U(1,:)' / h^2;
rhs(end-n+1:end) = -U(n,:)' / h^2;
rhs(1:n:n*(n-2)+1) = rhs(1:n:n*(n-2)+1) + U(:,1) / h^2;
rhs(n:n:n^2-n+1) = rhs(n:n:n^2-n+1) - U(:,n) / h^2;
Arrays have incompatible sizes for this operation.
% Solve the system of equations using matrix inversion
U_vec = A \ rhs;
U = reshape(U_vec, [n, n]);
% Plot the solution as a 3-D mesh
figure;
mesh(x, y, U');
xlabel('x');
ylabel('y');
zlabel('u(x,y)');
title('Solution of the Poisson equation using central differences');
Error using \

Answers (1)

Walter Roberson
Walter Roberson on 4 May 2023
You cannot add an 8 x 1 vector and a 9 x 1 vector.
% Define grid parameters
h = 1/4; % Grid spacing
x = -1:h:1; % Grid points in x direction
y = -1:h:1; % Grid points in y direction
n = length(x); % Number of grid points
% Initialize solution matrix
U = zeros(n,n);
% Set boundary conditions
U(1,:) = 0; % u(x,y)=0 along y=1
U(n,:) = 0; % u(x,y)=0 along y=-1
U(:,1) = 0.5*U(:,2); % u(x,y)=(1/2)u, x=-1,-1<y<1
U(:,n) = -0.5*U(:,n-1); % u(x,y)=-(1/2)u, x=1,-1<y<1
% Construct matrix A
B = [-6 2 0 0 0; 1 -6 1 0 0; 0 1 -6 1 0; 0 0 1 -6 1; 0 0 0 2 -6.25];
I = eye(5);
A = kron(eye(n-2), B) + kron(diag(ones(n-3,1), -1), I) + kron(diag(ones(n-3,1), 1), I);
% Construct right-hand side vector rhs from boundary conditions
rhs = zeros(n^2, 1);
rhs(1:n) = U(1,:)' / h^2;
rhs(end-n+1:end) = -U(n,:)' / h^2;
part1 = rhs(1:n:n*(n-2)+1);
part2 = U(:,1) / h^2;
size(part1)
ans = 1×2
8 1
size(part2)
ans = 1×2
9 1
rhs(1:n:n*(n-2)+1) = part1 + part2
Arrays have incompatible sizes for this operation.
rhs(n:n:n^2-n+1) = rhs(n:n:n^2-n+1) - U(:,n) / h^2;
  2 Comments
nana
nana on 4 May 2023
thank you for your calrification but how can i fix it please ?
Walter Roberson
Walter Roberson on 4 May 2023
I do not know.
I would suggest to you that it might be easier to costruct rhs as an n x n 2D array and then reshape it to vector afterwards.

Sign in to comment.

Categories

Find more on Resizing and Reshaping Matrices in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!