I'm plotting stress trajectories for simply supported rectangular beam under UDL. Attaching the code and it's output. Help me sort out the error & get correct trajectories.

data = readtable('data.xlsx');

Warning: Column headers from the file were modified to make them valid MATLAB identifiers before creating variable names for the table. The original column headers are saved in the VariableDescriptions property.
Set 'VariableNamingRule' to 'preserve' to use the original column headers as table variable names.

%here data is name of excel file from which i'm importing data of x,y grid & principal angles
x = data.x;
y = data.y;
theta1 = data.('Theta1');
theta2 = data.('Theta2');

%% ===== Build structured grid =====
xu = x;    % unique x-values
yu = y;    % unique y-values

Nx = length(xu);
Ny = length(yu);

% Create 2D mesh
[X, Y] = meshgrid(xu, yu);
size(theta1)

ans = 1×2

77 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Nx*Ny

ans = 5929

%% ===== Reshape angles into matrices =====
Theta1 = reshape(theta1, Ny, Nx);   % IMPORTANT: Ny rows, Nx cols

Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.

Theta2 = reshape(theta2, Ny, Nx);

%% ===== Compute direction vectors =====
U1 = cos(Theta1);  V1 = sin(Theta1);
U2 = cos(Theta2);  V2 = sin(Theta2);

%% ===== Create seeding points =====
[startX, startY] = meshgrid( linspace(min(xu),max(xu),8), ...
                              linspace(min(yu),max(yu),8) );

startX = startX(:);
startY = startY(:);

%% ===== Plot Trajectories =====
figure; hold on;

% Principal direction 1
streamline(X, Y, U1, V1, startX, startY);

% Principal direction 2
streamline(X, Y, U2, V2, startX, startY);

quiver(X, Y, U1, V1, 0.5, 'k');   % direction field (optional)

axis equal tight;
xlabel('X'); ylabel('Y');
title('Stress Trajectories For Simply Supported Steel Beam');
grid on;

[Syeda Hajra]

attaching data file in comment

[Torsten]

theta1 is 77x1, Nx*Ny equals 5929. No way to reshape because 77*1 doesn't equal 5929 (see above).

[Syeda Hajra]

i have attached the data sheet in comment below if you can help me how to plot stress trajectories using that data

Answer 1

[Star Strider]

According to the streamline documentation, 'startX' and 'startY' are supposed to be matrices rather than vectors.

Beyond that, I cannot run your code (no data).

[Syeda Hajra]

here is the data sheet i used in above code. x and y are coordinated of grid points and theta1 and theta2 are principal angles. need to plot stress trajectories using this data

[Star Strider]

I am not certain how to work with these data, so I need your help in interpreting them.

The and vectors can be made into matrices (with a bit of interpolation to make the shorter vectors have the same lengths as the longer vectors, although it would be easier to simply ignore the shorter vectors), however I am not certain what data should be included in the matrices.

How do 'x' and 'y' enter into this?

I began by exploring the data graphically --

data = readtable('data.xlsx', VariableNamingRule='preserve');

%here data is name of excel file from which i'm importing data of x,y grid & principal angles
x = data.x;
y = data.y;
theta1 = data.('Theta 1');
theta2 = data.('Theta 2');

disp(x)

     0
     0
     0
     0
     0
     0
     0
     2
     2
     2
     2
     2
     2
     2
     4
     4
     4
     4
     4
     4
     4
     6
     6
     6
     6
     6
     6
     6
     8
     8
     8
     8
     8
     8
     8
    10
    10
    10
    10
    10
    10
    10
    12
    12
    12
    12
    12
    12
    12
    14
    14
    14
    14
    14
    14
    14
    16
    16
    16
    16
    16
    16
    16
    18
    18
    18
    18
    18
    18
    18
    20
    20
    20
    20
    20
    20
    20

disp(y)

    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000
    2.0000
    1.3300
    0.6600
         0
   -0.6600
   -1.3300
   -2.0000

Lvp1 = islocalmax(theta1, FlatSelection='last');
pks1 = theta1(Lvp1);
plocs1 = find(Lvp1);
Lvv1 = islocalmin(theta1, FlatSelection='first');
vys1 = theta1(Lvv1);
vlocs1 = find(Lvv1);

Lvp2 = islocalmax(theta2, FlatSelection='last');
pks2 = theta2(Lvp2);
plocs2 = find(Lvp2);
Lvv2 = islocalmin(theta2, FlatSelection='first');
vys2 = theta2(Lvv2);
vlocs2 = find(Lvv2);

