vectorized vs for loop data allocation, they should come out to the same matrix but they do not.

1 次查看（过去 30 天）
Sam Blake
Sam Blake 2024-3-14
评论： Voss 2024-3-14
Hello, I am trying to do some numerical methods and while setting up my A matrix, i found something I couldnt figure out. I have the exact same equations leading to the same diagonals (at least they should be) in a vectorized and a for loop method. When I run this code, the two matrices, which should be the same just having different methods of data allocation, differ. I am not sure why as they should have the same values with the exception of the top two and bottom two rows as those are to be used for boundary conditions.Also, the vectorized method should be a sparse matrix due to spdiags fcn. here is the code, the first half is mostly set up, the section " Define Coefficients for Beam Bending" is where the matrix set up begins. so just to be explicitly clear, the matrix 'Mat' and 't' should be the same with the exception of the top and bottom two rows (zeroed out later for boundary conditions) but they differ. Thank you in advance to anyone willing to take their time to look it over.
%% discretize domains
n = 11 ; % Amount of Nodes
% Computational Domain
xi = linspace(0,L,n) ;
dxi = xi(2) - xi(1) ;
% Physical Domain
den = cos((7*pi)/8) - cos(pi/8) ;
term = pi*(((3*xi)/4)+(1/8));
x = (cos(term)-cos(pi/8)) / den ;
% Plotting Domains
% figure(1) ; plot(x,x.*0,'g*-') ; title('Physical Domain)') ;
% xlabel('Length Along Beam (m)') ;
% figure(2) ;plot(xi,xi.*0,'r*-') ; title('Computational Domain') ;
% xlabel('Length Along Beam (m)') ;
%% Define Mapping Function Derivatives
J = ((-3*pi)/4)*(sin(term) / den) ; % First Derivative
J2 = ((-9*pi^2)/16)*(cos(term) / den) ; % Second Derivative
J3 = ((27*pi^3)/64)*(sin(term) / den) ; % Third Derivative
J4 = ((81*pi^4)/256)*(cos(term) / den) ; % Fourth Derivative
%% Define Coefficients for Beam Bending Equation
A = J.^-4 ; % 4th deriv term
B = (-6.*(J.^-5)).*J2 ; % 3rd deriv term
C = ((15.*(J.^-6)).*(J2.^2)) - ((4.*(J.^-5)).*J3) ; % 2nd deriv term
% 1st Derivative Term
D = ((-15.*(J.^-7)).*(J2.^3))+((10.*(J.^-6)).*J2.*J3)-((J.^-5).*J4) ;
%% Set Up Matrix
diags = [ (A.*(1/dxi^4))-(B.*(1/(2*dxi^3)))
(A.*(-4/dxi^4))+(B.*(1/dxi^3))+(C.*(1/dxi^2))-(D*(1/(2*dxi)))
((A).*(6/dxi^4))-(C.*(2/dxi^2))
(A.*(-4/dxi^4))-(B.*(1/dxi^3))+(C.*(1/dxi^2))+(D.*(1/(2*dxi)))
(A.*(1/dxi^4))+(B.*(1/(2*dxi^3))) ]' ;
Mat = spdiags(diags, -2:2, n,n) ;
t = zeros(n,n) ;
for i = 3:n-2
t(i,i-2) = (A(i).*(1/dxi^4))-(B(i).*(1/(2*dxi^3))) ;
t(i,i-1) = (A(i).*(-4/dxi^4))+(B(i).*(1/dxi^3))+(C(i).*(1/dxi^2))-(D(i)*(1/(2*dxi))) ;
t(i,i) = ((A(i)).*(6/dxi^4))-(C(i).*(2/dxi^2)) ;
t(i,i+1) = (A(i).*(-4/dxi^4))-(B(i).*(1/dxi^3))+(C(i).*(1/dxi^2))+(D(i).*(1/(2*dxi))) ;
t(i,i+2) = (A(i).*(1/dxi^4))+(B(i).*(1/(2*dxi^3))) ;
end
  1 个评论
