Sử dụng Matlab trong đại số tuyến tính

,j));
 end;
 end;
 if l(g,g)==0
 break;
 end;
 end;
if l(g,g)==0
 break;
end;
end;
end;
E=l(2:N,1:m,:); % bien doi de dinh thuc ma tran khong thay doi
F=l(1,1:m);
F=((-1)^d)*F;
B=[F;E];%ma tran B thay the boi B ma xet su thay doi dau hang 1
hang=0;
for s=1:N %tim hang ma tran
 i=N+1-s;
 for j=1:m
 if B(i,j)~=0
 hang=i;
 break;
 end;
 end;
 if B(i,j)~=0
 break;
 end;
end;
if N==m % tim dinh thuc ma tran
 dinhthuc=1;
 for i=1:N
 dinhthuc=dinhthuc*B(i,i);
 end;
end;
2/ Giải hệ phương trình bằng phương pháp Cramer
function [X] = cramer(A,b)
if nargin <2 error ('phai nhap vao ma tran can phan tich');
end;
N = size(A,1);%1 la dong,2 la cot,i la hang,j la cot
m = size(A,2);
if N~=m error ('khong phai la he cramer,nhap lai ma tran he so A');
end;
n = size(b,1);
M = size(b,2);
if n>M
 b=b';
 if M~=1 error ('kich co he so tu do khong phu hop nhap lai b');
 end;
 if n~=N error ('kich co he so tu do khong phu hop nhap lai b');
 end;
end;
if M>n
 if n~=1 error ('kich co he so tu do khong phu hop nhap lai b');
 end;
 if M~=N error ('kich co he so tu do khong phu hop nhap lai b');
 end;
end; 
if det(A)==0 error ('det=0,khong the dung phuong phap nay');
end;
l=A;
X=0;
b=b';
for j=1:N
 if j==1
 l=A;
 l(:,j)=[];
 l=[b l];
 X(1,j)=det(l)/det(A);
 end;
 if j>=2
 if j<N
 l=A;
 a=l;
 x=a(1:N,1:(j-1),:);
 y=a(1:N,(j+1):N,:);
 l=[x,b];
 l=[l y];
 X(1,j)=det(l)/det(A);
 end; 
 if j==N
 l=A;
 l(:,j)=[];
 l=[l b];
 X(1,j)=det(l)/det(A);
 end;
 end;
end;
3/ Giải hệ bằng phương pháp khử Gauss
function [B,X] = gauss(A,b)
if nargin <1 error ('phai nhap vao ma tran can phan tich');
end;
N = size(A,1);%1 la dong,2 la cot,i la hang,j la cot
m = size(A,2)+1;
z = size(b,1);%kiem tra he so tu do
t = size(b,2);
c=b;
if t>z
 for j=1:t % chuyen vi ma tran
 for i=1:z
 b(j,i)=c(i,j);
 end;
 end;
b=b(:,1:z);
end;
if z>t
 b=c;
end;
z = size(b,1);
t = size(b,2);
%if z~=1
 % if t==N error ('kich co he so tu do khong phu hop');
 % end;
%end;
%if t~=1 error ('kich co he so tu do khong phu hop.');
%end;
l=[A b];
i=1;
d=0;
for h=1:N
k=l(h:N,:);%chon hang thu h den N
k=k(:,h:m); %chon cot thu h den m
for j=1:(m-h+1)
 for i=1:(N-h+1)
 if k(i,j)~=0
 e=i;
 if h==1
 f=j;
 end;
 if h~=1
 f=1;
 end;
 break;
 end;
 end;
 if k(i,j)~=0
 break;
 end;
end;%tim phan tu khac 0 dau tien de hoan doi
 k=circshift(k, [-(e-1), -(f-1)]);
 p=[l(1:(h-1),h:m);k];
 l=[l(1:N,1:(h-1)) p];
for g=h:(N-1)
 for i=(g+1):N
 if l(g,g)~=0
 c=l(i,g)/l(g,g);
 for j=g:m
 l(i,j)=l(i,j)-(c*l(g,j));
 end;
 end;
 if l(g,g)==0
 break;
 end;
 end;
if l(g,g)==0
 break;
end;
end;% bien doi so cap khu cac phan tu ben duoi
end;
B=l;
rA=0;
rB=0;
for s=1:N
 i=N+1-s;
 for j=1:(m-1)
 if B(i,j)~=0
 rA=i;
 break;
 end;
 end;
 if B(i,j)~=0
 break;
 end;
end;%tim hang ma tran he so
for s=1:N
 i=N+1-s;
 for j=1:m
 if B(i,j)~=0
 rB=i;
 break;
 end;
 end;
 if B(i,j)~=0
 break;
 end;
end;% tim hang ma tran mo rong
if rA<rB
 X= ('phuong trinh vo nghiem');
