列主元三角分解法在matlab中的实现

合集下载

MATLAB计算方法3解线性方程组计算解法

MATLAB计算方法3解线性方程组计算解法

直到(n-1) 原方程组化为
a11 x1 a12 x2 a1n xn a1,n1 a22 x2 a2 n xn a2 ,n1

ann xn an ,n1
(上三角方程组) (3.2) 以上为消元过程。
(n) 回代求解公式
a n ,n1 xn a nn n x k 1 [a k ,n1 a kj x j ] a kk j k 1 ( k n 1, n 2,...,1)
由矩阵乘法 (1) 1) l11 a11 l11
umj 1 ukj a kj ukj a kj l km umj
m 1
k 1
2 求L的第k列:用L的第i行 u的第k列
(i k 1, , n),即 ( l i 1 , , l ik , l kk , 0 0) ( u1k , u2 k , , ukk , 0 0)' a ik
( 2) 1)求u的第2行:用L的第2行 u的第j列 (j 2, , n) l 21 u1 j 1 u2 j a 2 j u2 j a 2 j l 21u1 j 2)求L的第2列:用L的第i行 u的第2列 (i 3,4, , n) l i 1 u12 l i 2 u22 a i 2 l i 2 (a i 2 l i 1 u12 ) / u22
m 1
l
k 1
im
umk l ik ukk a ik
k 1
l ik a ik l im umk ukk m 1
LU分解式: u1 j a1 j ( j 1,2, n) l i 1 a i 1 u11 ( i 2,3, , n) k 1 ukj a kj l km umj a kj m 1 ( j k , k 1, , n) k 1 l ik a ik l im umk ukk a ik m 1 ( i k 1, , n) ( k 2, 3, , n )

matlab线性方程组的矩阵解法

matlab线性方程组的矩阵解法
2 解 Ux = −1 7 2 0 0 1 -1 x1 1 x1 13 1 0 x 2 2 x2 2 ,得 2 = 13 x = − 13 ,得 x 0 −1 3 7 3 7 x x4 −1 0 11 4 − 11 7 7 0
function x=lupdsv(A,b) n=length(b); [LU,p]=lupd(A); y(1)=b(p(1)); for i=2:n y(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)'; end x(n)=y(n)/LU(n,n); for i=(n-1):-1:1 x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)')/LU(i,i); end
lupdsv.m %功能:调用列主元三角分解函数 [LU,p]=lupd(A) % 求解线性方程组Ax=b。 。 求解线性方程组
%解法:PA=LU, Ax=b←→PAx=Pb 解法: 解法 % % LUx=Pb, Ly=f=Pb, y=Ux f(i)=b(p(i))
%输入:方阵A,右端项 (行或列向量均可) 输入:方阵 ,右端项b(行或列向量均可) 输入 %输出:解x(行向量) 输出: 输出 (行向量)
Ax = d 用矩阵表示 应用追赶法求解三对角线性方程组。追赶法仍然 追赶法求解三对角线性方程组 应用追赶法求解三对角线性方程组。追赶法仍然 保持LU分解特性,它是一种特殊的LU分解。 LU分解特性 LU分解 保持LU分解特性,它是一种特殊的LU分解。充分利用 了系数矩阵的特点,而且使之分解更简单, 了系数矩阵的特点,而且使之分解更简单,得到对三对 角线性方程组的快速解法。 角线性方程组的快速解法。

数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告

数值线性代数第二版徐树方高立张平文上机知识题第一章实验报告

1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序;然后用你编写的程序求解84阶方程组;最后将你的计算结果与方程的精确解进行比较,并就此谈谈你对Gauss 消去法的看法。

Sol :(1)先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序,得到P U L ,,: 不选主元Gauss 消去法:[])(,A GaussLA U L =得到U L ,满足LU A = 列主元Gauss 消去法:[])(,,A GaussCol P U L =得到P U L ,,满足LU PA = (2)用前代法解()Pb or b Ly =,得y用回代法解y Ux =,得x求解程序为()P U L b A Gauss x ,,,,=(P 可缺省,缺省时默认为单位矩阵) (3)计算脚本为ex1_1 代码%算法1.1.3(计算三角分解:Gauss 消去法) function [L,U]=GaussLA(A) n=length(A); for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); endL=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end%算法1.2.2(计算列主元三角分解:列主元Gauss消去法)function [L,U,P]=GaussCol(A)n=length(A);for k=1:n-1[s,t]=max(abs(A(k:n,k)));p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); elsebreak;endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n)); P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A));endn=length(A);b=P*b;for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;endex1_1clc;clear;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1); b=[7;15*ones(82,1);14];%不选主元Gauss消去法[L,U]=GaussLA(A);x1_1=Gauss(A,b,L,U);%列主元Gauss消去法[L,U,P]=GaussCol(A);x1_2=Gauss(A,b,L,U,P);%解的比较subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss'); subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss'); subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精确解');结果为(其中Gauss 表示不选主元的Gauss 消去法,PGauss 表示列主元Gauss消去法,精确解为[]'⨯8411,,1 ):由图,显然列主元消去法与精确解更为接近,不选主元的Gauss 消去法误差比列主元消去法大,且不如列主元消去法稳定。

