点击下载
本文文档

当前位置：首页 - 正文

LU分解高斯消元列主元高斯消元matlab代码

来源：动视网责编：小OO 时间：2025-10-01 10:13:04

LU分解高斯消元列主元高斯消元matlab代码

数学实验作业一、矩阵LU分解：function[L,U,p]=lutx(A)[n,n]=size(A);p=(1:n)';fork=1:n-1[r,m]=max(abs(A(k:n,k)));m=m+k-1;if(A(m,k)~=0)if(m~=k)A([km],:)=A([mk],:);p([km])=p([mk]);endi=k+1:n;A(i,k)=A(i,k)/A(k,k);j=k+1:n;A(i,j)=A(i,j)-A(i,k)*A(k,j);endendL=tril(A,-1)+e

推荐度：

点击下载本文 文档为doc格式

导读数学实验作业一、矩阵LU分解：function[L,U,p]=lutx(A)[n,n]=size(A);p=(1:n)';fork=1:n-1[r,m]=max(abs(A(k:n,k)));m=m+k-1;if(A(m,k)~=0)if(m~=k)A([km],:)=A([mk],:);p([km])=p([mk]);endi=k+1:n;A(i,k)=A(i,k)/A(k,k);j=k+1:n;A(i,j)=A(i,j)-A(i,k)*A(k,j);endendL=tril(A,-1)+e

数学实验作业

一、矩阵LU分解：

function [L,U,p]=lutx(A)

[n,n]=size(A);

p=(1:n)';

for k=1:n-1

[r,m]=max(abs(A(k:n,k)));

m=m+k-1;

if (A(m,k)~=0)

if (m~=k)

A([k m],:)=A([m k],:);

p([k m])=p([m k]);

end

i=k+1:n;

A(i,k)=A(i,k)/A(k,k);

j=k+1:n;

A(i,j)=A(i,j)-A(i,k)*A(k,j);

end

L=tril(A,-1)+eye(n,n)

U=triu(A)

p

end

高斯消元法求解方程：

n=3;

a=[1 2 3 ;4 5 6 ;7 8 9 ];

b=[17 18 19];

l=eye(n);

y=1;

for i=1:(n-1)

y=y+1;

end

sum=0;

for j=1:n

end

sum=0;

for j=1:n

x(j)=0;

end

for k=n:-1:1

for j=1:n

sum=sum+x(j)*a(k,j)

end

x(k)=(b(k)-sum)/a(k,k)

sum=0;

end

列主元高斯消元法代码：

n=3;

a=[1 2 3 ;4 5 6 ;7 8 9 ];

b=[17 18 19];

l=eye(n);

p=eye(n);

ma=0

for i=1:(n-1)

for j=i:n

if a(j,i)>ma;

ma=a(j,i)

end

for k=i:n

if a(k,i)==ma

m=k;

end

for j=1:n

a1=a(m,j);

a(m,j)=a(i,j);

a(i,j)=a1

p1=p(m,j);

p(m(1),j)=p(i,j);

p(i,j)=p1;

end

b1=b(m);

b(m)=b(i);

b(i)=b1;

ma=0;

end

y=1;

for i=1:(n-1)

for j=1:(n-i)

if a(j+(i-1)*n+y)~=0

l(j+(i-1)*n+y)=a(j+(i-1)*n+y)/a(j+(i-1)*n+y-j)

for k=1:(n-i+1)

a(j+(i-1)*n+y+(k-1)*n)=a(j+(i-1)*n+y+(k-1)*n)-a(j+(i-1)*n+y+(k-1)*n-j)*l(j+(i-1)*n+y)

end

b(j+y-1)=b(j+y-1)-b(y)*l(j+(i-1)*n+y);

end

y=y+1;

end

sum=0;

for j=1:n

x(j)=0;

end

for k=n:-1:1

for j=1:n

sum=sum+x(j)*a(k,j)

end

x(k)=(b(k)-sum)/a(k,k)

sum=0;

end

全主元高斯消元法代码：

n=3;

a=[1 2 3 ;4 5 6 ;7 8 9 ];

b=[17 18 19];

l=eye(n);

p=eye(n);

q=eye(n);

max=0;

for i=1:(n-1)

for j=i:n

for k=i:n

if max max=abs(a(j,k));

end

for j=i:n

for k=i:n

if max==abs(a(j,k))

m=[j,k];

end

for j=1:n

a1=a(m(1),j);

a(m(1),j)=a(i,j);

a(i,j)=a1;

p1=p(m(1),j);

p(m(1),j)=p(i,j);

