搜档网
当前位置:搜档网 › 求解线性方程组LU分解法的Matlab程序#精选.

求解线性方程组LU分解法的Matlab程序#精选.

求解线性方程组LU分解法的Matlab程序#精选.
求解线性方程组LU分解法的Matlab程序#精选.

function [L,U,flag]=LU_decom(A)

[n,m]=size(A);

if n~=m

error('The rows and columns of matrix A must be equal!');

return;

end

L=eye(n);

U=zeros(n);

flag='OK';

for k=1:n

for j=k:n

z=0;

for q=1:k-1

z=z+L(k,q)*U(q,j);

end

U(k,j)=A(k,j)-z;

end

if abs(U(k,k))

flag='failure';

return;

end

for i=k+1:n

z=0;

for q=1:k-1

z=z+L(i,q)*U(q,k);

end

L(i,k)=(A(i,k)-z)/U(k,k);

end

end

最新文件仅供参考已改成word文本。方便更改如有侵权请联系网站删除

MatLab求解线性方程组

MatLab解线性方程组一文通 当齐次线性方程AX=0,rank(A)=r

MATLAB代码 解线性方程组的迭代法

解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多,可能不收敛!'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;

if(n>=M) disp('Warning:迭代次数太多,可能不收敛!'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多,可能不收敛!'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵

matlab解方程组

matlab解方程组 lnx表示成log(x) 而lgx表示成log10(x) 1-exp(((log(y))/x^0.5)/(x-1)) 1、解方程 最近有多人问如何用matlab解方程组的问题,其实在matlab中解方程组还是很方便的,例如,对于代数方程组Ax=b(A为系数矩阵,非奇异)的求解,MATLAB 中有两种方法: (1)x=inv(A)*b —采用求逆运算解方程组; (2)x=A\B —采用左除运算解方程组 PS:使用左除的运算效率要比求逆矩阵的效率高很多~ 例: x1+2x2=8 2x1+3x2=13 >>A=[1,2;2,3];b=[8;13]; >>x=inv(A)*b x = 2.00 3.00 >>x=A\B x = 2.00 3.00; 即二元一次方程组的解x1和x2分别是2和3。 对于同学问到的用matlab解多次的方程组,有符号解法,方法是:先解出符号解,然后用vpa(F,n)求出n位有效数字的数值解.具体步骤如下: 第一步:定义变量syms x y z ...; 第二步:求解[x,y,z,...]=solve('eqn1','eqn2',...,'eqnN','var1','var2',...'varN'); 第三步:求出n位有效数字的数值解x=vpa(x,n);y=vpa(y,n);z=vpa(z,n);...。 如:解二(多)元二(高)次方程组: x^2+3*y+1=0 y^2+4*x+1=0 解法如下: >>syms x y; >>[x,y]=solve('x^2+3*y+1=0','y^2+4*x+1=0'); >>x=vpa(x,4); >>y=vpa(y,4); 结果是:

MATLAB解线性方程组的直接方法

在这章中我们要学习线性方程组的直接法,特别是适合用数学软件在计算机上求解的方法. 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意:因为RA~=RB ,所以此方程组无解.') return end if RA==RB if RA==n disp('请注意:因为RA=RB=n ,所以此方程组有唯一解.') else disp('请注意:因为RA=RB> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意:因为RA=RB=n ,所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入 >>X=A\b, 运行后输出结果为 X =(0 0 0 0)’. (2) 在MATLAB 工作窗口输入程序 >> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b)

利用MATLAB求线性方程组

《MATLAB语言》课成论文 利用MATLAB求线性方程组 姓名:郭亚兰 学号:12010245331 专业:通信工程 班级:2010级通信工程一班 指导老师:汤全武 学院:物电学院 完成日期:2011年12月17日

