Visualising 3 times 3 linear algebraic systems and their intersections
2 次查看(过去 30 天)
显示 更早的评论
I want to interpret geometrically
linear systems. Consider for example the system
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1295220/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1295225/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1295230/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1295235/image.png)
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
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.
1 个评论
Bjorn Gustavsson
2023-2-14
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.
采纳的回答
Shushant
2023-3-14
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)
0 个评论
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!