xvi = 1:size(data,1);

figure
plot(xvi, theta1, DisplayName='\Theta_1')
hold on
hp2 = plot(xvi, theta2, DisplayName='\Theta_2');
plot(xvi(plocs1), pks1, '^b', DisplayName='\Theta_1 Peaks')
plot(xvi(vlocs1), vys1, 'vb', DisplayName='\Theta_1 Valleys')
plot(xvi(plocs2), pks2, '^', DisplayName='\Theta_2 Peaks', Color=hp2.Color)
plot(xvi(vlocs2), vys2, 'v', DisplayName='\Theta_2 Valleys', Color=hp2.Color)
hold off
grid
legend(Location='eastoutside')
xlabel("Index")
ylabel("\Theta (radians)")
title("See What The Data Look Like ...")

figure
plot(x, y)
grid
xlabel('x')
ylabel('y')

Both Θ vectors appear to be nearly symmetrical about element 39.

It would be easier to make this a 'one-off', however I would like to make it robust to other similar data sets that you might want to analyse.

I am working with these in MATLAB Online, so I have done more than I have illustrated here.

EDIT -- (02 Dec 2025 at 01:46)

I managed to get this far, although I had to edit your data to eliminate the short vectors in 'theta1' and 'theta2'. I have not used streamline in several years, and then with some data I understood. (I do not understand your data.)

This keeps throwing a 'Sample points must be unique.' error that I understand, however I do not understand it in the context of your data. I was hoping to make this work completely, however the error message is sufficiently uninformative that I have no idea what the problem actually is. The documentation appears to be lacking in describing any details of the argument requirements.

I am not enthusiastic about delving into the streamline code to see what is being interpolated. The row-wise and column-wise elements of 'theta1mtx' and 'theta2mtx' are not unique. That may be the problem.

I will continue to help you with this, however I will need your help in understanding how to apply streamline to your data.

This seems to work to create the matrices --

data = readtable('data.xlsx', VariableNamingRule='preserve');

%here data is name of excel file from which i'm importing data of x,y grid & principal angles
x = data.x;
y = data.y;
theta1 = data.('Theta 1');
theta2 = data.('Theta 2');

Lvp1 = islocalmax(theta1, FlatSelection='last');
pks1 = theta1(Lvp1);
plocs1 = find(Lvp1);
Lvv1 = islocalmin(theta1, FlatSelection='first');
vys1 = theta1(Lvv1);
vlocs1 = find(Lvv1);

Lvp2 = islocalmax(theta2, FlatSelection='last');
pks2 = theta2(Lvp2);
plocs2 = find(Lvp2);
Lvv2 = islocalmin(theta2, FlatSelection='first');
vys2 = theta2(Lvv2);
vlocs2 = find(Lvv2);

xvi = 1:size(data,1);

figure
plot(xvi, theta1, DisplayName='\Theta_1')
hold on
hp2 = plot(xvi, theta2, DisplayName='\Theta_2');
plot(xvi(plocs1), pks1, '^b', DisplayName='\Theta_1 Peaks')
plot(xvi(vlocs1), vys1, 'vb', DisplayName='\Theta_1 Valleys')
plot(xvi(plocs2), pks2, '^', DisplayName='\Theta_2 Peaks', Color=hp2.Color)
plot(xvi(vlocs2), vys2, 'v', DisplayName='\Theta_2 Valleys', Color=hp2.Color)
hold off
grid
% legend(Location='best')
xlabel("Index")
ylabel("\Theta (radians)")
title("See What The Data Look Like ...")

% text(xvi(vlocs1), vys1, compose('%d', vlocs1))
% text(xvi(plocs1), pks1, compose('%d', plocs1))

figure
hold on
for k = 2:numel(vlocs1)/2
    % k
    idxrng = plocs1(k):vlocs1(k);
    plot(xvi(idxrng), theta1(idxrng));
    theta1mtx(:,k-1) = theta1(idxrng);
end
% QQQ = idxrng(end)
for k = k:numel(plocs1)-2
    % k
    % idxends = [vlocs1(k) plocs1(k+1)]
    idxrng = vlocs1(k):plocs1(k+1);
    plot(xvi(idxrng), theta1(idxrng));
    theta1mtx(:,k-1) = theta1(idxrng);
end
hold off
grid
xlabel("Index")
ylabel("\Theta")
title("\Theta_1")

●

theta1mtx

theta1mtx = 7×7