matlab有限元三角形单元编程

matlab有限元三角形单元编程

matlab有限元三角形单元编程
在MATLAB中进行有限元分析,可以使用其自带的有限元分析工具箱(FEATool)进行编程。

以下是一个简单的例子,演示如何使用三角形单元进行有限元分析:
1. 打开MATLAB,进入FEATool环境。

2. 创建新的有限元模型。

选择“Model”菜单下的“New Model”选项,进入“Model Builder”界面。

3. 在“Model Builder”界面中,选择“2D Triangle”单元类型,并在绘图区域中绘制出三角形网格。

4. 在“Model Builder”界面中,设置材料属性、边界条件和载荷等参数。

5. 运行有限元分析。

选择“Model”菜单下的“Solve”选项,进行有限元求解。

6. 查看结果。

选择“Model”菜单下的“Results”选项,可以查看位移、应力、应变等结果。

以上是一个简单的例子,演示了如何使用三角形单元进行有限元分析。

在实际应用中,还需要根据具体问题进行详细的建模和计算。

计算方法3-4

计算方法3-4

( 2) 找r2,使 a r2 2 max a i 2 ,
2 i n
对调2 r2 行 .
消元:用a 22把a i 2消为0 ( i 3,4, , n) : ai 2 第2行 a 第i行,则 22 ai 2 a ij a 2 j a a ij 22 (i 3,4, , n;j 2,3, , n 1 )
3 3
Gauss 列主元消去法:
优点 ------ 计算结果更可靠; 缺点 ------ 挑主元花机时更多, 有变动,程序复杂。 次序
x1 ,, xn
用Matlab实现选列主元Gauss消去法解线性方程组 在Matlab程序编辑器中输入:
function x=nagauss2(a,b,flag) %a为系数矩阵;b为右 端列向量;flag若为0,则显示中间过程,否则不显示 if nargin<3,flag=0;end n=length(b); a=[a,b]; % 选主元 for k=1:(n-1) [ap,p]=max(abs(a(k:n,k)));p=p+k-1; if p>k,t=a(k,:); a(k,:)=a(p,:); a(p,:)=t; end % 消元 a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1)); a((k+1):n,k)=zeros(n-k,1);
消 元: 用a11将ai 1 ( i 2,, n)化为零;
ai 1 把 a 第1行,加到第i 行。 11

数值分析实验报告

数值分析实验报告
end
%消元过程
fori=k+1:n
m=A(i,k)/A(k,k);
forj=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
%回代过程
ifabs(A(n,n))<1e-10
flag='failure';return;
*x=(x0,x1….,xn),插值节点
*y=(y0,y1,…,yn);被插函数f(x)在插值节点处的函数值
*t求插值函数Pn(x)在t处的函数值
*返回值 插值函数Pn(x)在t处的函数值
*/
procedureNewton
forj=0to n
d1jyj;
endfor
forj=1to n
fori=j to n
[n,m]=size(A);nb=length(b)
%当方程组行与列的维数不相等时,停止计算,并输出出错信息
ifn~=m
error('The row and columns of matrix A must beepual!');
return;
end
%当方程组与右端项的维数不匹配时,停止计算,并输出错误信息
clear
fprintf('gauss-seidel迭代法')
x1_(1)=0;
x2_(1)=0;
x3_(1)=0;
fori=1:9
x1_(i+1)=7.2+0.1*x2_(i)+0.2*x3_(i);

matlab求解代数方程组解析

matlab求解代数方程组解析