利用MATLAB求解线性方程组 (郭亚兰 12010245331 2010 级通信一班) 【摘要】在高等数学及线性代数中涉及许多的数值问题,未知数的求解,微积分,不定积分,线性方程组的求解等对其手工求解都是比较复杂,而MATLAB语言正是处理线性方程组的求解的很好工具。线性代数是数学的一个分支,它的研究对象是向量,向量空间(或称线性空间),线性变换和有限维的线性方程组。因而,线性代数被广泛地应用于抽象代数和泛函分析中;由于科学研究中的非线性模型通常可以被近似为线性模型,使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是数学的一个分支,它的研究对象是向量,向量空间(或称线性空间),线性变换和有限维的线性方程组。因而,线性代数被广泛地应用于抽象代数和泛函分析中;由于科学研究中的非线性模型通常可以被近似为线性模型,使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是讨论矩阵理论、与矩阵结合的有限维向量空间及其线性变换理论的一门学科。 【关键字】线性代数MATLAB语言秩矩阵解 一、基本概念 1、N级行列式A:A等于所有取自不同性不同列的n个元素的积的代数和。 2、矩阵B:矩阵的概念是很直观的,可以说是一张表。 3、线性无关:一向量组(a1,a2,…,an)不线性相关,既没有不全为零的数 k1,k2,………kn使得:k1*a1+k2*a2+………+kn*an=0 4、秩:向量组的极在线性无关组所含向量的个数成为这个向量组的秩。 5、矩阵B的秩:行秩,指矩阵的行向量组的秩;列秩类似。记:R(B)

关于matlab矩阵分解

(1) LU分解 A是非奇异的,LU分解总是可以进行的。 [L,U]=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换),矩阵X必须是方阵。 [L,U,P]=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PX=LU。矩阵X必须是方阵。 实现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]'; [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。 实现QR分解后,线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。 例7-3 用QR分解求解例7-1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [Q,R]=qr(A); x=R\(Q\b) 或采用QR分解的第2种格式,命令如下: [Q,R,E]=qr(A); x=E*(R\(Q\b)) (3) Cholesky分解

线性方程组求解matlab实现

3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意:因为RA~=RB ,所以此方程组无解.') return end if RA==RB if RA==n disp('请注意:因为RA=RB=n ,所以此方程组有唯一解.') else disp('请注意:因为RA=RB> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意:因为RA=RB=n ,所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入 >>X=A\b, 运行后输出结果为 X =(0 0 0 0)’. (2) 在MATLAB 工作窗口输入程序 >> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果 请注意:因为RA=RB> A=[4 2 -1;3 -1 2;11 3 0]; b=[2;10;8]; [RA,RB,n]=jiepb(A,B) 运行后输出结果 请注意:因为RA~=RB ,所以此方程组无解. RA =2,RB =3,n =3 (4)在MATLAB 工作窗口输入程序

Matlab_中的矩阵分解函数

Matlab 中的矩阵分解函数 矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解(三角分解)、QR分解(正交变换)、Cholesky分解,以及Schur分解、Hessenberg分解、奇异分解等。 (1) LU分解(三角分解)lu函数 [L,U]=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换),使之满足X=LU。注意,这里的矩阵X必须是方阵。 实现LU分解后,线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb),这样可以大大提高运算速度。 [L,U,P]=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PX=LU。当然矩阵X同样必须是方阵。(设P 是一个m×n的(0,1) 矩阵,如m≤n且P*P′=E,则称P为一个m×n的置换矩阵。) 例1用LU分解求解例1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]';

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。 实现QR分解后,线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。例2用QR分解求解例1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]';

第4讲(4)Matlab中的矩阵分解命令

1 4—6 矩阵分解的Matlab 命令 2 (1)矩阵的LU 分解 (2)矩阵的QR 分解(3 )矩阵的Cholesky 分解(4) 矩阵的奇异值分解(5)矩阵的特征值分解(6) 矩阵的Schur 分解(7) 矩阵的Jordan 标准型分解 3(1)矩阵的LU 分解 在Matlab 中用函数lu 来实现矩阵的LU 分解,其命令格式为:[L, U]=lu(X) 说明:U 矩阵为上三角矩阵,满足X=L*U. 4 [L,U,P]=lu(X) 说明:返回的P 矩阵是置换矩阵,矩阵U 是上三角矩阵,矩阵L 满秩矩阵,满足L*U=P*X. 5例4.1 >>a=[1,2,4;2,-1,5;-1,10,4];>> [b,c,p]=lu(a)运行结果:b = 1.0000 0 0-0.5000 1.0000 00.5000 0.2632 1.0000 6 c = 2.0000 -1.0000 5.00000 9.5000 6.50000 0 -0.2105p =0 1 00 0 11 0 0