0.5880 0.3958 0.2375 -0.7854 -0.7854 -0.7854 -0.7854 0.3203 0.1524 0.0792 0 -0.1107 -0.2375 -0.3958 0 0 0 0 -0.0349 -0.0792 -0.1524 0 0 0 0 0 0 0 -0.1524 -0.0792 -0.0349 0 0 0 0 -0.3958 -0.2375 -0.1107 0.0349 0.0792 0.1524 0.3203 -0.7854 -0.7854 -0.7854 0.1107 0.2375 0.3958 0.5880

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
hold on
for k = 2:numel(vlocs2)/2
    % k
    idxrng = plocs2(k):vlocs2(k);
    plot(xvi(idxrng), theta2(idxrng));
    theta2mtx(:,k-1) = theta2(idxrng);
end
% QQQ = idxrng(end)
for k = k:numel(plocs2)-2
    % k
    % idxends = [vlocs2(k) plocs2(k+1)]
    idxrng = vlocs2(k):plocs2(k+1);
    plot(xvi(idxrng), theta2(idxrng));
    theta2mtx(:,k-1) = theta2(idxrng);
end
hold off
grid
xlabel("Index")
ylabel("\Theta")
title("\Theta_2")

theta2mtx

theta2mtx = 7×7

2.1588 1.9666 1.8083 0.7854 0.7854 0.7854 0.7854 1.8911 1.7232 1.6500 1.5708 1.4601 1.3333 1.1750 1.5708 1.5708 1.5708 1.5708 1.5359 1.4916 1.4184 1.5708 1.5708 1.5708 1.5708 1.5708 1.5708 1.5708 1.4184 1.4916 1.5359 1.5708 1.5708 1.5708 1.5708 1.1750 1.3333 1.4601 1.6057 1.6500 1.7232 1.8911 0.7854 0.7854 0.7854 1.6815 1.8083 1.9666 2.1588

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(x, y, 'p-')
grid
xlabel('x')
ylabel('y')
axis('padded')

% Sz_x = size(x)
% Sz_y = size(y)

xmtx = reshape(x, 7, []).';
ymtx = reshape(y, 7, []).';

% Sz_xm = size(xmtx)
% Sz_ym = size(ymtx)

xmtx = xmtx(1:7,:);  
ymtx = ymtx(1:7,:);

Sz_xm = size(xmtx)              % Clip 'xmtx' To Equal 'theta1'

Sz_xm = 1×2

7 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Sz_ym = size(ymtx)              % Clip 'ymtx' To Equal 'theta2'

Sz_ym = 1×2

7 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Sz_th1m = size(theta1mtx)
% Sz_th2m = size(theta2mtx)



%% ===== Build structured grid =====
xu = unique(x);    % unique x-values
yu = unique(y);    % unique y-values

Nx = length(xu);
Ny = length(yu);

% Create 2D mesh
% [X, Y] = meshgrid(xu, yu);
% size(theta1)
% Nx*Ny
%% ===== Reshape angles into matrices =====
% Theta1 = reshape(theta1, Ny, Nx);   % IMPORTANT: Ny rows, Nx cols
% Theta2 = reshape(theta2, Ny, Nx);

Theta1 = theta1mtx;
Theta2 = theta2mtx;

%% ===== Compute direction vectors =====
U1 = cos(Theta1);  V1 = sin(Theta1);
U2 = cos(Theta2);  V2 = sin(Theta2);

%% ===== Create seeding points =====
% [startX, startY] = meshgrid( linspace(min(xu),max(xu),8), ...
%                               linspace(min(yu),max(yu),8) );

% startX = startX(:);
% startY = startY(:);

startX = xmtx(:,1);
startY = ymtx(1,:);

X = xmtx;
Y = ymtx;

%% ===== Plot Trajectories =====
figure; hold on;

% Principal direction 1
streamline(X, Y, U1, V1, startX, startY);

Error using matlab.internal.math.interp1
Sample points must be unique.

Error in interp1 (line 188)
VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in stream2 (line 62)
sxi=interp1(xx(:),1:szu(2),sx(k));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in streamline (line 21)
verts = stream2(x,y,u,v,sx,sy,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% Principal direction 2
streamline(X, Y, U2, V2, startX, startY);

quiver(X, Y, U1, V1, 0.5, 'k');   % direction field (optional)

axis equal tight;
xlabel('X'); ylabel('Y');
title('Stress Trajectories For Simply Supported Steel Beam');
grid on;

.