Visualising 3 times 3 linear algebraic systems and their intersections

I want to interpret geometrically linear systems. Consider for example the system

With the following implementation of Rouché-Capelli Theorem one sees that the system has a unique solution

clc, clear, format compact
% matrix with coefficients of unknowns
A = [1/10, 1/3, 1/3; 1/2, 1/20, 1/2; 1/5, 1/10, 1/0.5];
% column matrix with constant terms
C = [1; 1; 1];
% Implementation of Rouché-Cappelli algorithm 
mat_augm = [A C];   % augmented matrix
% verify that if the ranks are equal (consistent system)
if (rank(mat_augm) == rank (A))  
    % verify that the rank equals the number of unknowns
    if rank(mat_augm) == size(A,2)   
        disp('Le système possède une solution unique.')
        X_sol = A\C  % solution with left division method
    else % rank is less than number of unknowns
        disp('Le système possède un nombre infini de solutions.')
        disp('La forme échelonnée du système est :')
        rref(mat_augm)  % row reduced echelon form of augmented matrix
    end
else  % Ranks are different; incosistent system
    disp('Le système est impossible.')
    disp('La forme échelonnée du système est :')
    rref(mat_augm)  % row reduced echelon form of augmented matrix
end

Le système possède une solution unique.

X_sol = 3×1

1.5385 2.3077 0.2308

I use the following instructions for visualising the three planes, their lines of intersection (of any two of them) and their point of intersection (of all of them).

clf
set(0,'defaultTextInterpreter','latex');
X_sol = [1.5385 2.3077 0.2308];
% visualisation des plans
% premier plan
x1 = 0:0.5:10; x2 = 0:0.5:3;
[X Y] = meshgrid(x1,x2);
Z = 3-Y-(3/10)*X;
H1 = surf(X,Y,Z);
set(H1,'facecolor', 'r','FaceAlpha',0.3, 'EdgeColor' ,'none') 
hold on
% deuxième plan
x1 = 0:0.5:2; x2 = 0:0.5:20;
[X Y] = meshgrid(x1,x2);
Z = 2-X-(1/10)*Y;
H2 = surf(X,Y,Z);
set(H2,'facecolor', 'b','FaceAlpha',0.3, 'EdgeColor','none') 
% troisième plan
x1 = 0:0.5:5; x2 = 0:0.5:10;
[X Y] = meshgrid(x1,x2);
Z = 0.5*(1-(1/5)*X-(1/10)*Y);
H3 = surf(X,Y,Z);
set(H3,'facecolor', 'g','FaceAlpha',0.3, 'EdgeColor' ,'none') 
% point d'intersection
plot3(X_sol(1), X_sol(2), X_sol(3), 'k.', 'MarkerSize', 30)
% lines of intersections
% line of intersection of first two planes
% calculate intersection in plane z = 0
A1 = [ 1/10 1/3 ; 1/2 1/20 ];
B1 = [ 1; 1 ];
X1 = A1\B1;
P1 = [ X1; 0 ];
% calculate intersection in plane x = 0
A2 = [ 1/3 1/3 ; 1/20 1/2 ];
B2 = [ 1; 1 ];
X2 = A2\B2;
P2 = [ 0; X2 ];
P1P2 = [P1, P2];
line(P1P2(1,:), P1P2(2,:), P1P2(3,:), 'Color',...
    'black', 'linewidth',1.5); view(3)
% line of interesection of first and third plane
% calculate intersection in plane z = 0
A1 = [ 1/10 1/3 ; 1/5 1/10 ];
B1 = [ 1; 1 ];
X1 = A1\B1;
P1 = [ X1; 0 ];
% calculate intersection in plane x = 0
A2 = [ 1/3 1/3 ; 1/10 1/0.5 ];
B2 = [ 1; 1 ];
X2 = A2\B2;
P2 = [ 0; X2 ];
P1P2 = [P1, P2];
line(P1P2(1,:), P1P2(2,:), P1P2(3,:), 'Color','black',...
    'linewidth',1.5); view(3)