end;
if rA==rB
 if rA<(m-1)
 X= ('phuong trinh co vo so nghiem');
 end;
 u=0;
 if rA==(m-1) %giai nghiem he theo phuong phap khu Gauss
 for s=1:(m-1)
 i=m-s;
 if i==(m-1);
 X(i)=B(i,m)/B(i,i);
 end;
 if i~=(m-1)
 for c=1:(s-1)
 y=m-c;
 u=u+X(y)*B(i,y);
 X(i)=(B(i,m)-u)/B(i,i);
 end;
 end;
 u=0;
 end;
 end;
end;
4/ Tìm ma trận nghịch đảo 
function [B] = nghichdao (A)% tim ma tran nghich dao theo ma tran phu hop
if nargin <1 error ('phai nhap vao ma tran can phan tich');
end;
N = size(A,1);%1 la dong,2 la cot,i la hang,j la cot
m = size(A,2);
d=0;
a=A;
i=1;
if N~=m
 error ('ma tran khong vuong');
end;
l=A;
for h=1:N %bien doi bac thang
 k=A(h:N,:);
 k=k(:,h:N);
 for j=1:(N-h+1)
 for i=1:(N-h+1)
 if k(i,j)~=0
 e=i;
 if h==1
 f=j;
 end;
 if h~=1
 f=1;
 end;
 break;
 end;
 end;
 if k(i,j)~=0
 break;
 end;
 end;
 d=d+(e-1)*(N-1)+(f-1)*(N-1);
 k=circshift(k, [-(e-1), -(f-1)]);
 p=[A(1:(h-1),h:N);k];
 A=[A(1:N,1:(h-1)) p];
 for g=h:(N-1)
 for i=(g+1):N
 if A(g,g)~=0
 c=A(i,g)/A(g,g);
 for j=g:N
 A(i,j)=A(i,j)-(c*A(g,j));
 end;
 end;
 if A(g,g)==0
 break;
 end;
 end;
 if A(g,g)==0
 break;
 end;
 end;
end;
E=A(2:N,1:N,:);
F=A(1,1:N);
F=((-1)^d)*F;
K=[F;E];
dinhthuccha=1;
for i=1:N %tinh dinh thuc
 dinhthuccha=dinhthuccha*K(i,i);
end;
for o=1:N
 l(o,:)=[];
 D=l;
 for t=1:N
 l(:,t)=[]; 
 N=N-1;
 for h=1:N
 k=l(h:N,:);
 k=k(:,h:N);
 for j=1:(N-h+1)
 for i=1:(N-h+1)
 if k(i,j)~=0
 e=i;
 if h==1
 f=j;
 end;
 if h~=1
 f=1;
 end;
 break;
 end;
 end;
 if k(i,j)~=0
 break;
 end;
 end;
 d=d+(e-1)*(N-1)+(f-1)*(N-1);
 k=circshift(k, [-(e-1), -(f-1)]);
 p=[l(1:(h-1),h:N);k];
 l=[l(1:N,1:(h-1)) p];
 for g=h:(N-1)
 for i=(g+1):N
 if l(g,g)~=0
 c=l(i,g)/l(g,g);
 for j=g:N
 l(i,j)=l(i,j)-(c*l(g,j));
 end;
 end;
 if l(g,g)==0
 break;
 end;
 end;
 if l(g,g)==0
 break;
 end;
 end;
 end;
 E=l(2:N,1:N,:);
 F=l(1,1:N);
 F=((-1)^d)*F;
 G=[F;E];
 dinhthuccon=1;
 for i=1:N
 dinhthuccon=dinhthuccon*G(i,i);
 end;
 B(o,t)=((-1)^(o+t))*dinhthuccon;
 dinhthuccon=1;
 l=D;
 N=N+1;
 end;
 l=a;
end;%tim ma tran dinh thuc con
C=B;
for j=1:N % chuyen vi ma tran
 for i=1:N
 B(i,j)=C(j,i);
 end;
end;
B=B/dinhthuccha;
5/Phân tích LU của ma trận A
function [l,u,dinhthuc] = LU(A)%tam giac duoi,tam giac tren,dinh thuc
if nargin <1 error ('phai nhap vao ma tran can phan tich');
end;
a=A;
N = size(A,1);%1 la dong,2 la cot
m = size(A,2);
if N~=m 
 A = ('khong la ma tran vuong');
end;
if N==m
 l=zeros(N);
 u=zeros(N);
 for i=1:N
 l(i,i)=1;
 end;
 for j=1:N
 u(1,j)=A(1,j);
 end;
 if u(1,1)== 0
 error('khong the phan tich duoc, nhap lai ma tran khac');
 end;
 for i=2:N
 l(i,1)=A(i,1)/u(1,1);
 end;
 for i=2:N-1
 for j=i:N
 sum=0;
 for k=1:i-1
 sum=sum+l(i,k)*u(k,j);
 end;
 u(i,j)=A(i,j)-sum;
 end;
 if u(i,i)==0
 error('khong the phan tich duoc, nhap lai ma tran khac'); 
 end;
 for j=i+1:N
 sum=0;
 for k=1:i-1
 sum=sum+l(j,k)*u(k,i);
 end;
 l(j,i)=(A(j,i)-sum)/u(i,i);
 end;
 end;
 sum=0;
 for k=1:N-1
 sum=sum+l(N,k)*u(k,N);
 end;
 u(N,N)=A(N,N)-sum;
 dinhthuc=1;
 for i=1:N
 dinhthuc = dinhthuc*u(i,i)*l(i,i);
 end;