第三讲 Matlab 求解代数方程组理论介绍:直接法+迭代法,简单介绍相关知识和应用条件及注意事项 软件求解:各种求解程序讨论如下表示含有n 个未知数、由n 个方程构成的线性方程组:11112211211222221122n n n n n n nn n na x a x a xb a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪+++=⎩ (1)一、直接法 1.高斯消元法:高斯消元法的基本原理: 在(1)中设110,a ≠将第一行乘以111,k a a -加到第(2,3,,),k k n = 得: (1)(1)(1)(1)11112211(2)(1)(2)22112(2)(2)(2)22n n n n n nn n n a x a x a x b a x a x b a x a x b ⎧+++=⎪++=⎪⎨⎪⎪++=⎩(2)其中(1)(1)1111,.k k aa b b ==再设(2)220,a ≠将(2)式的第二行乘以(2)2(2)22,(3,,)k a k n a -= 加到第k 行,如此进行下去最终得到:(1)(1)(1)(1)11112211(2)(1)(2)22112(1)(1)(1)1,111,1()()n n n n n n n n n n n n n n n n nn n n a x a x a x b a x a x b a x a x b a x b --------⎧+++=⎪++=⎪⎪⎨⎪+=⎪⎪=⎩(3) 从(3)式最后一个方程解出n x ,代入它上面的一个方程解出1n x -,并如此进行下去,即可依次将121,,,,n n x x x x - 全部解出,这样在()0(1,2,,)k kk a k n ≠= 的假设下,由上而下的消元由下而上的回代,构成了方程组的高斯消元法. 高斯消元法的矩阵表示:若记11(),(,,),(,,)T T ij n n n n A a x x x b b b ⨯=== ,则(1)式可表为.Ax b =于是高斯消元法的过程可用矩阵表示为:121121.n n M M M Ax M M M b --=其中:(1)21(1)111(1)1(1)11111n a a M a a ⎛⎫ ⎪ ⎪- ⎪=⎪ ⎪ ⎪ ⎪- ⎪⎝⎭ (2)32(2)222(2)2(2)221111n a a M a a ⎛⎫⎪⎪ ⎪-⎪=⎪ ⎪ ⎪⎪- ⎪⎝⎭高斯消元法的Matlab 程序: %顺序gauss 消去法,gauss 函数 function[A,u]=gauss(a,n) for k=1:n-1%消去过程 for i=k+1:n for j=k+1:n+1%如果a(k,k)=0,则不能削去 if abs(a(k,k))>1e-6 %计算第k 步的增广矩阵 a(i,j)=a(i,j)-a(i,k)/a(k,k)*a(k,j); else%a(k,k)=0,顺序gauss 消去失败 disp (‘顺序gauss 消去失败‘); pause; exit; end end end end%回代过程 x(n)=a(n,n+1)/a(n,n); for i=n-1:-1:1 s=0; for j=i+1:n s=s+a(i,j)*x(j); endx(i)=(a(i,n+1)-s)/a(i,i); end%返回gauss 消去后的增广矩阵 A=triu(a); %返回方程组的解 u=x ;练习和分析与思考: 用高斯消元法解方程组:12345124512345124512452471523814476192536x x x x x x x x x x x x x x x x x x x x x x ++++=⎧⎪+++=⎪⎪++++=⎨⎪+++=⎪+++=⎪⎩2.列主元素消元法在高斯消元法中进行到第k 步时,不论()k ik a 是否为0,都按列选择()||(,,)k ik a i k n = 中最大的一个,称为列主元,将列主元所在行与第k 行交换再按高斯消元法进行下去称为列主元素消元法。

数值分析MATLAB科学计算—线性方程组

数值分析MATLAB科学计算—线性方程组

科学计算—理论、方法及其基于MATLAB的程序实现与分析 三、 解线性方程组(线性矩阵方程)解线性方程组是科学计算中最常见的问题。

所说的“最常见”有两方面的含义:1) 问题的本身是求解线性方程组;2) 许多问题的求解需要或归结为线性方程组的求解。

关于线性方程组B A x B Ax 1-=⇒=(1)其求解方法有两类:1) 直接法:高斯消去法(Gaussian Elimination ); 2) 间接法:各种迭代法(Iteration )。

1、高斯消去法1) 引例考虑如下(梯形)线性方程组:()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路:把一般的线性方程组(1)化成(上或下)梯形的形式。

2)高斯消去法——示例考虑如下线性方程组:⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x 1) 第一个方程的两端乘12加到第二个方程的两端,第一个方程的两端乘-1加到第三个方程的两端,得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x2) 第二个方程的两端乘001.010-加到第三个方程的两端,得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x3) 从上述方程组的第三个方程依此求解,得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x 3)高斯消去法的不足及其改进——高斯(全、列)主元素消去法在上例中,由于建模、计算等原因,系数2.001而产生0.0005的误差,实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x注:数值稳定的算法高斯列主元素消去法就是在消元的每一步选取(列)主元素—一列中绝对值最大的元取做主元素,高斯列主元素消去法是数值稳定的方法。

线性方程组及MATLAB应用

线性方程组及MATLAB应用

数值实验 线性方程组与MATLAB 应用王1.实验目的:理解矩阵的范数与条件数。