% line of interesection of second and third plane
% calculate intersection in plane z = 0
A1 = [ 1/2 1/20 ; 1/5 1/10 ];
B1 = [ 1; 1 ];
X1 = A1\B1;
P1 = [ X1; 0 ];
% calculate intersection in plane y = 0
A2 = [1/2 1/2; 1/5 1/0.5];
B2 = [ 1; 1 ];
X2 = A2\B2;
P2 = [X2(1); 0; X2(2)];
P1P2 = [P1, P2];
line(P1P2(1,:), P1P2(2,:), P1P2(3,:), 'Color',...
    'black', 'linewidth',1.5); view(3)

% plot settings and labeling
axis square, grid on, grid minor, box on
xlabel('$x_1$'), ylabel('$x_2$'), zlabel('$x_3$')
hAxis = gca;
set(hAxis,'FontSize',14)
xlim([0 10]) , ylim([0 20]), zlim([0 3])
set(gca, 'Xdir', 'reverse', 'YDir', 'reverse')
xticks([0,5,10]), yticks([0, 10, 20]),  zticks([0,1,2,3])
hAxis.XTickLabels = 0:5:10;
hAxis.YTickLabels = 0:10:20;
hAxis.ZTickLabels = 0:1:3;
title({'$\frac{1}{10}x_1+\frac{1}{3}x_2+\frac{1}{3}x_3=1$ (red)',...
    '$\frac{1}{2}x_1+\frac{1}{20}x_2+\frac{1}{2}x_3=1$ (blue)',...
    '$\frac{1}{5}x_1+\frac{1}{10}x_2+\frac{1}{0.5}x_3=1$ (green)', blanks(1)})
hAxis.XRuler.FirstCrossoverValue  = 0; % X crossover with Y axis
hAxis.XRuler.SecondCrossoverValue  = 0; % X crossover with Z axis
hAxis.YRuler.FirstCrossoverValue  = 0; % Y crossover with X axis
hAxis.YRuler.SecondCrossoverValue  = 0; % Y crossover with Z axis
hAxis.ZRuler.FirstCrossoverValue  = 0; % Z crossover with X axis
hAxis.ZRuler.SecondCrossoverValue = 0; % Z crossover with Y axis
hold off

set(0,'defaultTextInterpreter','none'); 

I get the desired visualisation. But my code demands a lot of manual work in order to specify, for instance, the lines of intersection. I am wondering how I can make the processus more user-independent.

Thank you very much.

[Bjorn Gustavsson]

Isn't one of the main "possible stumbling blocks" (poor translation to English by me) that you have the origin as a part of your display. Maybe it becomes easier if you start the visualization from the unique slution, make the thre intersection-lines of suitable lengths and use them to make corner-coordinates for your planes? This obviously scraps the neat coordinate-axis-crossings of your plotted planes - but that might turn ugly in a hurry if someone selects one or two equations with large off-sets.

Answer 1

[Shushant]

Hi @Dimitrios Anagnostou,

According to my understanding, you are trying to remove the manual or redundant work of writing the same piece of code multiple times as in the case of calculating the lines of intersections. To achieve this, you can take advantage of MATLAB function. To demonstrate the same I have made a function called "calculateInterscetion".