7 (2)矩阵的QR 分解 在Matlab 中,矩阵的QR 分解可由函数qr 来实现,其常用的调用格式如下:①[B,C]=qr(A) 说明:返回的C 矩阵为上三角矩阵,矩阵B 为满秩矩阵。 [Q,R,E]=qr(A) 说明:返回的矩阵E 是置换矩阵,矩阵R 是上三角矩阵,矩阵Q 是满秩矩阵,上述矩阵满足关系A*E=Q*R. 8 例4.2 >> a=[1,2,4;2,-1,5;-1,10,4];>> [b,c,e]=qr(a)运行结果:b = -0.1952 -0.5068 -0.83970.0976 -0.8619 0.4976-0.9759 0.0152 0.2177 9 c = -10.2470 -4.1964 0.9759 0 -6.2762 -2.24580 0 -0.0622e =0 0 11 0 00 1 0 10 (3 )矩阵的Cholesky 分解 在Matlab 中用函数chol 对矩阵进行Cholesky 分解,函数chol 的调用格式为:1.R=chol(X) 说明:矩阵X 必须是正定矩阵,否则会返回错误信息,返回的矩阵R 是上三角矩阵。2.[R,p]=chol(X) 说明:此调用格式不管矩阵X 是否正定,都不会返回错误信息。如果矩阵X 正定,则返回上三角矩阵R, p 为零;如果矩阵X 非正定,则返回的矩阵R 也是上三角矩阵,但p 为正数。 11例4.4 >> a=[3,-1,1;-1,5,2;1,2,4];>> b=chol(a)运行结果:b = 1.7321 -0.5774 0.57740 2.1602 1.08010 0 1.5811 12 (4) 矩阵的奇异值分解 在Matlab 中,矩阵的奇异值分解则由函数svd 来实现,其调用格式为:[b,c,d]=svd(A) 说明:返回的矩阵b 为左奇异矩阵,矩阵d 为右奇异矩阵,矩阵c 为奇异值矩阵。例4.5 >> a=[1,2,4;2,-1,5;-1,10,4];>> [b,c,d]=svd(a)

实验一用matlab求解线性方程组

实验1.1 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref (A ) %将矩阵A 化为阶梯形的最简式 null (A ) %求满足AX =0的解空间的一组基,即齐次线性方程组的基 础解系 【例1】 求下列齐次线性方程组的一个基础解系,并写出通解: 我们可以通过两种方法来解: 解法1: >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果: ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵,得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x

取x2,x4为自由未知量,扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系,记 所以齐次方程组的通解为 解法2: clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用, 若去掉r ,是什么结果? 执行后可得结果: B= 1 0 1 0 0 1 0 1 ?? ?=-=-0 04321x x x x ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x +=

MATLAB 矩阵分解算法大全

(1)LU 分解法程序:function x=solvebyLU(A,b) % 该函数利用LU分解法求线性方程组Ax=b的解 flag=isexist(A,b); %调用第一小节中的isexist函数判断方程组解的情况if flag==0 disp('该方程组无解!'); x=[]; return; else r=rank(A); [m,n]=size(A); [L,U,P]=lu(A); y(1)=b(1); if m>1 for i=2:m y(i)=b(i)-L(i,1:i-1)*y(1:i-1)'; end end y=y'; % 解Ux=y得原方程组的一个特解 x0(r)=y(r)/U(r,r); if r>1 for i=r-1:-1:1 x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i); end end x0=x0'; if flag==1 %若方程组有唯一解 x=x0; return; else %若方程组有无穷多解 format rat; Z=null(A,'r'); %求出对应齐次方程组的基础解系 [mZ,nZ]=size(Z); x0(r+1:n)=0; for i=1:nZ t=sym(char([107 48+i])); k(i)=t; %取k=[k1,k2...,]; end x=x0; for i=1:nZ x=x+k(i)*Z(:,i); %将方程组的通解表示为特解加对应齐次通解形式 end end end (2)矩阵的QR分解法(c语言):