实验内容:已知矩阵⎪⎪⎪⎪⎪⎭⎫⎝⎛------=1111111111111111A 求1A ,2A ,∞A 和)(2A cond 。

解:编写了一个M 文件来求矩阵A 的范数与条件数:test3_1.m 如下:A=[1 1 1 1;-1 1 -1 1;-1 -1 1 1;1 -1 -1 1]; norm(A,1) norm(A,2) norm(A,inf) cond(A,2)计算结果依次是: 4 2 4 1.00002.实验目的:研究高斯消去法的数值稳定性(出现小主元)。

实验内容:设方程组b Ax =,其中两个矩阵如下,分别对以上两个方程组(1)⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--⨯=-11212592.1121130.6291.51314.59103.0151A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2178.4617.591b (2)⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=201015152699990999999999.23107102A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=1500019000000000.582b (1)计算矩阵的条件数,判断系数矩阵是良态的还是病态的?解: 本题编写了一个test3_21的M 文件如下:A1=[0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; A2=[10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2]; cond(A1) cond(A2)求得两个矩阵的条件数分别为68.4296和8.9939,易知这矩阵A1的条件数远远大于1,而矩阵A2的条件数刚大于1,故这,矩阵A1为病态矩阵,矩阵A2为良态矩阵。

(2)用列主元消去法求得L 和U 及解向量412,∈R x x ;解:本题利用Matlab 的列主元三角分解函数lu();具体求解如下: >> A1=[0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; >> A2=[10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2];>> b1=[59.17;46.78;1;2];>> b2=[8;5.0000000000001;5;1];>> [L1,U1]=lu(A1)L1 = 0.0000 1.0000 0 00.4724 -0.1755 1.0000 01.0000 0 0 00.0893 0.0202 -0.1738 1.0000 U1 = 11.2000 9.0000 5.0000 2.00000 59.1400 3.0000 1.00000 0 -2.8354 1.23070 0 0 1.0151 >> [L2,U2]=lu(A2)L2 =1.0000 0 0 0 -0.3000 -0.0000 1.0000 00.5000 1.0000 0 00 0.4000 -0.3333 1.0000 U2 =10.0000 -7.0000 0 1.00000 2.5000 5.0000 -1.50000 0 6.0000 2.30000 0 0 3.3667 >> y1=L1\b1;>> x1=U1\y1x1 =3.84571.6095-15.476110.4113>> y2=L2\b2;>> x2=U2\y2x2 =0.1337-0.82180.88420.9109用不选主元的高斯消去法求得L和U及解向量412, Rx x;解:编写一个LU_Fact的M文件储存不选主元的LU分解法然后调用求解。

线性方程组直接解法实验

线性方程组直接解法实验

实验一 线性方程组直接解法实验一、实验目的1.运用matlab 软件完成线性方程组的直接实验;2.通过实验,了解Doolittle 分解方法和列主元消去法解方程组的过程,并比较两种方法的优点。

二、实验题目分别用Doolittle 分解方法和列主元消去法解方程组123410-7018-3 2.09999962 5.9000015-15-1521021⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭x x x x . 输出A ,b ;Doolittle 分解方法的L 和U ;解向量x,det A ;列主元方法的行交换次序,解向量x,det A ;比较两种方法所得的结果。

三、实验原理1) Doolittle 分解方法的原理算法原理:应用高斯消去法解n 阶线性方程Ax b =经过1n -步消去后得出一个等价的上三角形方程组()()n n A x b =,对上三角形方程组用逐步回代就可以求出解来。

这个过程也可通过矩阵分解来实现。

将非奇异阵分解成一个下三角阵L 和上三角阵U 的乘积称为对矩阵A 的三角分解,又称LU 分解。

根据LU 分解,将Ax b =分解为Ly bUx y =⎧⎨=⎩形式,简化了求解问题。

程序框图:变量说明:ij a 为系数矩阵元素,i b 为常数矩阵系数,,ij ij l u 分别为下、上三角矩阵元素。

2)列主元消去法解方程组的原理算法原理:列选主元是当变换到第k步时,从k列的kk a及以下的各元素中选取绝对值a的位置上,然后再进行消元过程。

交换系数矩阵中最大的元素,通过行交换将其交换到kk的两行(包括常数项),相当于两个方程的位置交换了。

程序框图:Array变量说明:k表示消元到a为消元第k步时第k步,kk主对角线元素3)四、实验过程及结果1)Doolittle分解方法的输出结果----------计算实习题----------Page64 第1题用Doolittle分解方法解方程组A =10.0000 -7.0000 0 1.0000-3.0000 2.1000 6.0000 2.00005.0000 -1.0000 5.0000 -1.00002.0000 1.0000 0 2.0000b =8.00005.90005.00001.0000L =1.0e+006 *0.0000 0 0 0-0.0000 0.0000 0 00.0000 -2.5000 0.0000 00.0000 -2.4000 0.0000 0.0000 U =1.0e+007 *0.0000 -0.0000 0 0.00000 -0.0000 0.0000 0.00000 0 1.5000 0.57500 0 0 0.0000 X =-0.0000-1.00001.00001.0000det(A)值为-762.00009000----------输出完毕----------2)列主元消去法输出结果----------计算实习题----------Page64 第1题列主元消去法解方程组A =10.0000 -7.0000 0 1.0000-3.0000 2.1000 6.0000 2.00005.0000 -1.0000 5.0000 -1.00002.0000 1.0000 0 2.0000b =8.00005.90005.00001.0000X =0.0000-1.00001.00001.0000detA值为-762.00009000----------输出完毕----------五、实验分析1.运用LU分解法可以成批地解方程组,且速度快.若c先求LU=A3,再解(LU)x=b,则要重新计算,计算量增加;如果按照上述方法计算,能够减少运算次数,加快运算速度.3. ⑴无论当n=10、n=100、n=1000时,x1与x2的值都相等,且随着n的增大,变化的只是解的中间部分数字,头、前后几位数都没有变化.⑵高斯消去法应用于三对角方程组得到的就是所谓的“追赶法”.追赶法不需要对零元素计算,只有6n-5次乘除法计算量,求解速度快.且当系数矩阵对角占优时数值稳定,是解三对角方程组的优秀解法.⑶用LU分解法解此方程组速度慢.顺序高斯消去法实际上就是将方程组的系数矩阵分解成单位下三角矩阵与上三角矩阵的乘积.顺序高斯消去法的消元过程相当于LU分解过程和Ly=b的求解,回代过程则相当于解线性方程组Ux=y,故其求解速度慢.六、附原程序1)Doolittle分解方法原程序fprintf('----------计算实习题----------\n')fprintf('Page64 第1题用Doolittle分解方法解方程组\n')A=[10 -7 0 1 ; -3 2.099999 6 2 ;5 -1 5 -1 ; 2 1 0 2];b=[8;5.900001;5;1];n=length(A);U=zeros(n,n);L=eye(n,n);U(1,:)=A(1,:);L(2:n,1)=A(2:n,1)/U(1,1);for i=2:n;U(i,i:n)=A(i,i:n)-L(i,1:i-1)*U(1:i-1,i:n);L(i+1:n,i)=(A(i+1:n,i)-L(i+1:n,1:i-1)*U(1:i-1,i))/U(i,i); endY=zeros(n);Y(1)=b(1);for i=2:nY(i)=b(i)-L(i,1:i-1)*Y(1:i-1,1);endX=zeros(n,1);if det(U)==0;X=0;elseX(n)=Y(n)/U(n,n);for i=n-1:-1:1X(i)=(Y(i)-U(i,i+1:n)*X(i+1:n,1))/U(i,i);endendAbLUXfprintf('det(A)值为%9.8f\n',det(A))fprintf('----------输出完毕 ----------\n')2)列主元消去法原程序fprintf('----------计算实习题----------\n')fprintf('Page64 第1题列主元消去法解方程组\n')A=[10 -7 0 1 ; -3 2.099999 6 2 ;5 -1 5 -1 ; 2 1 0 2];b=[8;5.900001;5;1];C=[A b];n=length(A);D=zeros(n,n+1);l=zeros(n,1);for i=1:nD=C;k=min(find(C(i:n,i)==max(C(i:n,i))));C(i,i:n+1)=D(k+i-1,i:n+1);C(k+i-1,i:n+1)=D(i,i:n+1);l(i+1:n,1)=C(i+1:n,i)/C(i,i);C(i+1:n,i:n+1)= C(i+1:n,i:n+1)- l(i+1:n,1)*C(i,i:n+1); endX=zeros(n,1);X(n)=C(n,n+1)/C(n,n);for i=n-1:-1:1X(i)=(C(i,n+1)-C(i,i+1:n)*X(i+1:n,1))/C(i,i); endAbXfprintf('detA值为%9.8f\n',det(A))fprintf('----------输出完毕----------\n')。

