Matlab matrix calc - 《MatLab》

Gauss Elimination 函数方案1
Another approach Gauss Elimination 方案2
Cramer’s Rule
LU decomposition
三角追赶法
Gauss Seidel’s
inv() 逆矩阵
eig() 本征值
eigs(Matrix, k) 前k项最大本征值

A=[1 -1 1; 5 -4 3; 2 1 1];
B=[-4 -12 11]
%SAVE AS .M FILE

det（A）

ans =

     2

Gauss Elimination 函数方案1

function  y = gauss(A,b)
n = length(b);
y = zeros(1,n);
for i = 1:n-1
    if(A(i,i)==0)
        t = min(find(A(i+1:n,i)~=0)+i);
        if (isempty(t))
            disp('gauss ERROR: matrix is singular');
            return
        end;
        temp = A(i,:);
        tb = b(i);
        A(i,:)=A(t,:);
        b(i)=b(t);
        A(t,:)=temp;
        b(t)=tb;
    end;
    for j =i+1:n
        m=-A(j,i)/A(i,i);
        A(j,i)=0;
        A(j,i+1:n) = A(j,i+1:n) + m*A(i,i+1:n);
        b(j)=b(j)+m*b(i);
    end;
end;
y(n)=b(n)/A(n,n);
for i =n-1:-1:1
    y(i)=(b(i) - sum(y(i+1:n).*A(i,i+1:n)))/A(i,i);
end;


%SAVE AS gauss.M

Another approach Gauss Elimination 方案2

function x = solveGauss(A,b)
s = length(A);
for j = 1:(s-1)
    for i = s:-1:j+1
        m = A(i,j)/A(j,j);
        A(i,:) = A(i,:) - m*A(j,:);
        b(i) = b(i) - m*b(j);
    end
end 
x = zeros(s,1);
x(s) = b(s)/A(s,s);               
for i = s-1:-1:1                    
    sum = 0;
    for j = s:-1:i+1                
        sum = sum + A(i,j)*x(j);    
    end 
    x(i) = (b(i)- sum)/A(i,i);
end 

%SAVE AS solveGauss.m

When executing

gauss(A,B)

ans =

     3     6    -1


solveGauss(A,B)

ans =

     3
     6
    -1

结果发现一个行一个列，原因是两个程序在初始化答案矩阵的时候，一个给的n by 1（列），一个为1 by n（行）

Cramer’s Rule

function y = cramer(A,b)
[m,n] = size(A);
y = zeros(m,1);
for i = 1:m
    T = A;
    T(:,i) = b;
    y(i)  =det(T)/det(A);
end

%SAVE AS cramer.m

When executing;

cramer(A,B)

ans =

     3
     6
    -1

LU decomposition

function [L,U,x]=luDecom(A,b)
n=length(b);
L=eye(n);
U=zeros(n,n);
x=zeros(n,1);
y=zeros(n,1);

%LU DECOMPOSITION

for i=1:n
    U(1,:)=A(1,:);
    if i ==1
        L(i,1)=1;
    else
        L(i,1)=A(i,1)/U(1,1);
    end
end
for i=2:n
    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;
        if j~=n  % "~=" checks if j and n are NOT equal( j not euqal to n returns 1)
            sum=0;
            for k=1:i-1
                sum=sum+L(j+1,k)*U(k,i);
            end
            L(j+1,i)=(A(j+1,i)-sum)/U(i,i);
        end
    end
end

% Solve Ly=b
y(1)=b(1);
for k=2:n
    sum=0;
    for j=1:k-1
        sum=sum+L(k,j)*y(j);
    end
    y(k)=b(k)-sum;
end

% Solve Ux=y
x(n)=y(n)/U(n,n);
for k=n-1:-1:1
    sum=0;
    for j=k+1:n
        sum=sum+U(k,j)*x(j);
    end
    x(k)=(y(k)-sum)/U(k,k);
end

%SAVE AS luDecom.m

Assigning A=[1 -1 1; 5 -4 3; 2 1 1];
and b=[-4 -12 11].
When executing

[L,U,x]=luDecom(A,b)

L =

     1     0     0
     5     1     0
     2     3     1


U =

     1    -1     1
     0     1    -2
     0     0     5


x =

     3
     6
    -1

注意，在高版本的Matlab中，可以使用 [~,~,x]=luDecom(A,b) 来只显示x的结果，但在Matlab 7.0以下相对复杂，输出全部结果更为简单。

三角追赶法

function x=tri(a,b,c,f)
n=length(f);
x=zeros(1,n);
y=zeros(1,n);
d=zeros(1,n);
u=zeros(1,n);
d(1)=b(1);
for i=1:n-1
    u(i)=c(i)/d(i);
    d(i+1)=b(i+1)-a(i+1)*u(i);
end

%追
y(1)=f(1)/d(1);
for i=2:n
    y(i)=(f(i)-a(i)*y(i-1))/d(i);
end

%赶
x(n)=y(n);
for i=n-1:-1:1
    x(i)=y(i)-u(i)*x(i+1);
end

when executing

tri(a,b,c,f)

ans =

    1.0000    2.0000    3.0000    4.0000    5.0000

Gauss Seidel’s

function [x,k,flag]=gaussSeidel(A,b,w)
n=length(A);
k=0;
x=zeros(n,1);
y=zeros(n,1);
flag='OK';
delta=1e-5;
maxn=100;

while 1
    y=x;
    for i=1:n
        z=b(i);
        for j=1:n
            if j~=i
                z=z-A(i,j)*x(j);
            end
        end
        if abs(A(i,i))<1e-10|k==maxn
            flag='FAIL';
            return;
        end
        z=z/A(i,i);
        x(i)=z;
    end
    if norm(y-x,inf)<delta
        break;
    end
    k=k+1
end

Executing with A=[7 2 1 -2; 9 15 3 -2; -2 -2 11 5; 1 3 2 13]; b=[ 4 7 -1 0], giving

迭代误差随迭代步数的变化曲线: ???

for maxn=1:20
        err=x(1)-0.4979;
    plot(maxn,err);
end

inv() 逆矩阵

H = [2 -1 0; -1 2 -1; 0 -1 2]

inv(H)

When executing

ans =

    0.7500    0.5000    0.2500
    0.5000    1.0000    0.5000
    0.2500    0.5000    0.7500

用 matrixcalc.org) 验证，发现正确。

eig() 本征值

eig(H)

when executing

用 matrixcalc.org) 验证，发现正确

eigs(Matrix, k) 前k项最大本征值

eigs(H,1)

when executing

Iteration 1: a few Ritz values of the 3-by-3 matrix:
     0


ans =

    3.4142

显然合理。