how to add matrices to global matrix

Question

I have a problem with matrices in this code clc clear all dbstop if error % Initial data: TubeLength =2.38e-9; %%TubeLength = input ('TubeLength'); %m TubeDiameter=1.62e-9; %%%TubeDiameter = input ('TubeDiameter'); %m r=TubeDiameter/2; A=pi*r^2; D = 3; %degrees of freedom at each node presc_disp=1e-9; %m L=0.141e-9; %m Length of a C-C bond numIterations=1; %number of iterations numNaN=0; % %%%%% the ref . %%Force constants from Tserpes et al. (2005): % Kr is EA=6.52e-7 %N·nm^-1 % ktheta is EI=8.76e-10 %N·nm·rad^-2 % ktao is GJ=2.78e-10 %N·nm·rad^-2 EA=L*652 %N EI=L*8.76e-19 %N·m^2·rad^-2 GJ=L*2.78e-19; %N·m^2·rad^-2 % Transformation to nondimensional space: % [F]=EI/L^2; [L]=L. Then, multiply the resulting forces with % EI/L^2 and displacements with L to get forces and % displacements in nm. Lreal=L; EIreal=EI; presc_disp_real=presc_disp; EA=EA/EI*L^2; GJ=GJ/EI; presc_disp=presc_disp/L; L=1; EI=1; % Stiffness matrix of a C-C bond (adapted to 3 degree of freedom): KBond=[EA/L 0 0 -EA/L 0 0 0 12*EI/L^3 0 0 -12*EI/L^3 0 ; 0 0 12*EI/L^3 0 0 -12*EI/L^3; -EA/L 0 0 EA/L 0 0 ; 0 -12*EI/L^3 0 0 12*EI/L^3 0 ; 0 0 -12*EI/L^3 0 0 12*EI/L^3] % Calculations: [elementss1]=xlsread('elementss1.xlsx'); [nodess1]=xlsread('nodess1.xlsx'); %%nnodess = input ('number of nodes'); %%nElementss = input ('number of elements'); nnodess1= 200; %%% number of nodes nElementss2 = 200; %%%% number of elements nodesDef=zeros(1,nnodess1); k=zeros(nnodess1*D,nnodess1*D); % Defects: defectsPercentage=1; %per cent numDef=round((nnodess1*defectsPercentage)/100); for i=1:numDef def=round((nodess1-1)*rand)+1; %The addition of 1 and -1 is to avoid zero as the resulting random number. elements(D,:)=[]; nnodess1=nnodess1-1; end zeroMatrix = [0 0 0; 0 0 0; 0 0 0]; vi = [1 0 0]; vj = [0 1 0]; vk = [0 0 1]; for i=1:nnodess1-1 n = nodess1(2,i); m = nodess1(3,i) localZ = [1-2, 3-2, 3-4]; localZ = localZ/sqrt(localZ*localZ'); % Find local y if not(localZ(3)==0) y3 = (-localZ(1)-2*localZ(2))/localZ(3); localY = [1 2 y3]; elseif not(localZ(2)==0) y2 = (-localZ(1)-2*localZ(3))/localZ(2); localY = [1 y2 2]; elseif not(localZ(1)==0) y1 = (-localZ(2)-2*localZ(3))/localZ(1); localY = [y1 1 2]; end localY = localY/sqrt(localY*localY'); % Find local z localX = cross(localY, localZ); localX = localX/sqrt(localX*localX'); % Values for rotation matrix cosZI = localZ*vi'; cosZJ = localZ*vj'; cosZK = localZ*vk'; cosYI = localY*vi'; cosYJ = localY*vj'; cosYK = localY*vk'; cosXI = localX*vi'; cosXJ = localX*vj'; cosXK = localX*vk'; rotationSmall = [ cosZI cosZJ cosZK; cosYI cosYJ cosYK; cosXI cosXJ cosXK]; RR = [rotationSmall zeroMatrix;zeroMatrix rotationSmall]; KElement = abs(RR'*KBond*RR) n = nodess1(2,i); m = nodess1(3,i); k((n-1)*D+1:n*D,(n-1)*D+1:n*D)=k((n-1)*D+1:n*D, (n-1)*D+1:n*D)+KElement(1:D, 1:D); k((m-1)*D+1:m*D,(m-1)*D+1:m*D)=k((m-1)*D+1:m*D, (m-1)*D+1:m*D)+KElement(D+1:2*D, D+1:2*D); k((n-1)*D+1:n*D,(m-1)*D+1:m*D)=k((n-1)*D+1:n*D, (m-1)*D+1:m*D)+KElement(1:D, D+1:D*2); k((m-1)*D+1:m*D,(n-1)*D+1:n*D)=k((m-1)*D+1:m*D, (n-1)*D+1:n*D)+KElement(D+1:D*2, 1:D); end head = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494; 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ; 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ; 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 ]; fload = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 ]*presc_disp; dof=3;nNodes=100; save k_real.mat k; load=zeros(1,dof*nNodes); for i=1:size(head(1,:),2) for j=1:dof if head(j+1,i)~=0 % displacement prescribed for m=1:dof*nNodes load(m)=load(m)-k(m,(((head(1,i)-1)*dof+j)*fload(j+1,i))); end k((head(1,i))*(dof+j),:)=zeros(1,dof*2*nNodes); k(:,(head(1,i)-1)*dof+j)=zeros(dof*2*nNodes,1); k((head(1,i)-1)*dof+j,(head(1,i)-1)*dof+j)=1.0; load(1,(head(1,i)-1)*dof+j)=fload(j+1,i); if fload(j+1,i) ~= 0 kdof=j; else % force prescribed load(1,(head(1,i)-1)*dof+j)=fload(j+1,i); if fload(j+1,i) ~= 0 kdof=j; end end end end end clear RR KEelement; clear elements; disps=loads/k; clear k; loads k_real.mat forces=k*disps' clear k tip=[476 477 478 479 482 483 488 489 490 491 492 493 494]; % 'tip'contains the node numbers at one tip avdisp=0; for i=1:size(tip,2) tip_disps(i)=disps((tip(i)-1)*dof+kdof); avdisp=avdisp+disps((tip(i)-1)*dof+kdof); end fixed=[1 4 5 16 17 34 35 51 52 72 73 86 87]; % 'fixed'contains the node numbers at the other tip totforce=0; for i=1:size(tip,2) totforce=totforce+forces((tip(i)-1)*dof+kdof); end avdisp=avdisp/size(tip,2); avdisp=avdisp*Lreal; % in meters tip_disps=tip_disps*Lreal; % in meters totforce = totforce * EIreal/Lreal^2; % in Newtons YY=totforce/(TubeDiameter*avdisp/TubeLength*pi) % N/m=Pa*m Young(it)=YY; it=it+1; the error was Index in position 2 is invalid. Array indices must be positive integers or logical values. Error in Untitled1111111 (line 129) load(m)=load(m)-k(m,(((head(1,i)-1)*dof+j)*fload(j+1,i)));

John Michell · Accepted Answer

To add local element matrices into a global matrix in MATLAB, map local degrees of freedom (DOFs) to global index positions using submatrix indexing: K_global(dofs, dofs) = K_global(dofs, dofs) + k_local. For large-scale numerical simulations (such as Finite Element Analysis or Power Flow Y-bus matrices), assemble using index triplets with sparse(I, J, V, N, N) to prevent memory reallocation bottlenecks and achieve optimal computational speed.
Method 1: Direct Submatrix Indexing (Small to Medium Systems)
Preallocate the global matrix with zeros and add each local element matrix at its corresponding global index positions inside a loop:
% =========================================================================
% Direct Assembly of Global Stiffness Matrix (1D Bar/Spring Elements)
% =========================================================================
clc;
clear;

num_nodes = 4; % Total nodes (Global Matrix Size: 4x4)
num_elements = 3;

% Preallocate dense global matrix with zeros
K_global = zeros(num_nodes, num_nodes);

% Element Connectivity: [Node_Start, Node_End]
connectivity = [
    1, 2;
    2, 3;
    3, 4
];

% Local 2x2 stiffness matrix for each element: k * [1, -1; -1, 1]
k_elem = [ 100, -100;
          -100,  100];

% Loop over each element and add to global matrix
for elem = 1:num_elements
    % Extract global DOFs for the current element
    dofs = connectivity(elem, :);
    
    % In-place addition of local matrix into global coordinates
    K_global(dofs, dofs) = K_global(dofs, dofs) + k_elem;
end

disp('Assembled Global Matrix:');
disp(K_global);
Method 2: High-Performance Vectorized Sparse Assembly (Best for FEA & Large Grids)
In large models (10,000+ nodes), repeatedly updating K_global(dofs, dofs) in a loop triggers severe memory reallocations. The industry-standard MATLAB method collects index triplets $(I, J, V)$ and builds the global sparse matrix in a single call:
% =========================================================================
% High-Performance Sparse Assembly using sparse(I, J, V, N, N)
% =========================================================================
num_elements = 10000;
num_nodes = num_elements + 1;

% Preallocate triplet vectors (each 2x2 element contributes 4 entries)
entries_per_elem = 4;
total_entries = num_elements * entries_per_elem;

I_idx = zeros(total_entries, 1);
J_idx = zeros(total_entries, 1);
V_val = zeros(total_entries, 1);

ptr = 1;
for elem = 1:num_elements
    dofs = [elem, elem + 1]; % Global node indices
    k_local = [ 200, -200;
               -200,  200];
    
    % Generate meshgrid of local-to-global indices
    [J_grid, I_grid] = meshgrid(dofs, dofs);
    
    % Store in triplet vectors
    I_idx(ptr:ptr+3) = I_grid(:);
    J_idx(ptr:ptr+3) = J_grid(:);
    V_val(ptr:ptr+3) = k_local(:);
    
    ptr = ptr + 4;
end

% Build the entire global matrix instantly (automatically sums overlapping entries)
K_sparse_global = sparse(I_idx, J_idx, V_val, num_nodes, num_nodes);

fprintf('Global Sparse Matrix Assembled: %d x %d with %d non-zero entries.\n', ...
    size(K_sparse_global, 1), size(K_sparse_global, 2), nnz(K_sparse_global));
Method 3: Block Diagonal Assembly (Independent Subsystems)
If you want to combine multiple independent submatrices into a larger block-diagonal matrix without overlap, use blkdiag():
A = [1, 2; 3, 4];
B = [5, 6; 7, 8];
C = [9];

% Combine submatrices along diagonal
Global_Block = blkdiag(A, B, C);
disp(Global_Block);
Assembly Methods Comparison



Method
Best Suited For
Memory / Speed Performance




K(idx, idx) = K(idx, idx) + k
Small models (< 500 nodes), educational exercises.
Simple syntax, but slow for large systems due to quadratic memory operations.


sparse(I, J, V, N, N)
Large FEA simulations, structural grids, CFD solvers.
Fastest; automatically sums identical $(i,j)$ coordinate pairs with minimal RAM.


blkdiag(A, B, C, ...)
State-space modeling, decoupled dynamical systems.
Zero-overlap diagonal block assembly.



Key Best Practices

Always Preallocate: Never expand the global matrix dynamically inside a loop with unallocated matrices (e.g., K = []). Preallocate with zeros(N, N) or use triplet arrays.
Use Sparse Matrices for FEA: Real-world finite element matrices are typically 95%+ zero entries. Using sparse() reduces memory usage from gigabytes to megabytes.
Verify Symmetry: For conservative physical systems, verify stiffness symmetry using issymmetric(K_global).

how to add matrices to global matrix

Method 1: Direct Submatrix Indexing (Small to Medium Systems)

Method 2: High-Performance Vectorized Sparse Assembly (Best for FEA & Large Grids)

Method 3: Block Diagonal Assembly (Independent Subsystems)

Assembly Methods Comparison

Key Best Practices

Need a Custom Version or Complete Simulation for This Problem?

How to I use get param to get a value from an output of a bus?

How do I uninstall and reinstall MathWorks Service Host?

Why am I receiving error messages about singularities in my Simulink model?

Simulink turn off data logging

MatLab App Design - how to use default value for Edit Field function?

How do I extract data from MATLAB figures?

IEEE 13-BUS Hybrid AC & DC Power Distribution System using MATLAB Simulink

Adaptive Path Planning and Collision Avoidance for Autonomous Vehicles on Unstructured Indian Roads

Lithium-Ion Battery Electro-Thermal Modeling and Cooling Analysis in MATLAB

Get a Free Consultation or a Sample Assignment Review!

Method	Best Suited For	Memory / Speed Performance
`K(idx, idx) = K(idx, idx) + k`	Small models (< 500 nodes), educational exercises.	Simple syntax, but slow for large systems due to quadratic memory operations.
`sparse(I, J, V, N, N)`	Large FEA simulations, structural grids, CFD solvers.	Fastest; automatically sums identical \((i,j)\) coordinate pairs with minimal RAM.
`blkdiag(A, B, C, ...)`	State-space modeling, decoupled dynamical systems.	Zero-overlap diagonal block assembly.