void QR(double a[N][N],double q[N][N],double r1[N][N],int n) /*QR分解*/ { int i,j,k,r,m; double temp,sum,dr,cr,hr; double ur[N],pr[N],wr[N]; double q1[N][N],emp[N][N]; for(i=1;i=ZERO) { sum=0; for(k=r;kZERO)m=-1; else m=1; cr=m*dr; hr=cr*(cr-a[r][r]); for(i=1;ir) ur[i]=a[i][r]; }; for(i=1;i

Matlab线性方程组求解(Gauss消去法)

Matlab线性方程组求解 1. Gauss消元法: function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end 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 det=det*a(k,k); %计算行列式 end det=det*a(n,n); 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 Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 2. 列主元Gauss消去法: function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0; %选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n

Matlab求解线性方程组非线性方程组

求解线性方程组 solve,linsolve 例: A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开,元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff(fun,var,n):对表达式fun中的变量var求n阶导数。 例如:F=sym('u(x,y)*v(x,y)'); %sym()用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty():用习惯书写方式显示变量;ans是答案表达式 非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名;fun其中 x0位求解方程的初始向量或矩阵; option为设置命令参数 建立文件fun.m: function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注: ...为续行符 m文件必须以function为文件头,调用符为@;文件名必须与定义的函数名相同;fsolve()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中,求解线性方程组时,主要采用前面章节介绍的除法运算符“/”和“\”。如: X=A\B表示求矩阵方程AX=B的解; 的解。XA=B表示矩阵方程B/A=X. 对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。 如果矩阵A不是方阵,其维数是m×n,则有: m=n 恰定方程,求解精确解; m>n 超定方程,寻求最小二乘解; m

矩阵的LU分解(自编MATLAB)实验报告

LU 分解原理 定理:设A C n n ,如果 A 的顺序主子式 A 11 ≠0, |a 11 a 12 a 21 a 22|≠0,…,|a 11a 12a 21a 22…a 12…a 22??a n?11a n?12? ?a n?1n?1 |≠0 则存在唯一的主对角线上元素全为 1 的下三角矩阵L 与唯一的上三角矩阵 U ,使得 A =LU . 证明:对矩阵A 的阶数使用数学归纳法. 显然,当 n=1 时,A 11=1 ?A 11 就是唯一的分解式。现假定对 n-1 阶矩阵,定理的结论成立。对 A 进行分块 A =( A A ?A A A A A A A AA ) 其中A A ,A A ∈A A ?A .由于 n-1 阶矩阵 A A ?A 的 k 阶顺序主子式就是 A 的 k 阶主子式(k=1,2,…,n-2),故它们都不为零.从而由归纳法假设,A A ?A 有唯一的 LU 分解 A A ?A =A A ?A A A ?A 其中A A ?A 的主对角线上的元素都1.由于 |A A ?A |=| A 11A 12A 21A 22 …A 12…A 22 ? ? A A ?11A A ?12 ? ? A A ?1A ?1 |=|A A ?A A A ?A |≠0 所以A A ?A 及A A ?A 是n-1阶可逆矩阵 先假设已有 A =LU ,其中 L =( A A ?A 0A A 1 ), U= ( A A ?A A A A A AA ) A ,A ∈A A ?A 是待定向量。作乘积 AA = (A A ?A A A ?A A A ?A A A A A A ?A A AA +A A A ) =(A A ?A A A A A A A AA )=A 则A,A 必须满足 A A ?A A =A A ,A A A A ?A =A A A ,A AA +A A A =A AA

线性方程组求解Matlab程序(精.选)

线性方程组求解 1.直接法 Gauss消元法: function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end 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 det=det*a(k,k); end

det=det*a(n,n); 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 Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 列主元Gauss消去法: function x=detGauss(a,b) % Gauss列主元消去法

[n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0;% 选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; 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;det=-det; end

矩阵分解的MATALAB实现

5.3.3 矩阵分解的MATALAB实现 矩阵分解(decomposition, factorization)是多半将矩阵拆解为数个三角形矩阵(triangular matrix),依使用目的的不同,可分为三种矩阵分解法:1)三角分解法(Triangular Factorization),2)QR分解法(QR Factorization),3)奇异值分解法(Singular Value Decompostion)。 (1) 三角分解法 三角分解法是将原正方(square) 矩阵分解成一个上三角形矩阵或是排列(permuted) 的上三角形矩阵和一个下三角形矩阵,这样的分解法又称为LU 分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程,求反矩阵,和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵。 我们举以下二个矩阵为例: 利用三角分解法可将A和B二矩阵分别拆解为上下三角形矩阵 注意B分解的矩阵得到的第一个矩阵[LB]是排列的下三角形矩阵,如果第二、三列互换,则此变成完全的下三角形矩阵。 以MATLAB函数计算上述的LU分解法,其语法为[L,U]=lu(A),其中L代表下三角形矩阵U代表上三角形矩阵。我们来看一个例子。 >> A = [1 2 -1; -2 -5 3; -1 -3 0]; B=[1 3 2; -2 -6 1; 2 5 7]; >> [L1,U1] = lu(A); [L2,U2] = lu(B); >> L1; U1 L1 = % 注意这个矩阵L1和之前的[LA]不相同 -0.5 1 0 1 0 0 0.5 1 1

