could you help me in this code
1 view (last 30 days)
Show older comments
% 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 \
0 Comments
Answers (1)
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)
size(part2)
rhs(1:n:n*(n-2)+1) = part1 + part2
rhs(n:n:n^2-n+1) = rhs(n:n:n^2-n+1) - U(:,n) / h^2;
2 Comments
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.
See Also
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!