clf
set(0,'defaultTextInterpreter','latex');
X_sol = [1.5385 2.3077 0.2308];
% visualisation des plans
% premier plan
x1 = 0:0.5:10; x2 = 0:0.5:3;
[X Y] = meshgrid(x1,x2);
Z = 3-Y-(3/10)*X;
H1 = surf(X,Y,Z);
set(H1,'facecolor', 'r','FaceAlpha',0.3, 'EdgeColor' ,'none') 
hold on
% deuxième plan
x1 = 0:0.5:2; x2 = 0:0.5:20;
[X Y] = meshgrid(x1,x2);
Z = 2-X-(1/10)*Y;
H2 = surf(X,Y,Z);
set(H2,'facecolor', 'b','FaceAlpha',0.3, 'EdgeColor','none') 
% troisième plan
x1 = 0:0.5:5; x2 = 0:0.5:10;
[X Y] = meshgrid(x1,x2);
Z = 0.5*(1-(1/5)*X-(1/10)*Y);
H3 = surf(X,Y,Z);
set(H3,'facecolor', 'g','FaceAlpha',0.3, 'EdgeColor' ,'none') 
% point d'intersection
plot3(X_sol(1), X_sol(2), X_sol(3), 'k.', 'MarkerSize', 30)
% lines of intersections
% line of intersection of first two planes
% calculate intersection in plane z = 0
A1 = [ 1/10 1/3 ; 1/2 1/20 ];
% calculate intersection in plane x = 0
A2 = [ 1/3 1/3 ; 1/20 1/2 ];
calculateIntersection(A1, A2, "z", "x");
% line of interesection of first and third plane
% calculate intersection in plane z = 0
A1 = [ 1/10 1/3 ; 1/5 1/10 ];
% calculate intersection in plane x = 0
A2 = [ 1/3 1/3 ; 1/10 1/0.5 ];
calculateIntersection(A1, A2, "z", "x");
% line of interesection of second and third plane
% calculate intersection in plane z = 0
A1 = [ 1/2 1/20 ; 1/5 1/10 ];
% calculate intersection in plane y = 0
A2 = [1/2 1/2; 1/5 1/0.5];
calculateIntersection(A1, A2, "z", "y");
% plot settings and labeling
axis square, grid on, grid minor, box on
xlabel('$x_1$'), ylabel('$x_2$'), zlabel('$x_3$')
hAxis = gca;
set(hAxis,'FontSize',14)
xlim([0 10]) , ylim([0 20]), zlim([0 3])
set(gca, 'Xdir', 'reverse', 'YDir', 'reverse')
xticks([0,5,10]), yticks([0, 10, 20]),  zticks([0,1,2,3])
hAxis.XTickLabels = 0:5:10;
hAxis.YTickLabels = 0:10:20;
hAxis.ZTickLabels = 0:1:3;
title({'$\frac{1}{10}x_1+\frac{1}{3}x_2+\frac{1}{3}x_3=1$ (red)',...
    '$\frac{1}{2}x_1+\frac{1}{20}x_2+\frac{1}{2}x_3=1$ (blue)',...
    '$\frac{1}{5}x_1+\frac{1}{10}x_2+\frac{1}{0.5}x_3=1$ (green)', blanks(1)})
hAxis.XRuler.FirstCrossoverValue  = 0; % X crossover with Y axis
hAxis.XRuler.SecondCrossoverValue  = 0; % X crossover with Z axis
hAxis.YRuler.FirstCrossoverValue  = 0; % Y crossover with X axis
hAxis.YRuler.SecondCrossoverValue  = 0; % Y crossover with Z axis
hAxis.ZRuler.FirstCrossoverValue  = 0; % Z crossover with X axis
hAxis.ZRuler.SecondCrossoverValue = 0; % Z crossover with Y axis
hold off

function calculateIntersection(A1, A2, axis1, axis2)
    B = [1;1];
    X1 = A1\B;
    P1 = calculateP(X1,axis1);
    X2 = A2\B;
    P2 = calculateP(X2,axis2);
    P1P2 = [P1, P2];
    line(P1P2(1,:), P1P2(2,:), P1P2(3,:), 'Color',...
    'black', 'linewidth',1.5); view(3)
    function P = calculateP(X, axis)
        if(axis == "x")
            P = [0;X];
        elseif(axis == "y")
            P = [X(1); 0; X(2)];
        else
            P = [X; 0];
        end
    end
end

In a similar way, you can create your own custom function based on your requirements and understanding.

Check out this documentation to learn more about MATLAB function - Declare function name, inputs, and outputs - MATLAB function (mathworks.com)