MATLAB实验一 解线性方程组的直接法

MATLAB实验一 解线性方程组的直接法

输出 Ax b 中系数 A LU 分解的矩阵 L 和 U ,解向量 x 和 det(A) ;用列主元法 的行交换次序解向量 x 和求 det(A) ;比较两种方法所得结果。
2、用列主高斯消元法解线性方程组 Ax b 。
3.01 6.03 1.99 x1 1 4.16 1.23 x 2 1 (1) 、 1.27 0.987 4.81 9.34 x 1 3 3.00 6.03 1.99 x1 1 4.16 1.23 x 2 1 (2) 、 1.27 0.990 4.81 9.34 x 1 3
index = 1 3、在 MATLAB 窗口:
A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32 23 33 31]'; x=A\b b1=[32.1 22.9 33.1 30.9]'; x1=A\b1 A1=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; x2=A1\b delta_b=norm(b-b1)/norm(b) delta_A=norm(A-A1)/norm(A) delta_x1=norm(x-x1)/norm(x) delta_x2=norm(x-x2)/norm(x)
二. 实验要求 1、按照题目要求完成实验内容; 2、写出相应的 Matlab 程序; 3、给出实验结果(可以用表格展示实验结果); 4、分析和讨论实验结果并提出可能的优化实验。 5、写出实验报告。 三. 实验步骤 1、用 LU 分解及列主元高斯消去法解线性方程组
7 10 3 2.099999 a) 5 1 2 1 1 x1 8 6 2 x 2 5.900001 , 5 1 x3 5 0 2 1 x 4 0

线性代数方程组数值解法及MATLAB实现综述

线性代数方程组数值解法及MATLAB实现综述

线性代数方程组数值解法及MATLAB 实现综述廖淑芳 20122090 数计学院 12计算机科学与技术1班(职教本科) 一、分析课题随着科学技术的发展,提出了大量复杂的数值计算问题,在建立电子计算机成为数值计算的主要工具以后,它以数字计算机求解数学问题的理论和方法为研究对象。

其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。