matlab常用解方程及方程组函数

1. roots 求解多项式的根 r=roots(c) 注意: c 为一维向量,者返回指定多项式的所有根( 包括复根),poly 和roots 是互为反运算,还有就是roots 只能求解多项式的解 还有下面几个函数poly2sym、sym2poly 、eig >>syms x >>y=x A5+3*x A3+3; >>c=sym2poly(y);%求解多项式系数 >>r=roots(c); >>poly(r) 2. residue 求留数 [r, p, k] = residue(b,a) >>b = [ 5 3 -2 7] >>a = [-4 0 8 3] >>[r, p, k] = residue(b,a) 3. solve 符号解方程(组)——使用最多的 g = solve(eq1,eq2,...,eqn,var1,var2,...,varn) 注意:eqn 和varn 可以是符号表达式,也可以是字符串表达式,但是使用符号表达式时不能有“=号”,假如说varn 没有给出,使用findsym 函数找出默认的求解变量。返回的g 是个结构体,以varn 为字段。由于符号求解的局限性,好多情况下可能得到空矩阵,此时只能用数值解法 解方程A=solve('a*xA2 + b*x + c') 解方程组B=solve('a*uA2 + vA2', 'u - v = 1', 'aA2 - 5*a + 6') 4. fzero 数值求零点 [x,fval,exitflag,output]=fzero(fun,x0,options,p1,p2...) fun 是目标函数,可以是句柄(@)、inline 函数或M 文件名 x0 是初值,可以是标量也可以是长度为2 的向量,前者给定一个位置,后者是给定一个范围options 是优化参数,通过optimset 设置,optimget 获取,一般使用默认的就可以了,具体参照帮助 p1,p2...为需要传递的其它参数 假如说(x/1446)A2+p/504.1+(t/330.9)*(log(1-x/1446)+(1-1 /5.3)*x/1446)=0 的根,其中p,t 是已知

用matlab解线性方程组

用matlab解线性方程组 2008-04-12 17:00 一。高斯消去法 1.顺序高斯消去法 直接编写命令文件 a=[] d=[]' [n,n]=size(a); c=n+1 a(:,c)=d; for k=1:n-1 a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去 end x=[0 0 0 0]' %回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1 x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end 2.列主高斯消去法 *由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行,所以要通过修正来把m 改成真正的行的值。该程序只是演示程序,真正机器计算不需要算主元素所在列以下各行应为零的值。 直接编写命令文件 a=[] d=[] ' [n,n]=size(a); c=n+1 a(:,c)=d; %(增广) 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],:); %换行 end a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去end end x=[0 0 0 0]' %回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1 x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end

