matlab高斯列主元消去法
matlab高斯消元法
matlab高斯消元法
高斯消元法(Gaussian elimination)是一种求解线性方程组的方法,也常用于求解矩阵的逆以及计算矩阵的秩等问题。
在MATLAB 中,可以通过调用 `rref` 函数来实现高斯消元法。
`rref` 函数的用法如下:
```
X = rref(A)
```
其中,`A` 是要进行高斯消元的矩阵,`X` 是结果。
`rref` 函数会将矩阵 `A` 转化为行最简形,并返回结果矩阵 `X`。
以下是一个使用高斯消元法求解线性方程组的示例:
```matlab
A = [2, 1, -1, 8; -3, -1, 2, -11; -2, 1, 2, -3];
X = rref(A);
```
执行上述代码后,`X` 将会是一个行最简形的矩阵,表示线性方程组的解。
另外,如果你只想求解线性方程组的解而不需要得到行最简形矩阵,可以使用 `linsolve` 函数。
```matlab
A = [2, 1, -1; -3, -1, 2; -2, 1, 2];
b = [8; -11; -3];
X = linsolve(A, b);
```
上述代码中,`A` 是系数矩阵,`b` 是常数向量,`X` 是解向量。
执行代码后,`X` 将会是线性方程组的解。
matlab高斯-约旦消去法
matlab高斯-约旦消去法
高斯-约旦消去法是一种线性代数中的消元法,常用于求解线性方程组。
该方法通过矩阵的初等变换将方程组转化为阶梯型矩阵,从而求解出未知数的值。
具体步骤如下:
设有n个未知数,m个方程,方程组的系数矩阵为A,右端常数为b。
1. 将系数矩阵A和右端常数b组合成增广矩阵Ab。
2. 从第一行开始,将该行的第一个非零元素(称为主元)作为消元元素,用该元素将下面所有行的对应列元素消为零。
3. 重复以上步骤,依次将每一行的主元素作为消元元素,直到将整个矩阵消成阶梯型矩阵。
4. 倒序回代,求出每个未知数的值。
以上就是高斯-约旦消去法的主要步骤。
在实际应用中,需要注意判断矩阵是否可逆,以及主元素是否为零等情况,以保证求解的正确性。
高斯列主元消去法matlab
高斯列主元消去法matlab高斯列主元消去法是一种数值方法,用于求解线性方程组。
下面是使用MATLAB实现高斯列主元消去法的代码示例:```matlabfunction x = gauss_elimination(A, b)% 输入参数:% A: 系数矩阵% b: 右侧常数向量% 输出参数:% x: 方程组的解向量[m, n] = size(A);if m ~= nerror('系数矩阵必须是一个方阵!');end% 增广矩阵AB = [A, b];for k = 1:n-1% 找到主元所在的行[~, pivot_row] = max(abs(AB(k:end, k)));pivot_row = pivot_row + k - 1;% 交换当前行和主元所在行AB([k, pivot_row], :) = AB([pivot_row, k], :);% 消元for i = k+1:nfactor = AB(i, k) / AB(k, k);AB(i, k+1:end) = AB(i, k+1:end) - factor * AB(k, k+1:end);endend% 回代x = zeros(n, 1);x(n) = AB(n, n+1) / AB(n, n);for i = n-1:-1:1x(i) = (AB(i, n+1) - AB(i, i+1:n) * x(i+1:n)) / AB(i, i);endend```使用时,可以调用该函数并传入系数矩阵和右侧常数向量,例如:```matlabA = [1, 2, 3; 4, 5, 6; 7, 8, 9];b = [10; 20; 30];x = gauss_elimination(A, b);```该代码将返回线性方程组的解向量`x`。
gauss列主元素消去法matlab
高斯列主元素消去法是一种解线性方程组的常用方法,特别在数值分析和线性代数中应用广泛。
在Matlab中,我们可以使用该方法来解决大规模的线性方程组,包括矩阵的求解和矩阵的反转。
一、高斯列主元素消去法的基本原理高斯列主元素消去法是一种基于矩阵消元的方法,它通过一系列的矩阵变换将原始的线性方程组转化为上三角形式,然后再进行回代求解。
这个方法的核心就是通过矩阵的变换来简化原始的线性方程组,使得求解过程更加简单高效。
在Matlab中,我们可以利用矩阵运算和函数来实现高斯列主元素消去法,如`lu`分解函数和`\"`运算符等。
通过这些工具,我们能够快速地求解各种规模的线性方程组并得到准确的结果。
二、高斯列主元素消去法在Matlab中的实现在Matlab中,我们可以通过调用`lu`函数来实现高斯列主元素消去法。
该函数返回一个上三角矩阵U和一个置换矩阵P,使得PA=LU。
通过对U进行回代求解,我们可以得到线性方程组的解。
除了`lu`函数之外,Matlab还提供了一些其他的函数和工具来帮助我们实现高斯列主元素消去法,比如`\"`运算符和`inv`函数等。
通过这些工具的组合使用,我们能够更加灵活地进行线性方程组的求解,并且可以方便地处理特殊情况和边界条件。
三、高斯列主元素消去法的应用与局限性高斯列主元素消去法在实际应用中具有广泛的适用性,特别是对于大规模的线性方程组或者稀疏矩阵的求解。
通过Matlab中的工具和函数,我们可以快速地求解各种规模的线性方程组,并得到高精度的数值解。
然而,高斯列主元素消去法也存在一些局限性,比如对于奇异矩阵或者接近奇异矩阵的情况时,该方法的求解精度可能会下降。
在实际应用中,我们需要结合具体的问题和矩阵特性来选择合适的求解方法,以确保得到准确的结果。
四、个人观点和总结作为一种经典的线性方程组求解方法,高斯列主元素消去法在Matlab 中具有较好的实现和应用效果。
通过对其原理和实现细节的深入理解,我们能够更加灵活地应用该方法,并且能够更好地理解其适用性和局限性。
LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序
一、实验目的及题目1.1 实验目的:(1)学会用高斯列主元消去法,LU 分解法,Jacobi 迭代法和Gauss-Seidel 迭代法解线性方程组。
(2)学会用Matlab 编写各种方法求解线性方程组的程序。
1.2 实验题目:1. 用列主元消去法解方程组:1241234123412343421233234x x x x x x x x x x x x x x x ++=⎧⎪+-+=⎪⎨--+=-⎪⎪-++-=⎩2. 用LU 分解法解方程组,Ax b =其中4824012242412120620266216A --⎛⎫ ⎪- ⎪= ⎪ ⎪-⎝⎭,4422b ⎛⎫ ⎪ ⎪= ⎪- ⎪-⎝⎭3. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组:1232341231234102118311210631125x x x x x x x x x x x x x -+=-⎧⎪-+=-⎪⎨-+=⎪⎪-+-+=⎩二、实验原理、程序框图、程序代码等2.1实验原理2.1.1高斯列主元消去法的原理Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数,从而将方程组化为等价形式:1111221122222n n n n nn n nb x b x b x g b x b x g b x g +++=⎧⎪++=⎪⎨⎪⎪=⎩这个过程就是消元,然后再回代就好了。
具体过程如下: 对于1,2,,1k n =-,若()0,k kk a ≠依次计算()()(1)()()(1)()()/,,1,,k k ik ik kk k k k ij ij ik kjk k k i i ik k m a a a a m a b b m b i j k n++==-=-=+然后将其回代得到:()()()()()1/()/,1,2,,1n n n n nn n k k k k k kj j kk j k x b a x b a x a k n n =+⎧=⎪⎨=-=--⎪⎩∑以上是高斯消去。
列主元的高斯消去法
A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n);
b(i)=b(i)*b(k);
end
end
%回代过程
x(n)=b(n)/A(n,n);
for k=n-1:-1:1;
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)')/A(k,k);
end
x=x';
end
index = find(A(:,k)==-ak);
end
%交换列主元
temp = A(index,:);
A(index,:) = A(k,:);
A(k,:) = temp;
temp = b(index);b(index) = b(k); b(k) = temp;
%消元过程
for i=k+1:n
m=A(i,k)/A(k,k);
列主元的高斯消去法列主元高斯消去法列主元消去法高斯消去法matlab高斯消去法全主元消去法高斯消元法高斯消元高斯消元matlab高斯乔丹消元法
列主元的高斯消去法
利用列主元的高斯消去法matlab程序源代码:
首先建立一个gaussMethod.m的文件,用来实现列主元的消去方法。
function x=gaussMethod(A,b)
然后调用gaussMethod函数,来实现列主元的高斯消去法。在命令框中输入下列命令:
输出结果如下:
%高斯列主元消去法,要求系数矩阵非奇异的,%
n = size(A,1);
if abs(det(A))<= 1e-8
error('系数矩阵是奇异的');
return;
Gauss消去matlab程序
%消元过程 for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end end
贵州师范大学数学与计算机科学学院
%回代过程 if abs(A(n,n))<1e-10 index=0; return; end for k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); end
贵州师范大学数学与计算机科学学院
Gauss列主元消去法 列主元消去法Matlab程序 (2) 列主元消去法 程序
function X=Gausslzxq(A,B) % Input A is an N×N nonsingular matrix % B is an N×1 vector % Output X is an N×1 matrix containing the solution to AX=B % Initialize X and the temporary storage matrix C [N N]=size(A); X=zeros(N,1); C=zeros(1,N+1); Az=[A B]; % 形成增广矩阵: Az=[A | B]
贵州师范大学数学与计算机科学学院
if a_max<1e-10 index=0; end if r>k %交换两行 for j=k:n z=A(k,j); A(k,j)=A(r,j); A(r,j)=z; end z=b(k); b(k)=b(r); b(r)=z; end return;
贵州师范大学数学与计算机科学学院
高斯消去、追赶法matlab
⾼斯消去、追赶法matlab 1. 分别⽤Gauss消去法、列主元Gauss消去法、三⾓分解⽅法求解⽅程组程序:(1)Guess消去法:function x=GaussXQByOrder(A,b)%Gauss消去法N = size(A);n = N(1);x = zeros(n,1);for i=1:(n-1)for j=(i+1):nif(A(i,i)==0)disp('对⾓元不能为0');return;endm = A(j,i)/A(i,i);A(j,i:n)=A(j,i:n)-m*A(i,i:n);b(j)=b(j)-m*b(i);endendx(n)=b(n)/A(n,n);for i=n-1:-1:1x(i)=(b(i)-sum(A(i,i+1:n)*x(i+1:n)))/A(i,i);end命令⾏输⼊:A=[1 -1 2 1;-1 3 0 -3 ;2 0 9 -6;1 -3 -6 19];b=[1 3 5 7];b=b';x=GaussXQByOrder(A,b)运算结果:x =-8.0000000000000000.3333333333333333.6666666666666672.000000000000000(2)列主元Gauss消去法程序:function x=GaussXQLineMain(A,b)%列主元Gauss消去法N = size(A);n = N(1);x = zeros(n,1);zz=zeros(1,n);for i=1:(n-1)[~,p]=max(abs(A(i:n,i)));zz=A(i,:);A(i,:)=A(p+i-1,:);A(p+i-1,:)=zz;temp=b(i);b(i)=b(i+p-1);b(i+p-1)=temp;for j=(i+1):nm = A(j,i)/A(i,i);A(j,i:n)=A(j,i:n)-m*A(i,i:n);b(j)=b(j)-m*b(i);endendx(n)=b(n)/A(n,n);for i=n-1:-1:1x(i)=(b(i)-sum(A(i,i+1:n)*x(i+1:n)))/A(i,i);end命令⾏:A=[1 -1 2 1;-1 3 0 -3 ;2 0 9 -6;1 -3 -6 19];b=[1 3 5 7];b=b';x=GaussXQLineMain(A,b)运⾏结果:x =-8.0000000000000050.3333333333333323.6666666666666682.000000000000000(3)三⾓分解⽅法程序:function x = LU(A,b)%三⾓分解N = size(A);n = N(1);L = eye(n,n);U = zeros(n,n);x = zeros(n,1);y = zeros(n,1);U(1,1:n) = A(1,1:n);L(1:n,1) = A(1:n,1)/U(1,1);for k=2:nfor i=k:nU(k,i) = A(k,i)-L(k,1:(k-1))*U(1:(k-1),i);endfor j=(k+1):nL(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k); endendy(1)=b(1)/L(1,1);for i=2:ny(i)=b(i)-sum(L(i,1:i-1)*y(1:i-1));endx(n)=y(n)/U(n,n);for i=n-1:-1:1x(i)=(y(i)-sum(U(i,i+1:n)*x(i+1:n)))/U(i,i);end命令⾏:A=[1 -1 2 1;-1 3 0 -3 ;2 0 9 -6;1 -3 -6 19];b=[1 3 5 7];b=b';x=LU(A,b)运⾏结果:x =-8.0000000000000000.3333333333333333.6666666666666672.000000000000000程序:function [times,wucha]=zhuiganfa(a,b,c,f)%追赶法:x为所求解,times为所有乘除运算次数(即时间),wucha为误差的2-范数。
高斯消元法(含MATLAB编程)
第2次选列主元后的增广矩阵
1 6 5 6 6 0 11/ 3 4 13/ 3 11 0 5/ 3 5 11/ 3 4 2 0 1 0 0
第2次消元后的增广矩阵
1 6 5 6 6 0 11/ 3 4 13 / 3 11 0 0 75 /11 62 /11 9 0 0 24 /11 37 /11 6
(1)输入增广矩阵A=[-3 2 6 4;10 -7 0 7;5 -1 5 6] 第1次选列主元后的增广矩阵 10 -3 -7 2 6 0 4 6 7 61/10 5/2 7
第1次消元后的增广矩阵 5 -1 5 10 -7 0 0 0 -1/10 5/2 6 5
第2次选列主元后的增广ห้องสมุดไป่ตู้阵 10 -7 0 0 0 5/2 -1/10 5 6
3 2 1 1 4 3 2 1 1 4 3/4 4 3 2 1 3/4 7/4 3/2 5/4 1/4 ( 2) A : b 1/2 3 4 3 1 1/2 6/7 4 3 1 4 1 1/4 2 3 4 1 1/4 5/7 3 3 2 1 1 4 3 2 1 1 4 3/4 7/4 3/2 5/4 1/4 3/4 7/4 3/2 5/4 1/4 , 1/2 6/7 12/7 10/7 -12/7 1/2 6/7 12/7 10/7 -12/7 1 1/4 5/7 5/6 5/3 0 1/4 5/7 5/6 4 2 1 1 4 3 3/ 4 1 7 / 4 3/ 2 5/ 4 ,U L 1/ 2 6 / 7 1 12 / 7 10 / 7 1 5/ 3 1/ 4 5 / 7
matlab 列Gauss消去法
列Gauss 消去法进行列主元Gauss 消去法的矩阵表示形式 U A Q L Q L Q L n n =⋅⋅⋅--112211对于某电路的分析,归结为求解线性方程组RI=V 。
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=10; 29 2- 0 0 0 9- 0 0 0 7 2- 27 0 0 5- 0 0 0 0 7- 0 0 41 30- 0 0 0 0 0 120 0 30- 47 7- 0 0 0 0 20- 0 5- 0 7- 57 30- 0 0 0 0 9- 0 0 0 30- 79 10- 0 0 23- 0 0 0 0 0 10- 31 9- 0 27 0 0 0 0 11- 0 9- 35 13- 15- 0 0 0 10- 0 0 0 13- 31R T T V ]10,7,7,12,20,0,23,27,15[----=(1)编制解n 阶线性方程组Ax=b 的列主元三角分解法的通用程序;(2)用所编制的程序解线性方程组RI=V ,并打印出解向量,保留五位有效数;(3)本编程之中,你提高了哪些编程能力?编程思路列主元Gauss 消去法的基本步骤如下:1.从第一行到第n 行找出增广矩阵R(i,1)中最大的数R(k,1),然后将第一行和第k 行的所有元素互相交换,对第2行到第n 行消去每行第一个数值。
2.从第2行到第n 行找出增广矩阵R(i,2)中最大的数R(k,2),然后将第二行和第k 行的所有元素互相交换,对第2行到第n 行消去每行第二个数值。
3.重复这个过程,直到第n 行,完成全部消元。
4.从第n 行开始计算出I(n),然后回代求解出所有解。
M 文件源代码:function R=fGauss();R=[31 -13 0 0 0 -10 0 0 0 -15;-13 35 -9 0 -11 0 0 0 0 27;⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=----)1()1(1)1(3)1(2)1(1)1(3)1(13)1(33)1(32)1(31)1(2)1(12)1(23)1(22)1(21)1(1)1(11)1(13)1(12)1(11nn nn n n n n n n n n n a a a a a a a a a a a a a a a a a a a a A0 -9 31 -10 0 0 0 0 0 -23;0 0 -10 79 -30 0 0 0 -9 0;0 0 0 -30 57 -7 0 -5 0 -20;0 0 0 0 -7 47 -30 0 0 12;0 0 0 0 0 -30 41 0 0 -7;0 0 0 0 -5 0 0 27 -2 7;0 0 0 -9 0 0 0 -2 29 10; ];m=9; %行数n=m+1; %扩展矩阵的列数%列主元法将矩阵转化成上三角矩阵for i=1:m;k=i; %当前计算到第k行% for count=k:n; %对k行下面的矩阵进行处理%a=abs(R(k,k));j=k+1;kk=k;while j<nb=abs(R(j,k));if b>aa=b;kk=j; %kk是最大的数所在的行else j=j+1;endendif kk~=kfor ii=k:n; %交换两行aa=R(k,ii);R(k,ii)=R(kk,ii);R(kk,ii)=aa;endend% end%l=0;pp=k+1;for p=pp:m;q=k;l=R(p,k)/R(k,k);while q<n+1R(p,q)=R(p,q)-l*R(k,q);q=q+1;end% fprintf('l=%f\n',l);endend% fprintf('%f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\n',R);% fprintf('\n \n \n');%回代过程,将结果存放在I矩阵中I=[0 0 0 0 0 0 0 0 0]; %存放结果的矩阵s=0;t=0;I(1,m)=R(m,n)/R(m,m);for s=(m-1):-1:1t=s+1;sum=0;while t<nsum=sum+R(s,t)*I(1,t);t=t+1;%fprintf('sum=%f\n',sum);endI(1,s)=(R(s,n)-sum)/R(s,s);%fprintf('sum=%f\t R(ss)=%f\n',sum,R(s,s));endfprintf('%f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\n',I);(1)方程的解向量为[]T-0.3454360.7128120.289234=-0.220609--I0.2902290.0578230.2010540.4304000.154309-(2)编写此题的感想:在本题的编程中,使我对for、while等循环结构的分析使用更为熟练,对于循环里面套循环的分析思路更为清晰。
高斯列主元消去法
问题提出:采用高斯列主元消去法解线性方程组。
算法(公式)推导:高斯顺序消去法有一个最大的缺点就是一旦对角元素为0,就进行不下去了,为了解决这个问题就有了高斯主元消去法。
如果在高斯顺序消去法消去过程进行到第i 步时,先选取a ri ()n r i ≤≤中(即第i 列)绝对值最大的元素,设为第j 行的元素aji ,然后将第i+1行至第n 行中的每一行减去第i 行乘以ii kj a a (k 代表行号),依次进行消元,这样得到的算法叫高斯按列主元消去法。
高斯按列主元消去法的算法步骤介绍如下:1. 将方程组写成以下的增广矩阵的形式: 432144434241343332312423222114131211b b b b a a a a a a a a a a a a a a a a 2. 对k=1,2,3,…..,n-1,令∑==nk s sk pk a a max ,交换增广矩阵的第k 行与第p 行;对j=k+1,K+2,……..,n,计算*km jkjm jm kk a a a a a =-(m=k,k+1,....n)kk jk k j j a a b b b *-=算法结束。
3. 在MATLABE 中编程实现的高斯按列主元消去法函数为:GaussXQLineMain功能:高斯按列主元消去法求线性方程组Ax=b 的解调用格式:[x,XA]=GaussXQLineMain(A,b)其中,A :线性方程组的系数矩阵;B:线性方程组中的常数向量;x:线性方程组的解:XA:消元后的系数矩阵(可选的输出参数)。
高斯列主元消去法用MATLAB实现如下所示:4.其中用到上三角矩阵求解函数:在MATLABE中编程实现的上三角系数矩阵求解函数为:SolveUPTriangle 功能:求上三角系数矩阵的线性方程组Ax=b的解调用格式:x=SolveUpTriangel(A,b)其中,A :线性方程组的系数矩阵;b :线性方程组中的常数向量; X :线性方程组的解;上三角系数矩阵求解函数用MATLAB 实现如下所示:高斯按列主元消去法解线性方程组应用实例:用高斯按列主元消去法求解下列线性方程组的解。
LU分解高斯消元列主元高斯消元matlab代码
数学实验作业一、矩阵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]);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)+eye(n,n)U=triu(A)pend高斯消元法求解方程: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)for j=1:(n-i)if a(j+(i-1)*n+y)~=0l(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) endb(j+y-1)=b(j+y-1)-b(y)*l(j+(i-1)*n+y);endendy=y+1;endsum=0;for j=1:nsum=sum+x(j)+a(k,j);endsum=0;for j=1:nx(j)=0;endfor k=n:-1:1for j=1:nsum=sum+x(j)*a(k,j)endx(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=0for i=1:(n-1)for j=i:nif a(j,i)>ma;ma=a(j,i)endendfor k=i:nif a(k,i)==mam=k;endendfor j=1:na1=a(m,j);a(m,j)=a(i,j);a(i,j)=a1p1=p(m,j);p(m(1),j)=p(i,j);p(i,j)=p1;endb1=b(m);b(m)=b(i);b(i)=b1;ma=0;endy=1;for i=1:(n-1)for j=1:(n-i)if a(j+(i-1)*n+y)~=0l(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)endb(j+y-1)=b(j+y-1)-b(y)*l(j+(i-1)*n+y);endendy=y+1;endsum=0;for j=1:nx(j)=0;endfor k=n:-1:1for j=1:nsum=sum+x(j)*a(k,j)endx(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:nfor k=i:nif max<abs(a(j,k)) max=abs(a(j,k));endendendfor j=i:nfor k=i:nif max==abs(a(j,k)) m=[j,k];endendendfor j=1:na1=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;endb1=b(m(1));b(m(1))=b(i);b(i)=b1;for j=1:na1=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;endmax=0;endy=1;for i=1:(n-1)for j=1:(n-i)if a(j+(i-1)*n+y)~=0l(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);endb(j+y-1)=b(j+y-1)-b(y)*l(j+(i-1)*n+y);endendy=y+1;endsum=0;for j=1:nx(j)=0;endfor k=n:-1:1for j=1:nsum=sum+x(j)*a(k,j)endx(k)=(b(k)-sum)/a(k,k)sum=0;end解:编写矩阵:0 1 0 0 0 -1 0 0 0 0 0 0 00 0 1 0 0 0 0 0 0 0 0 0 0a 0 0 -1 -a 0 0 0 0 0 0 0 0a 0 -1 0 -a 0 0 0 0 0 0 0 00 0 0 1 0 0 0 -1 0 0 0 0 00 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 00 0 0 0 a 0 1 0 -a 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 -10 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 1 a 0 0 -a 00 0 0 0 0 0 0 0 a 0 1 a 00 0 0 0 0 0 0 0 0 0 0 a 1B= (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 1 0 0 0 -1 0 0 0 0 0 0 00 0 1 0 0 0 0 0 0 0 0 0 0a 0 0 -1 -a 0 0 0 0 0 0 0 0a 0 -1 0 -a 0 0 0 0 0 0 0 00 0 0 1 0 0 0 -1 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 00 0 0 0 a 1 0 0 -a -1 0 0 00 0 0 0 a 0 1 0 -a 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 -10 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 1 a 0 0 -a 00 0 0 0 0 0 0 0 a 0 1 a 00 0 0 0 0 0 0 0 0 0 0 a -1]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 =101520>> F=A\BF =10*2^(1/2)-101010-1010-15*2^(1/2)520-5*2^(1/2)5LU分解:>> a=2^(1/2)/2a =0.7071>> A=[0 1 0 0 0 -1 0 0 0 0 0 0 00 0 1 0 0 0 0 0 0 0 0 0 0a 0 0 -1 -a 0 0 0 0 0 0 0 0a 0 -1 0 -a 0 0 0 0 0 0 0 00 0 0 1 0 0 0 -1 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 00 0 0 0 a 1 0 0 -a -1 0 0 00 0 0 0 a 0 1 0 -a 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 -10 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 1 a 0 0 -a 00 0 0 0 0 0 0 0 a 0 1 a 00 0 0 0 0 0 0 0 0 0 0 a 1]A =Columns 1 through 80 1.0000 0 0 0 -1.0000 0 00 0 1.0000 0 0 0 0 00.7071 0 0 -1.0000 -0.7071 0 0 00.7071 0 -1.0000 0 -0.7071 0 0 00 0 0 1.0000 0 0 0 -1.00000 0 0 0 0 0 1.0000 00 0 0 0 0.7071 1.0000 0 00 0 0 0 0.7071 0 1.0000 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 01.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.7071 -1.0000 0 0 0-0.7071 0 0 0 00 1.0000 0 0 -1.00000 0 1.0000 0 00.7071 0 0 -0.7071 00.7071 0 1.0000 0.7071 00 0 0 0.7071 1.0000>> [L,U,P]=lu(A)L =Columns 1 through 81.0000 0 0 0 0 0 0 00 1.0000 0 0 0 0 0 00 0 1.0000 0 0 0 0 01.0000 0 -1.0000 1.0000 0 0 0 00 0 0 0 1.0000 0 0 00 0 0 0 1.0000 1.0000 0 00 0 0 0 0 0 1.0000 00 0 0 1.0000 0 0 01.00000 0 0 0 0 0 0 -1.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 01.0000 0 0 0 00 1.0000 0 0 00 0 1.0000 0 01.0000 0 1.0000 1.0000 00 0 0 0.5000 1.0000U =Columns 1 through 80.7071 0 0 -1.0000 -0.7071 0 0 00 1.0000 0 0 0 -1.0000 0 00 0 1.0000 0 0 0 0 00 0 0 1.0000 0 0 0 00 0 0 0 0.7071 1.0000 0 00 0 0 0 0 -1.0000 1.0000 00 0 0 0 0 0 1.0000 00 0 0 0 0 0 0 -1.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.7071 -1.0000 0 0 00 1.0000 0 0 00 0 0 0 00 0 0 0 00.7071 0 0 -0.7071 00 1.0000 0 0 -1.00000 0 1.0000 0 00 0 0 1.4142 00 0 0 0 1.0000P =0 0 1 0 0 0 0 0 0 0 0 0 01 0 0 0 0 0 0 0 0 0 0 0 00 1 0 0 0 0 0 0 0 0 0 0 00 0 0 1 0 0 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 00 0 0 0 0 0 0 1 0 0 0 0 00 0 0 0 0 1 0 0 0 0 0 0 00 0 0 0 1 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 0 1 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 00 0 0 0 0 0 0 0 0 0 0 1 00 0 0 0 0 0 0 0 0 0 0 0 1>> F=U\(L\B)F =10.606637.500010.606627.5000-15.000010.606627.5000-10.60667.5000解:程序及运算结果: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=0for 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 zeroif (A(m,k) ~= 0)% Swap pivot rowif (m ~= k)A([k m],:) = A([m k],:);p([k m]) = p([m k]);w=w+1;end% Compute multipliersi = k+1:n;A(i,k) = A(i,k)/A(k,k);% Update the remainder of the matrixj = k+1:n;A(i,j) = A(i,j) - A(i,k)*A(k,j);endend% Separate resultL = tril(A,-1) + eye(n,n)U = triu(A)psig=(-1)^wmydet.m文件:function det=mydet(A)[L,U,p,sig] = lutx(A)det=sig*prod(diag(U))运行结果:>> a=2^(1/2)/2a =0.7071>> A=[0 1 0 0 0 -1 0 0 0 0 0 0 00 0 1 0 0 0 0 0 0 0 0 0 0a 0 0 -1 -a 0 0 0 0 0 0 0 0a 0 -1 0 -a 0 0 0 0 0 0 0 00 0 0 1 0 0 0 -1 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 00 0 0 0 a 1 0 0 -a -1 0 0 00 0 0 0 a 0 1 0 -a 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 -10 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 1 a 0 0 -a 00 0 0 0 0 0 0 0 a 0 1 a 00 0 0 0 0 0 0 0 0 0 0 a 1]A =Columns 1 through 80 1.0000 0 0 0 -1.0000 0 00 0 1.0000 0 0 0 0 00.7071 0 0 -1.0000 -0.7071 0 0 00.7071 0 -1.0000 0 -0.7071 0 0 00 0 0 1.0000 0 0 0 -1.00000 0 0 0 0 0 1.0000 00 0 0 0 0.7071 1.0000 0 00 0 0 0 0.7071 0 1.0000 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 01.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.7071 -1.0000 0 0 0-0.7071 0 0 0 00 1.0000 0 0 -1.00000 0 1.0000 0 00.7071 0 0 -0.7071 00.7071 0 1.0000 0.7071 00 0 0 0.7071 1.0000>> mydet(A)w =L =Columns 1 through 81.0000 0 0 0 0 0 0 00 1.0000 0 0 0 0 0 00 0 1.0000 0 0 0 0 01.0000 0 -1.0000 1.0000 0 0 0 00 0 0 0 1.0000 0 0 00 0 0 0 1.0000 1.0000 0 00 0 0 0 0 0 1.0000 00 0 0 1.0000 0 0 01.00000 0 0 0 0 0 0-1.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 01.0000 0 0 0 00 1.0000 0 0 00 0 1.0000 0 01.0000 0 1.0000 1.0000 00 0 0 0.5000 1.0000U =Columns 1 through 80.7071 0 0 -1.0000 -0.7071 0 0 00 1.0000 0 0 0 -1.0000 0 00 0 1.0000 0 0 0 0 00 0 0 1.0000 0 0 0 00 0 0 0 0.7071 1.0000 0 00 0 0 0 0 -1.0000 1.0000 00 0 0 0 0 0 1.0000 00 0 0 0 0 0 0 -1.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.7071 -1.0000 0 0 00 1.0000 0 0 00 0 0 0 00 0 0 0 00.7071 0 0 -0.7071 00 1.0000 0 0 -1.00000 0 1.0000 0 00 0 0 1.4142 00 0 0 0 1.0000p =31247865119101213sig =-1L =Columns 1 through 81.0000 0 0 0 0 0 0 00 1.0000 0 0 0 0 0 00 0 1.0000 0 0 0 0 01.0000 0 -1.0000 1.0000 0 0 0 00 0 0 0 1.0000 0 0 00 0 0 0 1.0000 1.0000 0 00 0 0 0 0 0 1.0000 00 0 0 1.0000 0 0 01.00000 0 0 0 0 0 0 -1.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 01.0000 0 0 0 00 1.0000 0 0 00 0 1.0000 0 01.0000 0 1.0000 1.0000 00 0 0 0.5000 1.0000U =Columns 1 through 80.7071 0 0 -1.0000 -0.7071 0 0 00 1.0000 0 0 0 -1.0000 0 00 0 1.0000 0 0 0 0 00 0 0 1.0000 0 0 0 00 0 0 0 0.7071 1.0000 0 00 0 0 0 0 -1.0000 1.0000 00 0 0 0 0 0 1.0000 00 0 0 0 0 0 0 -1.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.7071 -1.0000 0 0 00 1.0000 0 0 00 0 0 0 00 0 0 0 00.7071 0 0 -0.7071 00 1.0000 0 0 -1.00000 0 1.0000 0 00 0 0 1.4142 00 0 0 0 1.0000 p =31247865119101213sig =-1ans =-0.5000解:程序及运算结果: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]);endfor i=k+1:nA(i,k)=A(i,k)/A(k,k);endfor j=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endendendL=tril(A,-1)+eye(n,n)U=triu(A)运行结果:>> lutx1(A)L =1 0 0 0 0 0 0 0 0 0 0 0 00 1 0 0 0 0 0 0 0 0 0 0 00 0 1 0 0 0 0 0 0 0 0 0 00 0 0 1 0 0 0 0 0 0 0 0 01 0 -1 0 1 0 0 0 0 0 0 0 00 0 0 0 -1 1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 00 0 0 0 0 0 0 1 0 0 0 0 00 0 0 0 -1 0 1 0 1 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 00 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 0 -1 0 1 1 00 0 0 0 0 0 0 0 0 0 0 1 1 U =Columns 1 through 80.7071 0 0 -1.0000 -0.7071 0 0 00 1.0000 0 0 0 -1.0000 0 00 0 1.0000 0 0 0 0 00 0 0 1.0000 0 0 0 -1.00000 0 0 0 -0.7071 0 0 00 0 0 0 0 1.0000 0 00 0 0 0 0 0 1.0000 00 0 0 0 0 0 01.00000 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0Columns 9 through 130 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.7071 -1.0000 0 0 00 0 0 0 00.7071 0 0 -0.7071 0-0.7071 0 0 0 00 1.0000 0 0 -1.00000 0 1.0000 0 00 0 0 0.7071 00 0 0 0 1.0000ans =1 0 0 0 0 0 0 0 0 0 0 0 00 1 0 0 0 0 0 0 0 0 0 0 00 0 1 0 0 0 0 0 0 0 0 0 00 0 0 1 0 0 0 0 0 0 0 0 01 0 -1 0 1 0 0 0 0 0 0 0 00 0 0 0 -1 1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 00 0 0 0 0 0 0 1 0 0 0 0 00 0 0 0 -1 0 1 0 1 0 0 0 00 0 0 0 0 0 0 0 0 1 0 0 00 0 0 0 0 0 0 0 0 0 1 0 00 0 0 0 0 0 0 0 -1 0 1 1 00 0 0 0 0 0 0 0 0 0 0 1 1。
高斯消元法,列主元素消元法及LU分解法的matlab程序
§2.2.1高斯消元法的MATLAB程序f unction [RA,RB,n,X]=gaus(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.') endend运行命令及结果a=[2.51 1.48 4.53;1.48 0.93 -1.30 ;2.68 3.04 -1.48];b=[0.05;1.03;-0.53];[RA,RB,n,X] =gaus (A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.RA =3RB =3n =3X =1.4531-1.5892-0.2749§2.2.2 列主元素消元法的MATLAB程序function [RA,RB,n,X]=liezhu(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1[Y,j]=max(abs(B(p:n,p))); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.') endend运行命令及结果a=[2.51 1.48 4.53;1.48 0.93 -1.30 ;2.68 3.04 -1.48];b=[0.05;1.03;-0.53];[RA,RB,n,X]=liezhu(A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.RA =3RB =3n =3X =1.4531-1.5892-0.2749§2.2.3 LU分解法的MATLAB程序function hl=zhjLU(A)[n n] =size(A); RA=rank(A);if RA~=ndisp('请注意:因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);returnendif RA==nfor p=1:nh(p)=det(A(1:p, 1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:'), hl;RAreturnendendif h(1,i)~=0disp('请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend运行命令及结果a=[2.51 1.48 4.53;1.48 0.93 -1.30 ;2.68 3.04 -1.48];hl=zhjLU(A)请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:RA =3U =2.5100 1.4800 4.53000 0.9300 -3.97110 0 -0.0837L =1.0000 0 00.5896 1.0000 01.0677 1.5696 1.0000hl =2.5100 0.1439 13.6410>> U=[2.5100 1.4800 4.53000 0.9300 -3.97110 0 -0.0837];>>L= [1.0000 0 00.5896 1.0000 01.0677 1.5696 1.0000];>> b=[0.05;1.03;-0.53];U1=inv(U); L1=inv(L); X=U1*L1*b,x=A\bX =-111.8440110.953125.7324x =1.4531-1.5892-0.2749例2.1: 用高斯消元法求解下面的非齐次线性方程组。