因此在各种数据处理中,线性代数方程组的求解是最常见的问题之一。

关于线性代数方程组的数值解法一般分为两大类:直接法和迭代法。

直接法就是经过有限步算术运算,可求的线性方程组精确解的方法(若计算过程没有舍入误差),但实际犹如舍入误差的存在和影响,这种方法也只能求得近似解,这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。

直接法包括高斯消元法,矩阵三角分解法、追赶法、平方根法。

迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。

迭代法具有需要计算机的存储单元少,程序设计简单,原始系数矩阵在计算过程始终不变等优点,但存在收敛性级收敛速度问题。

迭代法是解大型稀疏矩阵方程组(尤其是微分方程离散后得到的大型方程组)的重要方法。

迭代法包括Jacobi 法SOR 法、SSOR 法等多种方法。

二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss 消元法通过一系列的加减消元运算,也就是代数中的加减消去法,以使A 对角线以下的元素化为零,将方程组化为上三角矩阵;然后,再逐一回代求解出x 向量。

1.1消元过程1. 高斯消元法(加减消元):首先将A 化为上三角阵,再回代求解。

11121121222212n n n n nn n a a a b a a a b a a a b ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭L L M M O M M L (1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000nn n n n nn n a a a a b a a a b a a b a b ⎛⎫ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎝⎭L L L M M M OM M L 步骤如下:第一步:1111,2,,i a i i n a -⨯+=L 第行第行 11121121222212n nn n nn n a a a b a a a b a a a b ⎛⎫⎪ ⎪ ⎪⎪⎝⎭L L M M O M M L111211(2)(2)(2)2222(2)(2)(2)200n nn nn n a a a b a a b a a b ⎛⎫⎪⎪ ⎪ ⎪⎝⎭LL M M O M M L第二步:(2)2(2)222,3,,i a i i n a -⨯+=L 第行第行111211(2)(2)(2)2222(2)(2)(2)200nnn nn n a a a b a a b a a b ⎛⎫ ⎪⎪ ⎪ ⎪⎝⎭L L M M O M M L11121311(2)(2)(2)(2)222322(3)(3)(3)3333(3)(3)(3)300000n n nn nn n a a a a b a a a b a a b a a b ⎛⎫⎪ ⎪ ⎪⎪⎪ ⎪⎝⎭LL LM M M O M M L 类似的做下去,我们有:第k 步:()()k ,1,,k ikk kka i i k n a -⨯+=+L 第行第行。

matlab计算方法实验指导误差分析

matlab计算方法实验指导误差分析

实验一 误差分析实验1(病态问题)实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。

对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。

通过本实验可获得一个初步体会。

数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。

病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。

问题提出:考虑一个高次的代数多项式)1.1()()20()2)(1()(201∏=-=---=k k x x x x x p显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。

现考虑该多项式的一个扰动)2.1(0)(19=+x x p ε其中ε是一个非常小的数。

这相当于是对(1.1)中19x 的系数作一个小的扰动。

我们希望比较(1.1)和(1.2)根的差别,从而分析方程(1.1)的解对扰动的敏感性。

实验内容:为了实现方便,我们先介绍两个MATLAB 函数:“roots ”和“poly ”。

roots(a)u =其中若变量a 存储n+1维的向量,则该函数的输出u 为一个n 维的向量。

设a 的元素依次为121,,,+n a a a ,则输出u 的各分量是多项式方程01121=+++++-n n n n a x a x a x a的全部根;而函数 poly(v)b =的输出b 是一个n+1维向量,它是以n 维向量v 的各分量为根的多项式的系数。

可见“roots ”和“poly ”是两个互逆的运算函数。

))20:1((;)2();21,1(;000000001.0ve poly roots ess ve zeros ve ess +===上述简单的MATLAB 程序便得到(1.2)的全部根,程序中的“ess ”即是(1.2)中的ε。

实验要求:(1)选择充分小的ess ,反复进行上述实验,记录结果的变化并分析它们。

第三章 MATLAB线性方程组及矩阵特征值

第三章 MATLAB线性方程组及矩阵特征值

求解线性方程组的主要方法有:
直接法:理论,无舍入误差,有限步,精确解 迭代法:格式,无穷序列 解向量 x
§3.1 解线性方程组的直接法 一、 Gauss消去法
设有 a11 x1 a12 x2 a1n xn a1,n1 a21 x1 a22 x2 a2 n xn a2,n1 ai1 x1 ai 2 x2 ain xn ai ,n1 an1 x1 an 2 x2 ann xn an ,n1
(n) 回代求解公式
a xn n,n 1 ann n 1 xk [ak ,n1 akj x j ] akk j k 1 (k n 1, n 2,...,1
(3.3)
(3.3) 是回代过程。
说明: (1)也可采用无回代的列主元消去法(叫Gauss--Jordan消去法),该法同时消去对角元上 下的元素,且仍旧需要选主元,但比有回 代的列主元消去法的乘除运算次数多。 Gauss-Jordan消去法的优点之一是用它来 计算逆矩阵的算法非常容易解释。 (2)有回代的列主元消去法所进行的乘除运算 次数为 1 n3 n 2 1 n ,量很小。
回 代: an,n 1 xn a
nn n 1 xk [ak ,n1 akj x j ] akk j k 1 (k n 1, n 2,...,1
问 题:aii 0或 aii 0?
二、列主元素消去法---计算结果可靠
(1)找行号r1 使 ar11 max ai1 ,对调1 r1行:
(3.1)
线性代数:方法不好时工作量非常大, 工作量小的方法是 Gauss 消去法。

1、克劳特(Crout)(LU)分解法求解线性方程组的matlab实现

1、克劳特(Crout)(LU)分解法求解线性方程组的matlab实现

1、克劳特(Crout)(LU)分解法求解线性方程组的matlab实现1、克劳特(Crout)(LU)分解法求解线性方程组function [x,L,U]=Crout(A,b)%Crout分解法求解线性方程组%系数矩阵:AN=size(A);n=N(1);L=zeros(n,n); %下三角矩阵U=eye(n,n); %上三角矩阵L(1:n,1)=A(1:n,1); %L的第一列U(1,1:n)=A(1,1:n)/L(1,1); %U的第一行?for k=2:nfor i=k:nL(i,k)=A(i,k)-L(i,1:(k-1))*U(1:(k-1),k);%L的第k列endfor j=(k+1):nU(k,j)=(A(k,j)-L(k,1:(k-1))*U(1:(k-1),j))/(L(k,k));%U的第k行endend%y=inv(L)*b;%x=inv(U)*y;y=SolveDownTriangle(L,b);x=SolveUpTriangle(U,y); %求解线性方程组的解x%x=U\(L\b);function x=SolveUpTriangle(A,b)%求解上三角矩阵Ax=b的解N=size(A);n=N(1);for i=n:-1:1if(i<n)< p="">s=A(i,(i+1):n)*x((i+1):n,1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);endfunction x=SolveDownTriangle(A,b) %求解下三角矩阵Ax=b的解N=size(A);n=N(1);for i=1:nif(i>n)s=A(i,1:(i-1))*x(1:(i-1),1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end%求解线性方程组的解clcclearA=[12 -3 3;-16 3 -1;1 1 1];b=[15;-13;6];%x=A\b[x,L,U]=Crout(A,b)解:123L =12.0000 0 0-16.0000 -1.0000 01.0000 1.2500 4.5000U =1.0000 -0.2500 0.25000 1.0000 -3.00000 0 1.0000列主元LU分解function [L,U,x]=lux(A,b)%LU 分解法解线性方程组(列主元LU分解)[n,n]=size(A);p=eye(n);%p记录了选择主元时候所进行的行变换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);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);%解下三角矩阵 Ly=bnewb=p*b;y=zeros(n,1);for k=1:nj=1:k-1;y(k)=(newb(k)-L(k,j)*y(j))/L(k,k); end %解上三角方程组 Ux=yx=zeros(n,1);for k=n:-1:1j=k+1:n;x(k)=(y(k)-U(k,j)*x(j))/U(k,k);end</n)<>。

matlab中的矩阵分解

matlab中的矩阵分解

matlab 中的矩阵分解matlab 中的矩阵分解矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。

常见的矩阵分解有LU分解(三角分解)、QR分解(正交变换)、Cholesky分解,以及Schur分解、Hessenberg分解、奇异分解等。

(1) LU分解(三角分解)矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。

线性代数中已经证明,只要方阵A是非奇异(即行列式不等于0)的,LU分解总是可以进行的。

MATLAB提供的lu函数用于对矩阵进行LU分解,其调用格式为:[L,U]=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换),使之满足X=LU。

注意,这里的矩阵X必须是方阵。

[L,U,P]=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PX=LU。

当然矩阵X同样必须是方阵。

(设P 是一个m×n 的(0,1) 矩阵,如m≤n且P*P′=E,则称P为一个m×n的置换矩阵。

)实现LU分解后,线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb),这样可以大大提高运算速度。

