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

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

[Sam Blake]

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.

Answer 1

[Voss]

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

[Sam Blake]

Thank you for your quick response! okay I think I understand. Im pretty sure the matrix yielded by the for loop is actually giving me the correct matrix, So if i wanted to achieve the inverse of the example you gave, I would index the A,B,C,D coefficients from (3:n-2) in the 'diags' matrix?

[Voss]

If you want to use spdiags to generate a matrix like t (except the first two and last two rows), you can do this:

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))) ]' ;
disp(diags)

   1.0e+05 *

    4.7730   -1.4369   -7.3210    5.2308   -1.2459
    0.6251   -1.2070    0.7607   -0.1874    0.0086
    0.1814   -0.4826    0.4911   -0.2347    0.0448
    0.0816   -0.2544    0.3075   -0.1731    0.0384
    0.0495   -0.1715    0.2310   -0.1441    0.0351
    0.0378   -0.1429    0.2101   -0.1429    0.0378
    0.0351   -0.1441    0.2310   -0.1715    0.0495
    0.0384   -0.1731    0.3075   -0.2544    0.0816
    0.0448   -0.2347    0.4911   -0.4826    0.1814
    0.0086   -0.1874    0.7607   -1.2070    0.6251
   -1.2459    5.2308   -7.3210   -1.4369    4.7730

diags(:,1) = [diags(3:n,1); 0; 0];
diags(:,2) = [diags(2:n,2); 0];
diags(:,4) = [0; diags(1:n-1,4)];
diags(:,5) = [0; 0; diags(1:n-2,5)];
disp(diags)

   1.0e+05 *

    0.1814   -1.2070   -7.3210         0         0
    0.0816   -0.4826    0.7607    5.2308         0
    0.0495   -0.2544    0.4911   -0.1874   -1.2459
    0.0378   -0.1715    0.3075   -0.2347    0.0086
    0.0351   -0.1429    0.2310   -0.1731    0.0448
    0.0384   -0.1441    0.2101   -0.1441    0.0384
    0.0448   -0.1731    0.2310   -0.1429    0.0351
    0.0086   -0.2347    0.3075   -0.1715    0.0378
   -1.2459   -0.1874    0.4911   -0.2544    0.0495
         0    5.2308    0.7607   -0.4826    0.0816
         0         0   -7.3210   -1.2070    0.1814

Mat = spdiags(diags, -2:2, n, n);
disp(full(Mat))

   1.0e+05 *

   -7.3210    5.2308   -1.2459         0         0         0         0         0         0         0         0
   -1.2070    0.7607   -0.1874    0.0086         0         0         0         0         0         0         0
    0.1814   -0.4826    0.4911   -0.2347    0.0448         0         0         0         0         0         0
         0    0.0816   -0.2544    0.3075   -0.1731    0.0384         0         0         0         0         0
         0         0    0.0495   -0.1715    0.2310   -0.1441    0.0351         0         0         0         0
         0         0         0    0.0378   -0.1429    0.2101   -0.1429    0.0378         0         0         0
         0         0         0         0    0.0351   -0.1441    0.2310   -0.1715    0.0495         0         0
         0         0         0         0         0    0.0384   -0.1731    0.3075   -0.2544    0.0816         0
         0         0         0         0         0         0    0.0448   -0.2347    0.4911   -0.4826    0.1814
         0         0         0         0         0         0         0    0.0086   -0.1874    0.7607   -1.2070
         0         0         0         0         0         0         0         0   -1.2459    5.2308   -7.3210

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
disp(t)

   1.0e+04 *

         0         0         0         0         0         0         0         0         0         0         0
         0         0         0         0         0         0         0         0         0         0         0
    1.8142   -4.8265    4.9107   -2.3467    0.4483         0         0         0         0         0         0
         0    0.8162   -2.5438    3.0745   -1.7309    0.3840         0         0         0         0         0
         0         0    0.4949   -1.7151    2.3099   -1.4410    0.3513         0         0         0         0
         0         0         0    0.3782   -1.4289    2.1013   -1.4289    0.3782         0         0         0
         0         0         0         0    0.3513   -1.4410    2.3099   -1.7151    0.4949         0         0
         0         0         0         0         0    0.3840   -1.7309    3.0745   -2.5438    0.8162         0
         0         0         0         0         0         0    0.4483   -2.3467    4.9107   -4.8265    1.8142
         0         0         0         0         0         0         0         0         0         0         0
         0         0         0         0         0         0         0         0         0         0         0

isequal(Mat(3:end-2,:),t(3:end-2,:))

ans = logical
   1

[Sam Blake]

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!

[Voss]

Shorten and move up/down, yeah.

You're welcome!