Can someone please help me with coding a a general expression of n-stencils (where n
Please ignore the code above.
function D = myfinD()
x = chebpts(6);
%x = linspace(-1,1,6)';
stencils = 3;
x_ = chebpts(stencils);
%x_ = linspace(-1,1,stencils)';
n = numel(x);
d2L = zeros(n,stencils); % Defining an empty matrix
for i = 1:stencils
dL = zeros(n,1);
for k = 1:n
if (x_(i)~=x(k))
for l = 1:n
if (x_(i)~=x(l) && l~=k)
Lprod = ones(n,1);
for j = 1:n
if (x_(i)~=x(j) && j~= k && j~=l)
Lprod = Lprod.*(x-x(j))/(x_(i)-x(j));
end
end
dL = dL + Lprod/((x_(i)-x(k))*(x_(i)-x(l)));
end
end
end
end
d2L(:,i) = dL;
end
full(d2L)
if (rem(stencils,2)~=0)
if (stencils ~= 1)
m = (stencils-1)/2;
D = spdiags(d2L, -m:m, n, n);
else
D = spdiags(d2L, 0, n, n);
end
full(D)
else
m = stencils/2;
D = spdiags([d2L(:,1:m),zeros(n,1),d2L(:,m+1:end)], -m:m, n, n);
full(D)
end
end
Hello, check something called Kronecker product. I upload my code for 5-point stencil FDM (C is a matrix where every point defines conductivity of the domain at that point):
function U = FDM2D(C)
[nx,ny] = size(C);
dx = 1/(nx-1);
dy = 1/(ny-1);
a = ones(nx,1);
b = -2*ones(nx,1);
c = ones(nx,1);
sq = [a,b,c];
sq = spdiags(sq,-1:1,nx,nx);
sq([nx+1 nx*ny-nx]) = [2 2];
Dxx = kron(sq,eye(ny))/(dx)^2;
Dyy = kron(eye(ny),sq)/(dy)^2;
A = Dxx + Dyy;
A = A.*repmat(U2vek(C)',nx*ny,1);
A(1,:) = [1,zeros(1,nx*ny-1)];
A(end,:) = [zeros(1,nx*ny-1),1];
b = [-100;zeros(nx*ny-2,1);100];
u = A\b;
U = u2mat(u,nx,ny);
% contourf(flip(U),20); axis square; colorbar; title('Pressure')
end
% U2vek(C) function convert conductivity C matrix to a vector form c
function u = U2vek(U)
[nx,ny] = size(U);
u = zeros(nx*ny,1);
for ii=1:ny
c = 1+(ii-1)*nx; % First element of ith U line
l = U(ii,:)'; % Take the ith U line
u(end-ii*nx+1:end-(ii-1)*nx) = l; % Add ith line on top of U
end
end
% u2mat(u,n,m) function convert solution u vector back to its matrix form U
function U = u2mat(u,n,m)
U = [];
for i=1:m
c = 1+(i-1)*n; % First element of ith U line
l = u(c:c+(n-1)).'; % Take the ith U line
U = [l; U]; % Add ith line on top of U
end
end