例7-2 用LU分解求解例7-1中的线性方程组。

命令如下:A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]&#39;;[L,U]=lu(A);x=U\(L\b)或采用LU分解的第2种格式,命令如下:[L,U ,P]=lu(A);x=U\(L\P*b)(2) QR分解(正交变换)对矩阵X进行QR分解,就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。

QR 分解只能对方阵进行。

MATLAB的函数qr可用于对矩阵进行QR分解,其调用格式为:[Q,R]=qr(X):产生一个一个正交矩阵Q和一个上三角矩阵R,使之满足X=QR。

[Q,R,E]=qr(X):产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E,使之满足XE=QR。

列主元消去法matlab实验报告

列主元消去法matlab实验报告

列主元消去法matlab实验报告列主元消去法是一种常用的线性方程组求解方法,它通过选取主元元素来消去其他元素,从而简化方程组的求解过程。

本文将以Matlab为工具,对列主元消去法进行实验研究,并给出相应的实验报告。

我们需要明确列主元消去法的基本原理。

列主元消去法的核心思想是选取每一列的主元素,通过消去其他元素,从而将方程组转化为上三角形或下三角形的形式。

具体来说,通过选取第一列的主元素,将第一列下方的元素消去;然后选取第二列的主元素,将第二列下方的元素消去;依此类推,直到最后一列。