矩阵的Cholesky分解的Matlab实现

矩阵的Cholesky分解的Matlab实现 下次再补上改进的%Cholesky分解法, %Cholesky分解法 function [X]=m_chol(A,b) [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1); for i=1:N A(i,i)=sqrt(A(i,i)-A(i,1:i-1)*A(i,1:i-1)'); if A(i,i)==0 'A is singular. no unique solution' break end for j=i+1:N A(j,i)=(A(j,i)-A(j,1:i-1)*A(i,1:i-1)')/A(i,i); end end A b %前代法 for j=1:N Y(j)=(b(j)-A(j,1:j-1)*Y(1:j-1))/A(j,j); end Y % A=A' for k=N:-1:1 X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); end Cholesky分解 如果矩阵X是对称正定的,则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R,则下三角矩阵为其转置,即X=R'R。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解,其调用格式为: R=chol(X):产生一个上三角阵R,使R'R=X。若X为非对称正定,则输出一个出错信息。

[R,p]=chol(X):这个命令格式将不输出出错信息。当X为对称正定的,则p=0,R与上述格式得到的结果相同;否则p为一个正整数。如果X 为满秩矩阵,则R为一个阶数为q=p-1的上三角阵,且满足 R'R=X(1:q,1:q)。 实现Cholesky分解后,线性方程组Ax=b变成R‘Rx=b,所以x=R\(R’\b)。例7-4 用Cholesky分解求解例7-1中的线性方程组。 命令如下: A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; R=chol(A) ??? Error using ==> chol Matrix must be positive definite 命令执行时,出现错误信息,说明A为非正定矩阵。 转自: http://203.208.37.104/search?q=cache:EfIWKymoWB8J:www.math.o https://www.sodocs.net/doc/fb11444691.html,/forums/index.php%3Fact%3DAttach%26type%3Dpost%26id%3 D213920+matlab+%E7%9F%A9%E9%98%B5LU%E5%88%86%E8 %A7%A3&hl=zh-CN&ct=clnk&cd=15&gl=cn&client=firefox-a&st_usg =ALhdy2-W-UGiapmEd7-JkiCNACw5NK2Gew

雅可比解线性方程组matlab

雅可比迭代 使用雅可比迭代法求解线性方程组的步骤 步骤1:输入系数矩阵A和方程组右端向量B; 步骤2:将矩阵A分解为下三角阵L对角阵D和上三角阵U 可分解为(D+L+U)X=B for o=1:n d(o,o)=a(o,o); u(o,o+1:n)=-a(o,o+1:n); end for p=2:n l(p,1:p-1)=-a(p,1:p-1); end; 步骤3:将上式化简为x=B0x+f,其中B0=-D-1(L+U),f=D-1B for i=1:n b0(i,i+1:n)=u(i,i+1:n)/a(i,i); f(i,:)=b(i,:)/a(i,i); end for p=2:n b0(p,1:p-1)=l(p,1:p-1)/a(p,p);;

步骤4:采用迭代公式在允许误差范围e=1e-7内求得解向量x x0=x; x=b0*x+f 雅可比迭代法matlab程序: function [x,k]=jacobi(a,b) n=length(a); e=1e-7; m=100; x0=zeros(n,1); x=x0; k=0; d=zeros(n); l=zeros(n); u=zeros(n); b0=zeros(n); f=zeros(n,1);

x0=x+2*e; for o=1:n d(o,o)=a(o,o); u(o,o+1:n)=-a(o,o+1:n); end for p=2:n l(p,1:p-1)=-a(p,1:p-1); end for i=1:n b0(i,i+1:n)=u(i,i+1:n)/a(i,i); f(i,:)=b(i,:)/a(i,i); end for p=2:n b0(p,1:p-1)=l(p,1:p-1)/a(p,p); end while max(abs(x0-x))>e&k

相关主题