end;
A=a;
6/ Dùng quá trình trực giao Gram-Schmidt đưa họ độc lập tuyến tính về họ trực chuẩn
function orthonormal_basis = Gram_Schmidt_process(A)
matrix_size = size(A);
m = matrix_size(1,1);
n = matrix_size(1,2);
if A == zeros(m,n)
 error('There does not exist any type of basis for the zero vector space.');
elseif n == 1
 orthonormal_basis = A(1:m,1)/norm(A(1:m,1));
else
 flag = 0;
 if is_orthonormal_set(A) == 1
 orthonormal_basis = A;
 flag = 1;
 end
 if flag == 0;
 if rank(A) ~= n
 A = basis_col(A);
 end
 matrix_size = size(A);
 m = matrix_size(1,1);
 n = matrix_size(1,2);
 orthonormal_basis = A(1:m,1)/norm(A(1:m,1));
 for i = 2:n
 u = A(1:m,i);
 v = zeros(m,1);
 for j = 1:(i - 1)
 v = v - dot(u,orthonormal_basis(1:m,j))*orthonormal_basis(1:m,j);
 end
 v_ = u + v;
 orthonormal_basis(1:m,i) = v_/norm(v_);
 end
 end
end
7/ Kiểm tra họ vécto có phải là họ trực giao
function result = is_orthogonal_set(A)
matrix_size = size(A);
m = matrix_size(1,1);
n = matrix_size(1,2);
tolerance = 10^-10;
if n == 1
 result = 1;
else
 orthogonal_counter = 0;
 for i = 1:n
 for j = 1:n
 if i == j
 else
 if abs(dot(A(1:m,i),A(1:m,j))) <= tolerance
 orthogonal_counter = orthogonal_counter + 1;
 end
 end
 end
 end
 if orthogonal_counter == factorial(n)/factorial(n - 2);
 result = 1;
 else
 result = 0;
 end
end
8/ Bài toán bình phương cực tiểu.
function c = fit(n, t, y)
% The least-squares approximating polynomial of degree n (n>=0).
% Coordinates of points to be fitted are stored in the column vectors
% t and y. Coefficients of the approximating polynomial are stored in
% the vector c. Graphs of the data points and the least-squares
% approximating polynomial are also generated.
if ( n >= length(t))
error('Degree is too big')
end
v = fliplr(vander(t));
v = v(:,1:(n+1));
c = v\y;
c = fliplr(c');
x = linspace(min(t),max(t));
w = polyval(c, x);
plot(t,y,'ro',x,w);
title(sprintf('The least-squares polynomial of degree n = %2.0f',n))
legend('data points','fitting polynomial')
Chạy thử chương trình trên bằng cách nhập:
a/ t = linspace(0, pi/2, 10); t = t';
b/ y = sin(2*t);
c/ c = fit(3, t, y) 
9/ Ứng dụng trong lý thuyết mật mã.
Đoạn code dùng mã hóa thông tin
function B = code(s, A)
% String s is coded using a nonsingular matrix A.
% A coded message is stored in the vector B.
p = length(s);
[n,n] = size(A);
b = double(s);
r = rem(p,n);
if r ~= 0
b = [b zeros(1,n-r)]';
end
b = reshape(b,n,length(b)/n);
B = A*b;
B = B(:)';
Đoạn code giải mã thông tin
function s = dcode(B, A)
% Coded message, stored in the vector B, is
% decoded with the aid of the nonsingular matrix A
% and is stored in the string s.
[n,n]= size(A);
p = length(B);
B = reshape(B,n,p/n);
d = A\B;
s = char(d(:)');
Chạy thử chương trình trên:
a/ s=’Toi hoc dai so tuyen tinh’
b/ A=pascal(4) tạo ra ma trận khả nghịch tùy ý( trong trường hợp này dùng mt pascal cấp 4)
c/B=code(s,A)
d/ giải mã dcode(B,A)
10/ Ứng dụng trong lý thuyết mật mã.
function [xt, yt] = rot2d(t, x, y)
% Rotation of a two-dimensional object that is represented by two
% vectors x and y. The angle of rotation t is in the degree measure.
% Transformed vectors x and y are saved in xt and yt, respectively.
t1 = t*pi/180;
r = [cos(t1) -sin(t1);sin(t1) cos(t1)];
x = [x x(1)];
y = [y y(1)];
hold on
grid on
axis equal
fill(x, y,'b')
z = r*[x;y];
xt = z(1,:);
yt = z(2,:);
fill(xt, yt,'r');
title(sprintf('Plane rotation through the angle of %3.2f degrees',t))
hold off
Chạy thử chương trình bằng cách
a/ x=[1 2 3 2]
b/ y=[3 1 2 4]
c/ [xt, yt] = rot2d(75, x, y)