这样,我们就得到了一个上(下)三角形的方程组,可以通过回代(代入)的方法求解。

接下来,我们使用Matlab编写代码,实现列主元消去法。

首先,我们需要输入一个线性方程组的系数矩阵A和常数向量b,其中A 是一个n×n的矩阵,b是一个n×1的向量。

然后,我们通过选取主元素的方式进行消去操作,得到一个上三角形的方程组。

最后,我们通过回代(代入)的方法求解方程组的解。

具体实现的代码如下所示:```matlabfunction x = gauss_elimination(A, b)n = size(A, 1); %方程组的个数% 消元过程for k = 1:n-1[~, p] = max(abs(A(k:n, k))); %选取主元素 p = p + k - 1;% 交换第k行和第p行temp = A(k, :);A(k, :) = A(p, :);A(p, :) = temp;temp = b(k);b(k) = b(p);b(p) = temp;% 消去操作for i = k+1:nfactor = A(i, k) / A(k, k);A(i, :) = A(i, :) - factor * A(k, :);b(i) = b(i) - factor * b(k);endend% 回代(代入)过程x = zeros(n, 1);x(n) = b(n) / A(n, n);for i = n-1:-1:1x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i, i);endend```接下来,我们将使用一个具体的例子来说明列主元消去法的求解过程。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

列主元三角分解法在matlab中的实现
摘要:介绍了M atlab语言并给出用M atlab语言实现线性方程组的列主元三角分解法,其有效性已在计算机实现中得到了验证。

关键词:M atlab语言;高斯消去法;列主元三角分解法
0前言
M atlab是M atrix Laboratory(矩阵实验室)的缩写,它是由美国M athwork公司于1967年推出的软件包,现已发展成为一种功能强大的计算机语言。

它编程简单,使用方便,在M a tlab环境下数组的操作与数的操作一样简单,进行数学运算可以像草稿纸一样随心所欲,使计算机兼备高级计算器的优点。

M atlab语言具有强大的矩阵和向量的操作功能,是Fo rtran和C语言无法比拟的;M a tlab语言的函数库可任意扩充;语句简单,内涵丰富;还具有二维和三维绘图功能且使用方便,特别适用于科学和工程计算。

在科学和工程计算中,应用最广泛的是求解线性方程组的解,一般可用高斯消去法求解,如果系数矩阵不满足高斯消去法在计算机上可行的条件,那么消元过程中可能会出现零主元或小主元,消元或不可行或数值不稳定,解决办法就是对方程组进行行交换或列交换来消除零主元或小主元,这就是选主元的思想。

1 定义
列主元三角分解:如果A为非奇异矩阵,则存在排列矩阵P,使PA=LU,其中L为单位下三角矩阵,U为上三角阵。

列主元三角分角法是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,
就是避免小数作为分母项.
2 算法概述
列主元三角分解法和普通三角分解法基本上类似,所不同的是在构造Gauss 变换前,先在对应列中选择绝对值最大的元素(称为列主元),然后实施初等行交换将该元素调整到矩阵对角线上。

例如第)1,,2,1(-=n k 步变换叙述如下:
选主元:确定p 使{}1)1(
max -≤≤-=k ik n
i k k pk a a ; 行交换:将矩阵的第k 行和第p 行上的元素互换位置,即

实施Gauss 变换:通过初行变换,将列主对角线以下的元素消为零.即
3 列主元三角分解在matlab 中的实现
其程序如下:
function [l,u,p]=mylu(A)
[m,n]=size(A);
if m~=n
error('矩阵不是方阵')
return
end
if det(A)==0
error('矩阵不能被三角分解')
end
u=A;p=eye(m);l=eye(m);
for i=1:m
for j=i:m
t(j)=u(j,i);
for k=1:i-1
t(j)=t(j)-u(j,k)*u(k,i);
end
end
a=i;b=abs(t(i));
for j=i+1:m
if b<abs(t(j))
b=abs(t(j));
a=j;
end
end
if a~=i
for j=1:m
c=u(i,j);
u(i,j)=u(a,j);
u(a,j)=c;
end
for j=1:m
c=p(i,j);
p(i,j)=p(a,j);
p(a,j)=c;
end
c=t(a);
t(a)=t(i);
t(i)=c;
end
u(i,i)=t(i);
for j=i+1:m
u(j,i)=t(j)/t(i);
end
for j=i+1:m
for k=1:i-1
u(i,j)=u(i,j)-u(i,k)*u(k,j);
end
end
end
l=tril(u,-1)+eye(m);
u=triu(u,0);。

相关文档
最新文档