Sam Blake
Sam Blake 2024-3-14
Also I just want to say, dont get caught up in the equation. the two equations are copy and pastes of eachother so the issue is very likely not that. I think it might have to do with something in the spdiags fcn but I really cant be sure.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Voss
Voss 2024-3-14
Ran in:
Inside the loop that where t is constructed, the indexing on the right side of the assignments is not the same as what spdiags does.
Consider the first iteration of that loop: when i = 3, in the first assignment you are assigning to t(3,1), based on A(3) and B(3). But spdiags assigns to Mat(3,1) the value diags(1,1), which is calculated from A(1) and B(1). So, to be the same as what spdiags does, the assignment to t(3,1) should be based on A(1) and B(1), hence A(i-2) and B(i-2).
Similarly the other indices on the right side of those assignments (except the third one for the main diagonal), need to be adjusted. See below:
L = 1;
%% discretize domains
n = 11 ; % Amount of Nodes
% Computational Domain
xi = linspace(0,L,n) ;
dxi = xi(2) - xi(1) ;
% Physical Domain
den = cos((7*pi)/8) - cos(pi/8) ;
term = pi*(((3*xi)/4)+(1/8));
x = (cos(term)-cos(pi/8)) / den ;
% Plotting Domains
% figure(1) ; plot(x,x.*0,'g*-') ; title('Physical Domain)') ;
% xlabel('Length Along Beam (m)') ;
% figure(2) ;plot(xi,xi.*0,'r*-') ; title('Computational Domain') ;
% xlabel('Length Along Beam (m)') ;
%% Define Mapping Function Derivatives
J = ((-3*pi)/4)*(sin(term) / den) ; % First Derivative
J2 = ((-9*pi^2)/16)*(cos(term) / den) ; % Second Derivative
J3 = ((27*pi^3)/64)*(sin(term) / den) ; % Third Derivative
J4 = ((81*pi^4)/256)*(cos(term) / den) ; % Fourth Derivative
%% Define Coefficients for Beam Bending Equation
A = J.^-4 ; % 4th deriv term
B = (-6.*(J.^-5)).*J2 ; % 3rd deriv term
C = ((15.*(J.^-6)).*(J2.^2)) - ((4.*(J.^-5)).*J3) ; % 2nd deriv term
% 1st Derivative Term
D = ((-15.*(J.^-7)).*(J2.^3))+((10.*(J.^-6)).*J2.*J3)-((J.^-5).*J4) ;
%% Set Up Matrix
diags = [ (A.*(1/dxi^4))-(B.*(1/(2*dxi^3)))
(A.*(-4/dxi^4))+(B.*(1/dxi^3))+(C.*(1/dxi^2))-(D*(1/(2*dxi)))
((A).*(6/dxi^4))-(C.*(2/dxi^2))
(A.*(-4/dxi^4))-(B.*(1/dxi^3))+(C.*(1/dxi^2))+(D.*(1/(2*dxi)))
(A.*(1/dxi^4))+(B.*(1/(2*dxi^3))) ]' ;
Mat = spdiags(diags, -2:2, n,n) ;
t = zeros(n,n) ;
for i = 3:n-2
t(i,i-2) = (A(i-2).*(1/dxi^4))-(B(i-2).*(1/(2*dxi^3))) ;
t(i,i-1) = (A(i-1).*(-4/dxi^4))+(B(i-1).*(1/dxi^3))+(C(i-1).*(1/dxi^2))-(D(i-1)*(1/(2*dxi))) ;
t(i,i) = ((A(i)).*(6/dxi^4))-(C(i).*(2/dxi^2)) ;
t(i,i+1) = (A(i+1).*(-4/dxi^4))-(B(i+1).*(1/dxi^3))+(C(i+1).*(1/dxi^2))+(D(i+1).*(1/(2*dxi))) ;
t(i,i+2) = (A(i+2).*(1/dxi^4))+(B(i+2).*(1/(2*dxi^3))) ;
end
Now they are the same, except the first two rows and the last two rows:
isequal(Mat(3:end-2,:),t(3:end-2,:))
ans = logical
1
  4 个评论
显示 2更早的评论
Sam Blake
Sam Blake 2024-3-14
oh okay I see now. I needed to shorten the two left and right terms to for the terms cut off by the top and bottom of the boundary conditions! That makes a lot of sense, thank you so much!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Oil, Gas & Petrochemical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by