p(i,j)=p1;

end

b1=b(m(1));

b(m(1))=b(i);

b(i)=b1;

for j=1:n

a1=a(j,m(2));

a(j,m(2))=a(j,i);

a(j,i)=a1;

q1=q(j,m(2));

q(j,m(2))=q(j,i);

q(j,i)=q1;

end

max=0;

end

y=1;

for i=1:(n-1)

for j=1:(n-i)

if a(j+(i-1)*n+y)~=0

l(j+(i-1)*n+y)=a(j+(i-1)*n+y)/a(j+(i-1)*n+y-j)

for k=1:(n-i+1)

a(j+(i-1)*n+y+(k-1)*n)=a(j+(i-1)*n+y+(k-1)*n)-a(j+(i-1)*n+y+(k-1)*n-j)*l(j+(i-1)*n+y);

end

b(j+y-1)=b(j+y-1)-b(y)*l(j+(i-1)*n+y);

end

y=y+1;

end

sum=0;

for j=1:n

x(j)=0;

end

for k=n:-1:1

for j=1:n

sum=sum+x(j)*a(k,j)

end

x(k)=(b(k)-sum)/a(k,k)

sum=0;

end

解：编写矩阵：

0 1 0 0 0 -1 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0 0

a 0 0 -1 -a 0 0 0 0 0 0 0 0

a 0 -1 0 -a 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 -1 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0

A= 0 0 0 0 a 1 0 0 -a -1 0 0 0

0 0 0 0 a 0 1 0 -a 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0 -1

0 0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 1 a 0 0 -a 0

0 0 0 0 0 0 0 0 a 0 1 a 0

0 0 0 0 0 0 0 0 0 0 0 a 1

B= (0 10 0 0 0 0 0 15 0 20 0 0 0)’

F= (f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 f13)’

AF=B F=A\B

程序及运算结果：

>> a=sym(1/sqrt(2))

a =

2^(1/2)/2

>> A=[0

0

a

0

]

A =

[ 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0]

