Chương 3: Ma trận
A - M-file
1. Mydet: tinh dinh thuc cua 1 ma tran vuong theo dinh nghia quy nap --- P1
2. X=CramerM(A,B) --- P1
3. x= Gauss(A,b) – b la vecto cot phai --- P1
4. GaussM : B la ma tran ve phai --- P2
5. Gauss Seidel --- P3
6. singIter --- P4
7. MyInv - Ma trận nghịch đảo --- P5
8. [L, U] = crout --- P5
B - Hàm nội trú --- P6
1. Tính giá trị (radian)
2. Tính góc alpha
3. PP phân tích LU
4. Giải hệ pt
--------------------------------------------------------------------------------
1. Mydet: tinh dinh thuc cua 1 ma tran vuong theo dinh nghia quy nap
% tinh dinh thuc cua 1 ma tran vuong A theo dinh nghia qui nap
function d = Mydet(A)
[m n] = size(A);
if m ~= n
disp('???Error using ==> Mydet');
disp('Input Matrix must be square');
d = NaN;
else if n == 1
d = A;
else d = 0;
sig = 1;
for k = 1:n
d = d + sig * A(1,k) * Mydet(A(2:n,[1:k-1,k+1:n]));
sig = -sig;
end
end
end
2. X=CramerM(A,B)
% Giai he phuong trinh Ax=B bang phuong phap Cramer
function X = Cramer(A,B)
[m n] = size(A)
if m ~= n
disp('???Error using ==> Cramer');
disp('Input Matrix must be Square');
else
[p q] = size(B)
if p ~= n
disp('???Error using ==> Cramer');
disp('Dimensions of A & B must be agree');
else
detA = det(A);
for k = 1:n
D = A;
D(:,k) = B;
X(k) = det(D) / detA;
end
X = X.';
end
end
3. x= Gauss(A,b) – b la vecto cot phai
%tinh nghiem phuong trinh tuyen tinh su dung gauss
function X=gauss(A,b)
[m,n]=size(A);
X=NaN*ones(size(b));
if(m~=n)
disp(' A phai la mang vuong');
else
[p,q]=size(b);
if(p~=n)
disp('mang b phai co so dong bang so dong mang A')
else
%% du lieu da chuan
D=[A b];
for k=1:n
%%tim dong co vi tri thu k la lon nhat
[m,r]=max(abs(D(k:n,k)));
r=r+k-1; %%chuyen r(la vi tri trong D(..)ve vi tri trong D
if(r~=k)%% doi cho dong r va k
temp=D(k,:);
D(k,:)=D(r,:);
D(r,:)=temp;
end
D(k,:)=D(k,:)/D(k,k); %%dua so hang (k,k) =1
for i=[1:k-1,k+1:n]
D(i,:)=D(i,:)-D(i,k)*D(k,:);
end
X=D(:,n+1:n+q);
end
end
end
4. GaussM : B la ma tran ve phai
function x=Gauss1(A,B);
[ha ca] = size(A);
[hb cb] = size(B);
if ca~=hb
disp('Error using Gauss!');
x = NaN*ones(size(B));
else
if ha~=ca
disp('Error using Gauss! Matrix A is not square');
x = NaN*ones(size(B));
end;
P = [A B]; [rp cp] = size(P);
for i=1:rp
for j=i+1:rp
P(j,:) = P(j,:)*P(i,i) - P(i,:)*P(j,i);
end;
end;
A = P(:,1:ca);
B = P(:,ca+1:cp);
x=[];
for i = 1:cb
C = [A B(:,i)];
j = ha;
while j>0
if j==ha
X(j) = C(j,ca+1)/C(j,j);
else
S = 0; k = ca;
while k>j
S = S + X(k)*C(j,k);
k = k - 1;
end;
X(j) = (C(j,ca+1)-S)/C(j,j);
end;
j = j - 1;
end;
x=[x X'];
end;
end;
end
5. Gauss Seidel
function x = GaussSeidel(A,b)
[m,n] = size(A);
if(m~=n)
disp('ma tran A phai vuong');
else
[p,q] = size(b);
if(p~=n)
disp('so hang cua ma tran b phai bang so cot cua ma tran A');
else
d1 = 0; d2 = 0;
for i = 1:n
V1 = sum(A(i,:)) - A(i,i);
if V1 < A(i,i)
d1 = d1+1;
end
V2 = sum(A(:,i)) - A(i,i);
if V2 < A(i,i)
d2 = d2+1;
end
end
if d1 ~= n && d2 ~= n
disp('Ma tran khong phai ma tran cheo troi!');
elseif d1 == n && d2 ~= n
for i = 1:n
for j = 1:n
if i ~= j
B(i,j) = -A(i,j)/A(i,i);
else
B(i,j) = 0;
end
end
g(i) = b(i)/A(i,i);
end
q = norm(B);
g = g.';
x = g;
x1 = realmax*x;
while norm(x - x1)*q/(1-q) > 1e-12
x1 = x;
for k = 1:n
x(k) = g(k);
for j = 1:n
x(k) = x(k) + B(k,j)*x(j);
end
end
end
elseif d2 == n
for i = 1:n
for j = 1:n
if i ~= j
B(i,j) = -A(i,j)/A(j,j);
else
B(i,j) = 0;
end
end
g(i) = b(i);
end
q = norm(B);
g = g.';
x = g;
x1 = realmax*x;
while norm(x - x1)*q/(1-q) > 1e-12
x1 = x;
for k = 1:n
x(k) = g(k);
for j = 1:n
x(k) = x(k) + B(k,j)*x(j);
end
end
end
for i = 1:n
x(i) = x(i)/A(i,i);
end
end
end
end
end
6. x= singIter(A,b,p) – Ma trận hệ số A là chéo trội
function x = singiter(A,b)
[m,n] = size(A);
if(m~=n)
disp('ma tran A phai vuong');
else
[p,q] = size(b);
if(p~=n)
disp('so hang cua ma tran b phai bang so cot cua ma tran A');
else
for i=1:n
for j=1:n
B(i,j) = -A(i,j)/A(i,i);
end
B(i,i) = 0;
g(i) = b(i)/A(i,i);
end
end
end
q = norm(B);
g = g.';
x = g;
x1 = realmax * x;
while norm(x - x1)*q/(1-q) > 1e-12
x1 = x;
x = B*x1 + g;
end
x=x';
end
7. X=MyInv(A) Tính ma trận nghịch đảo của 1 ma trận A
function X=MyInv(A)
[m,n]=size(A);
X=NaN;%*ones(size(b));
if(m~=n)
disp(' A phai la ma tran vuong');
else
b=eye(n);
%% du lieu da chuan
D=[A b];
for k=1:n
%%tim dong co vi tri thu k la lon nhat
[m,r]=max(abs(D(k:n,k)));
r=r+k-1; %%chuyen r(la vi tri trong D(..)ve vi tri trong D
if(r~=k)%% doi cho dong r va k
temp=D(k,:);
D(k,:)=D(r,:);
D(r,:)=temp;
end
D(k,:)=D(k,:)/D(k,k); %%dua so hang (k,k) =1
for i=[1:k-1,k+1:n]
D(i,:)=D(i,:)-D(i,k)*D(k,:);
end
X=D(:,n+1:n+n);
end
end
8. [L,U] = Crout(A) – Phân tích ma trận A dưới dạng tích ma trận tam giác trên và dưới
%[L,U]=Crout(A) Phan tich ma tran vuong A thanh tich 2 ma tran tam giac
%duoi L va tren U theo thu tuc Crout: A = L x U
function [L,U]=Crout(A)
[m n]=size(A);
L=NaN;U=NaN;
if m~=n
disp('??? Error using ==>Crout');
disp('Input Matrix must be square');
else
p=eye(n);B=A;
for k=2:n
[M,r]=max(abs(B(k:n,k)));
r=r+k-1;
if r~=k
Lk=B(k,:);
B(k,:)=B(r,:);
B(r,:)=Lk;
Lk=p(k,:);
p(k,:)=p(r,:);
p(r,:)=Lk;
end
B(k+1:n,k)=B(k+1:n,k)/B(k,k);
B(k+1:n,k+1:n)=B(k+1:n,k+1:n)-B(k+1:n,k)*B(k,k+1:n);
end
U=triu(B);
L=tril(B,-1)+eye(n);
End
B – Hàm nội trú
1. Tính giá trị (radian) của góc nhọn alpha tạo bởi 2 vector
>> x=[2 4 5];
>> y=[-2 -5 7];
>> alpha=acos(x*y')/(norm(x)*norm(y))
alpha =
0 + 0.0521i
2. Tính tan alpha với alpha là góc nhọn tạo bởi 2 vector
>> x=[1 3 5];
>> y=[2 5 8];
>> alpha=acos(x*y')/(norm(x)*norm(y))
alpha =
0 + 0.0830i
>> tg=tan(alpha)
tg =
0 + 0.0828i
3. Giải hệ phương trình bằng phương pháp phân tích LU
>> A = [1 2 3 -2; 2 -1 -2 -3; 3 2 -1 2; 2 -3 2 1];
>> B = [6; 8; 4; -8];
>> X = A\B
X =
1.0000
2.0000
-1.0000
-2.0000
4. Giải hệ pt
>> A=[2 -1 0 3;1 1 2 -1;-1 2 3 1;0 1 2 1];
>> b=[1;2;-1;-2];
>> [n m]=size(b);
>> [Q,R,P]=qr(A);
>> b=Q'*P*b;
>> for k=n:-1:1
x(k)=b(k);
for j=k+1:n
x(k)=x(k)-R(k,j)*x(j);
end
x(k)=x(k)/R(k,k);
end
>> x=x
x =
14.5000 -2.5000 -9.0000 -24.5000
Bạn đang đọc truyện trên: AzTruyen.Top