[ 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

[ 2^(1/2)/2, 0, 0, -1, -2^(1/2)/2, 0, 0, 0, 0, 0, 0, 0, 0]

[ 2^(1/2)/2, 0, -1, 0, -2^(1/2)/2, 0, 0, 0, 0, 0, 0, 0, 0]

[ 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0]

[ 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0]

[ 0, 0, 0, 0, 2^(1/2)/2, 1, 0, 0, -2^(1/2)/2, -1, 0, 0, 0]

[ 0, 0, 0, 0, 2^(1/2)/2, 0, 1, 0, -2^(1/2)/2, 0, 0, 0, 0]

[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1]

[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]

[ 0, 0, 0, 0, 0, 0, 0, 1, 2^(1/2)/2, 0, 0, -2^(1/2)/2, 0]

[ 0, 0, 0, 0, 0, 0, 0, 0, 2^(1/2)/2, 0, 1, 2^(1/2)/2, 0]

[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2^(1/2)/2, 1]

>> B=[0;10;0;0;0;0;0;15;0;20;0;0;0]

B =

>> F=A\\B

F =

10*2^(1/2)

-15*2^(1/2)

-5*2^(1/2)

LU分解：

>> a=2^(1/2)/2

a =

>> A=[0

0

a

0

]

A =

Columns 1 through 8

0 1.0000 0 0 0 -1.0000 0 0

0 0 1.0000 0 0 0 0 0

0.7071 0 0 -1.0000 -0.7071 0 0 0

0.7071 0 -1.0000 0 -0.7071 0 0 0

0 0 0 1.0000 0 0 0 -1.0000

0 0 0 0 0 0 1.0000 0

0 0 0 0 0.7071 1.0000 0 0

0 0 0 0 0.7071 0 1.0000 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

-0.7071 -1.0000 0 0 0

-0.7071 0 0 0 0

0 1.0000 0 0 -1.0000

0 0 1.0000 0 0

0.7071 0 0 -0.7071 0

0.7071 0 1.0000 0.7071 0

0 0 0 0.7071 1.0000

>> [L,U,P]=lu(A)

L =

Columns 1 through 8

1.0000 0 0 0 0 0 0 0

0 1.0000 0 0 0 0 0 0

0 0 1.0000 0 0 0 0 0

1.0000 0 -1.0000 1.0000 0 0 0 0

0 0 0 0 1.0000 0 0 0

0 0 0 0 1.0000 1.0000 0 0

0 0 0 0 0 0 1.0000 0

0 0 0 1.0000 0 0 0 1.0000

0 0 0 0 0 0 0 -1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

1.0000 0 0 0 0

0 1.0000 0 0 0

0 0 1.0000 0 0

1.0000 0 1.0000 1.0000 0

0 0 0 0.5000 1.0000

U =

Columns 1 through 8

0.7071 0 0 -1.0000 -0.7071 0 0 0

0 1.0000 0 0 0 -1.0000 0 0

0 0 1.0000 0 0 0 0 0

0 0 0 1.0000 0 0 0 0

0 0 0 0 0.7071 1.0000 0 0

0 0 0 0 0 -1.0000 1.0000 0

0 0 0 0 0 0 1.0000 0

0 0 0 0 0 0 0 -1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

-0.7071 -1.0000 0 0 0

0 1.0000 0 0 0

0 0 0 0 0

0.7071 0 0 -0.7071 0

0 1.0000 0 0 -1.0000

0 0 1.0000 0 0

0 0 0 1.4142 0

0 0 0 0 1.0000

P =

0 0 1 0 0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 0 0 0 0 0 1

>> F=U\\(L\\B)

F =

-15.0000

-10.6066

解：

程序及运算结果：

Lutx.m文件：

function [L,U,p,sig] = lutx(A)

%LU Triangular factorization

% [L,U,p,sig] = lutx(A) computes a unit lower triangular

% matrix L, an upper triangular matrix U, a permutation

% vector p, and a scalar sig, so that L*U = A(p,:) and

% sig = +1 or -1 if p is an even or odd permutation.

[n,n] = size(A);

p = (1:n)';

w=0

for k = 1:n-1

% Find largest element below diagonal in k-th column

[r,m] = max(abs(A(k:n,k)));

m = m+k-1;

% Skip elimination if column is zero

if (A(m,k) ~= 0)

% Swap pivot row

if (m ~= k)

A([k m],:) = A([m k],:);

p([k m]) = p([m k]);

w=w+1;

end

% Compute multipliers

i = k+1:n;

A(i,k) = A(i,k)/A(k,k);

% Update the remainder of the matrix

j = k+1:n;

A(i,j) = A(i,j) - A(i,k)*A(k,j);

end

% Separate result

L = tril(A,-1) + eye(n,n)

U = triu(A)

p

sig=(-1)^w

mydet.m文件：

function det=mydet(A)

[L,U,p,sig] = lutx(A)

det=sig*prod(diag(U))

运行结果：

>> a=2^(1/2)/2

a =

>> A=[0

0

a

0

]

A =

Columns 1 through 8

0 1.0000 0 0 0 -1.0000 0 0

0 0 1.0000 0 0 0 0 0

0.7071 0 0 -1.0000 -0.7071 0 0 0

0.7071 0 -1.0000 0 -0.7071 0 0 0

0 0 0 1.0000 0 0 0 -1.0000

0 0 0 0 0 0 1.0000 0

0 0 0 0 0.7071 1.0000 0 0

0 0 0 0 0.7071 0 1.0000 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

-0.7071 -1.0000 0 0 0

-0.7071 0 0 0 0

0 1.0000 0 0 -1.0000

0 0 1.0000 0 0

0.7071 0 0 -0.7071 0

0.7071 0 1.0000 0.7071 0

0 0 0 0.7071 1.0000

>> mydet(A)

w =

L =

Columns 1 through 8

1.0000 0 0 0 0 0 0 0

0 1.0000 0 0 0 0 0 0

0 0 1.0000 0 0 0 0 0

1.0000 0 -1.0000 1.0000 0 0 0 0

0 0 0 0 1.0000 0 0 0

0 0 0 0 1.0000 1.0000 0 0

0 0 0 0 0 0 1.0000 0

0 0 0 1.0000 0 0 0 1.0000

0 0 0 0 0 0 0 -1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

1.0000 0 0 0 0

0 1.0000 0 0 0

0 0 1.0000 0 0

1.0000 0 1.0000 1.0000 0

0 0 0 0.5000 1.0000

U =

Columns 1 through 8

0.7071 0 0 -1.0000 -0.7071 0 0 0

0 1.0000 0 0 0 -1.0000 0 0

0 0 1.0000 0 0 0 0 0

0 0 0 1.0000 0 0 0 0

0 0 0 0 0.7071 1.0000 0 0

0 0 0 0 0 -1.0000 1.0000 0

0 0 0 0 0 0 1.0000 0

0 0 0 0 0 0 0 -1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

-0.7071 -1.0000 0 0 0

0 1.0000 0 0 0

0 0 0 0 0

0.7071 0 0 -0.7071 0

0 1.0000 0 0 -1.0000

0 0 1.0000 0 0

0 0 0 1.4142 0

0 0 0 0 1.0000

p =

sig =

L =

Columns 1 through 8

1.0000 0 0 0 0 0 0 0

0 1.0000 0 0 0 0 0 0

0 0 1.0000 0 0 0 0 0

1.0000 0 -1.0000 1.0000 0 0 0 0

0 0 0 0 1.0000 0 0 0

0 0 0 0 1.0000 1.0000 0 0

0 0 0 0 0 0 1.0000 0

0 0 0 1.0000 0 0 0 1.0000

0 0 0 0 0 0 0 -1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

1.0000 0 0 0 0

0 1.0000 0 0 0

0 0 1.0000 0 0

1.0000 0 1.0000 1.0000 0

0 0 0 0.5000 1.0000

U =

Columns 1 through 8

0.7071 0 0 -1.0000 -0.7071 0 0 0

0 1.0000 0 0 0 -1.0000 0 0

0 0 1.0000 0 0 0 0 0

0 0 0 1.0000 0 0 0 0

0 0 0 0 0.7071 1.0000 0 0

0 0 0 0 0 -1.0000 1.0000 0

0 0 0 0 0 0 1.0000 0

0 0 0 0 0 0 0 -1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

-0.7071 -1.0000 0 0 0

0 1.0000 0 0 0

0 0 0 0 0

0.7071 0 0 -0.7071 0

0 1.0000 0 0 -1.0000

0 0 1.0000 0 0

0 0 0 1.4142 0

0 0 0 0 1.0000

p =

sig =

ans =

解：

程序及运算结果：

function [L,U,p] = lutx1(A)

[n,n] = size(A);

p = (1:n)';

for k = 1:n-1

[r,m] = max(abs(A(k:n,k)));

m = m+k-1;

if (A(m,k) ~= 0)

if (m ~= k)

A([k m],:) = A([m k],:);

p([k m]) = p([m k]);

end

for i=k+1:n

A(i,k)=A(i,k)/A(k,k);

end

for j=k+1:n

A(i,j)=A(i,j)-A(i,k)*A(k,j);

end

L=tril(A,-1)+eye(n,n)

U=triu(A)

运行结果：

>> lutx1(A)

L =

1 0 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0 0

1 0 -1 0 1 0 0 0 0 0 0 0 0

0 0 0 0 -1 1 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 -1 0 1 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 -1 0 1 1 0

0 0 0 0 0 0 0 0 0 0 0 1 1

U =

Columns 1 through 8

0.7071 0 0 -1.0000 -0.7071 0 0 0

0 1.0000 0 0 0 -1.0000 0 0

0 0 1.0000 0 0 0 0 0

0 0 0 1.0000 0 0 0 -1.0000

0 0 0 0 -0.7071 0 0 0

0 0 0 0 0 1.0000 0 0

0 0 0 0 0 0 1.0000 0

0 0 0 0 0 0 0 1.0000

0 0 0 0 0 0 0 0

Columns 9 through 13

0 0 0 0 0

-0.7071 -1.0000 0 0 0

0 0 0 0 0

0.7071 0 0 -0.7071 0

-0.7071 0 0 0 0

0 1.0000 0 0 -1.0000

0 0 1.0000 0 0

0 0 0 0.7071 0

0 0 0 0 1.0000

ans =

1 0 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0 0

1 0 -1 0 1 0 0 0 0 0 0 0 0

0 0 0 0 -1 1 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 -1 0 1 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 -1 0 1 1 0

0 0 0 0 0 0 0 0 0 0 0 1 1

LU分解高斯消元列主元高斯消元matlab代码

数学实验作业一、矩阵LU分解：function[L,U,p]=lutx(A)[n,n]=size(A);p=(1:n)';fork=1:n-1[r,m]=max(abs(A(k:n,k)));m=m+k-1;if(A(m,k)~=0)if(m~=k)A([km],:)=A([mk],:);p([km])=p([mk]);endi=k+1:n;A(i,k)=A(i,k)/A(k,k);j=k+1:n;A(i,j)=A(i,j)-A(i,k)*A(k,j);endendL=tril(A,-1)+e

推荐度：

点击下载本文 文档为doc格式

热门焦点

LU分解高斯消元列主元高斯消元matlab代码

LU分解高斯消元列主元高斯消元matlab代码

LU分解高斯消元列主元高斯消元matlab代码

最新推荐

猜你喜欢

热门推荐