线性方程组求解Matlab程序
第三章MATLAB线性方程组及矩阵特征值
![第三章MATLAB线性方程组及矩阵特征值](https://img.taocdn.com/s3/m/9ca18b9985254b35eefdc8d376eeaeaad1f316a9.png)
情形2:m<n(不定方程)
- 0.8000
情形3:m>n(超定方程),多用于曲线拟合。
解线性方程组的一般函数文件如下:
function [x,y]=line_solution(A,b)
[m,n]=size(A);y=[];
if norm(b,1)>0
%非齐次方程组
if rank(A)==rank([A,b]) %方程组相容
for i=1:3
if i~=2, a(i,:)=a(i,:)-a(i,2)*a(2,:); end
end
a
a(3,:)=a(3,:)/a(3,3)
for i=1:3;
if i~=3, a(i,:)=a(i,:)-a(i,3)*a(3,:); end;
end; a A_inv = a(:,4:6) A*A_inv
2in
ai 2
,对调2
r2行.
消元:用a22把ai2消为0 (i 3, 4, , n) :
第2 行
ai 2 a22
第i行,则
aij
a2 j
ai 2 a22
aij (i
3, 4,
, n;j
2, 3,
, n 1)
到此原方程组化为
a11x1 a12 x2 a13 x3
a22 x2 a23 x3
2, 3,
, n;j 1, 2,
, n 1)
到此原方程组化为
a11 x1 a12 x2 a22 x2
ai2 x2
an2 x2
a1n xn a1,n1 a2n xn a2,n1
ain xn ai,n1
ann xn an,n1
(2) 找r2,使 ar2 2
用matlab求解线性方程组
![用matlab求解线性方程组](https://img.taocdn.com/s3/m/48a44e0290c69ec3d5bb7524.png)
用matlab 解线性方程组电子科技大学摘要:利用matlab 软件编写程序,分别利用雅克比迭代法和高斯赛德尔迭代法、列主元高斯消去法,改进平方根法求解不同方程组,其中对于雅克比迭代法和高斯赛德尔迭代法在收敛条件相同的情况下,比较两者的迭代次数,对于,列主元高斯消去法和改进平方根法,要求解出方程组的根。
关键词:雅克比迭代法;高斯赛德尔迭代法;列主元高斯消去法;改进平方根法引言:众所周知,在数学物理方程中,当涉及到解方程组的时候,按照常规的计算方法计算量很大,这样,就涉及到了计算方法的问题,算法里面,很多涉及到矩阵转换,经过处理,可以让我们简便的计算根,而matlab 是一个处理矩阵方程组很便利的软件,下面就是用几种不同的方法解方程组。
正文:一、雅克比迭代法和高斯赛德尔迭代法 1雅可比迭代法原理: 设线性方程组b Ax =的系数矩阵A 可逆且主对角元素nn a ,...,a ,a 2211均不为零,令()nn a ,...,a ,a diag D 2211=并将A 分解成()D D A A +-= 从而(1)可写成 ()b x A D Dx +-= 令 11f x B x +=其中b D f ,A D I B 1111--=-=. 以1B 为迭代矩阵的迭代法(公式)()()111f x B x k k +=+称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,则为⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+其中()()()()()Tn x ,...x ,x x 002010=为初始向量.由此看出,雅可比迭代法公式简单,每迭代一次只需计算一次矩阵和向量的乘法.在电算时需要两组存储单元,以存放()k x 及()1+k x . 2高斯赛德尔迭代法原理由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经算出最新的分量,但没被利用。
MATLAB程序设计第六讲
![MATLAB程序设计第六讲](https://img.taocdn.com/s3/m/22281c0416fc700abb68fc1c.png)
MATLAB程序设计杨凯2010 . 11主要内容自学))*MATLAB解方程与函数极值解方程与函数极值((自学(自学)自学)线性方程组求解(一、线性方程组求解二、非线性方程组求解三、函数极值四、常微分方程初值问题的数值解法*MATLAB符号计算一、符号计算基础二、微积分三、简化方程表达式四、解方程一、线性方程组求解(自学)1.1 直接解法1.利用左除运算符的直接解法对于线性方程组Ax=b,可以利用左除运算符“\”求可以利用左除运算符“解:x=A\b例*:用直接解法求解下列线性方程组用直接解法求解下列线性方程组。
命令如下命令如下::A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b2.利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积算法将一个矩阵分解成若干个矩阵的乘积。
常见的矩阵分解有LU 分解分解、、QR 分解分解、、Cholesky 分解分解,,以及Schur 分解分解、、Hessenberg 分解分解、、奇异分解等奇异分解等。
(1) LU 分解矩阵的LU 分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式角矩阵和一个上三角矩阵的乘积形式。
线性代数中已经证明,只要方阵A 是非奇异的是非奇异的,,LU 分解总是可以进行的分解总是可以进行的。
MATLAB 提供的lu 函数用于对矩阵进行LU 分解分解,,其调用格式为式为::[L,U]=lu(X):产生一个上三角阵U 和一个变换形式的下三角阵L(行交换),使之满足X=LU 。
注意注意,,这里的矩阵X 必须是方阵是方阵。
[L,U,P]=lu(X):产生一个上三角阵U 和一个下三角阵L 以及一个置换矩阵P ,使之满足PX=LU 。
当然矩阵X 同样必须是方阵方阵。
实现LU 分解后分解后,,线性方程组Ax=b 的解x=U\(L\b)或x=U\(L\Pb),这样可以大大提高运算速度这样可以大大提高运算速度。
MatLab解线性方程组
![MatLab解线性方程组](https://img.taocdn.com/s3/m/bfb6e664b84ae45c3b358ca9.png)
MatLab解线性方程组一文通(转帖)当齐次线性方程AX=0,rank(A)=r<n时,该方程有无穷多个解,怎样用MATLAB求它的一个基本解呢?用matlab 中的命令x=null(A, r )即可.其中:r=rank(A)A=[ 1 1 1 1 -3 -1 11 0 0 0 1 1 0-2 0 0 -1 0 -1 -2]用matlab 求解程序为:A=[1 1 1 1 -3 -1 1;1 0 0 0 1 1 0;-2 0 0 -1 0 -1 -2];r=rank(A);y=null(A, r )得到解为:y=[ 0 -1 -1 0-1 2 1 11 0 0 00 2 1 -20 1 0 00 0 1 00 0 0 1]其列向量为Ay=0的一个基本解MatLab解线性方程组一文通!-------------------作者:liguoy(2005-2-3)写在阅读本文前的引子。
一:读者对线性代数与Matlab 要有基本的了解;二:文中的通用exp.m文件,你须把具体的A和b代进去。
一:基本概念1.N级行列式A:|A|等于所有取自不同行不同列的n个元素的积的代数和。
2.矩阵B:矩阵的概念是很直观的,可以说是一张表。
3.线性无关:一向量组(a ,a ,…. a )不线性相关,即没有不全为零的数k ,k ,……kn使得:k1* a +k2* a +…..+kn*an=04. 秩:向量组的极在线性无关组所含向量的个数称为这个向量组的秩。
5.矩阵B的秩:行秩,指矩阵的行向量组的秩;列秩类似。
记:R(B)6.一般线性方程组是指形式: (1)其中x1,x2,…….xn为n个未知数,s为方程个数。
记:A*X=b7.性方程组的增广矩阵:=8. A*X=0 (2)二:基本理论三种基本变换:1,用一非零的数乘某一方程;2,把一个方程的倍数加到另一个方程;3互换两个方程的位置。
以上称初等变换。
消元法(理论上分析解的情况,一切矩阵计算的基础)首先用初等变换化线性方程组为阶梯形方程组,把最后的一些恒等式”0=0”(如果出现的话)去掉,1:如果剩下的方程当中最后的一个等式是零等于一非零数,那么方程组无解;否则有解,在有解的情况下,2:如果阶梯形方程组中方程的个数r等于未知量的个数,那么方程组有唯一的解,3:如果阶梯形方程组中方程的个数r小于是未知量的个数,那么方程组就有无穷个解。
专题4 使用MATLAB求解线性方程组的不同方法
![专题4 使用MATLAB求解线性方程组的不同方法](https://img.taocdn.com/s3/m/d439c56e0b4c2e3f5727637b.png)
Z = null(A) 求出 Ax=0 的基础解系后,将基础解系的向量正交单位化,存储在 Z 中. MATLAB 源代码如下: A=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3]
A= 12 2 1 2 1 -2 -2 1 -1 -4 -3
Rank(A) Ans= 2 A=sym(A) A= [1,2,2,1] [2,1,-2,-2] [1,-1,-4,-3] null(A) ans= [2, 5/3] [-2,-4/3] [1, 0] [0, 1]
运行结果为: rank_A =
2 rank_B =
2 S_H =
-2
1
1
0
0
2
0
1
S_P =
0
1.7500
0
-0.5000
则该线性方程组有无穷多解为:
2 1 0
x
k1
1 0
k2
0 2
7 0
/
4
,
k1
,Leabharlann k2R 0 1 1/ 2
nulla?r?求系数矩阵为a的齐次线性方程组ax0的基础解系结果为有理数bnulla求出ax0的基础解系后将基础解系的向量正交单位化存储在zmatlab源代码如下
专题 4 使用 MATLAB 求解线性方程组的方法
x1 2 x2 2x3 x4 0
【例
1】求齐次方程组
2 x1
end end
x1 2x2 2x3 3x4 2 【例 1.3】使用 Matlab 求解方程组 2x1 4x2 3x3 4x4 5
matlab线性方程组的矩阵解法
![matlab线性方程组的矩阵解法](https://img.taocdn.com/s3/m/b0eba5c65fbfc77da269b102.png)
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分解。充分利用 了系数矩阵的特点,而且使之分解更简单, 了系数矩阵的特点,而且使之分解更简单,得到对三对 角线性方程组的快速解法。 角线性方程组的快速解法。
matlab矩阵分解与线性方程组求解
![matlab矩阵分解与线性方程组求解](https://img.taocdn.com/s3/m/41dd73b6fd0a79563c1e7229.png)
格式
[Q, R] = rsf2csf(q, r) 例4-7
A=[1 1 1 3;1 2 1 1;1 1 3 1;-2 1 1 4]; [q, r]=schur (A) [Q, R]=rsf2csf(q, r)
4.2 秩与线性相关性
4.2.1
汪远征
矩阵和向量组的秩与向量组的线性相关性
矩阵 A 的秩是指矩阵 A 中最高阶非零子式的阶数,或
是矩阵线性无关的行数与列数;向量组的秩通常由该
向量组构成的矩阵来计算。 k = rank(A) 返回矩阵A的行(或列)向量中线性无关个数 k = rank(A,tol) tol为给定误差
在 MATLAB 中,求矩阵秩的函数是 rank 。其格式为:
4.2 秩与线性相关性
4.2.1
汪远征
矩阵和向量组的秩与向量组的线性相关性
4.2 秩与线性相关性
4.2.2
汪远征
求行阶梯矩阵及向量组的基
Matlab 将矩阵化成行最简形的命令是 rref或 rrefmovie 。
其格式为:
R = rref(A) R 是A的行最简行矩阵 [R,jb] = rref(A) jb 是一个向量,其含义为: r = length(jb) 为 A 的秩; A(:, jb)为A的列向量基;jb中元素表示基向量所在的 列。
阵。
4.1 矩阵分解
4.1.2
汪远征
Cholesky分解
例4-2
A=pascal(4) %产生4阶pascal矩阵 [R,p]=chol(A)
4.1 矩阵分解
4.1.3
汪远征QBiblioteka 分解将矩阵A分解成一个正交矩阵Q与一个上三角矩阵R的
乘积A=QR,称为QR分解。
MATLAB-平方根法和改进平方根法求解线性方程组例题与程序演示教学
![MATLAB-平方根法和改进平方根法求解线性方程组例题与程序演示教学](https://img.taocdn.com/s3/m/ce85aab1f705cc17552709a7.png)
M A T L A B-平方根法和改进平方根法求解线性方程组例题与程序(2)设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945x x x x x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎥----⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢---⎢⎥⎢⎥⎢--⎢⎥⎢⎢⎥⎣⎦⎣⎦⎣⎦⎥⎥⎥⎥ (1,1,0,2,1,1,0,2)T x *=--二、数学原理 1、平方根法解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法,这里对系数矩阵A 是有要求的,需要A 是对称正定矩阵,根据数值分析的相关理论,如果A 对称正定,那么系数矩阵就可以被分解为的T A=L L •形式,其中L 是下三角矩阵,将其代入Ax=b 中,可得:T LL x=b 进行如下分解:T L xL by y ⎧=⎨=⎩ 那么就可先计算y,再计算x ,由于L 是下三角矩阵,是T L 上三角矩阵,这样的计算比直接使用A 计算简便,同时你应该也发现了工作量就转移到了矩阵的分解上面,那么对于对称正定矩阵A 进行Cholesky 分解,我再描述一下过程吧: 如果你对原理很清楚那么这一段可以直接跳过的。
设T A=L L •,即1112111112112122221222221212....................................n n n n n n nn n n nn nn a a a l l l l aa a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步,由矩阵乘法,211111111,i i a l a l l ==g 故求得111111,2,3,...i i a l l i n a === 一般的,设矩阵L 的前k-1列元素已经求出 第k 步,由矩阵乘法得112211k k kk kmkkik im km ik kkm m a l l a l l l l --===+=+∑∑, 于是11(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩∑ 2、改进平方根法在平方根的基础上,为了避免开方运算,所以用TLDL A =计算;其中,11122.........n d D D D d ⎤⎤⎡⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎢⎥⎣⎦⎣⎣;得1121212212111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦L L M MO O O M L按行计算的L 元素及D 对元素公式 对于n i ,,2,1Λ=11(1,21)j ij ij ik jk k t a t l j i -==-=-∑…,./(1,2,)ij ij j l t d j ==…,i-1.11i i ii ik ikk d a t l -==-∑计算出LD T =的第i 行元素(1,2,i-1)ij t j =…,后,存放在A 的第i 行相置,然后再计算L 的第i 行元素,存放在A 的第i 行.D 的对角元素存放在A 的相应位置.对称正定矩阵A 按T LDL 分解和按T LL 分解计算量差不多,但T LDL 分解不需要开放计算。
利用matlab解线性方程组
![利用matlab解线性方程组](https://img.taocdn.com/s3/m/13f8ba13c281e53a5802ffbd.png)
数值计算实验——解线性方程组西南交通大学2012级茅7班20123257 陈鼎摘要本报告主要介绍了基于求解线性方程组的高斯消元法和列主消元法两种数值分析方法的算法原理及实现方法。
运用matlab数学软件辅助求解。
实验内容1.编写用高斯消元法解线性方程组的MATLAB程序,并求解下面的线性方程组,然后用逆矩阵解方程组的方法验证。
2.编写用列主消元法解线性方程组的MATLAB程序,并求解下面的线性方程组,然后用逆矩阵解方程组的方法验证。
给定方程组如下:①0.325x1+2.564x2+3.888x3+5x4=1.521②-1.548x1+3.648x2+4.214x3-4.214x4=2.614③-2.154x1+1.647x2+5.364x3+x4=3.978④0x1+2.141x2-2.354x3-2x4=4.214A.高斯消元法一、算法介绍高斯消元法是一种规则化的加减消元法。
基本思想是通过逐次消元计算把需要求解的线性方程组转化成为上三角方程组,即把现形方程组的系数矩阵转化为上三角矩阵,从而使一般线性方程组的求解转化为等价的上三角方程组的求解。
二、matlab程序function [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和常系数向量b;输出的量:系数矩阵A和增广矩阵B的秩RA、RB,方程中未知量的个数n和有关方程组解X及其解的信息。
matlab超松弛迭代法求方程组
![matlab超松弛迭代法求方程组](https://img.taocdn.com/s3/m/1e4f485c11a6f524ccbff121dd36a32d7375c7b8.png)
一、介绍MATLAB(Matrix Laboratory)是一种用于数值计算和数据可视化的专业软件。
在MATLAB中,超松弛迭代法是解决线性方程组的一种有效算法。
本文将介绍MATLAB中超松弛迭代法的基本原理和实现方法,并给出一个具体的例子进行演示。
二、超松弛迭代法的基本原理超松弛迭代法是一种逐步迭代的算法,用于求解线性方程组。
它的基本原理是通过不断迭代更新方程组的解,直到达到满足精度要求的解。
超松弛迭代法的公式如下:X(k+1) = (1-w)X(k) + w*(D-L)⁻¹*(b+U*X(k))其中,X(k)代表第k次迭代的解向量,X(k+1)代表第k+1次迭代的解向量,D、L和U分别代表方程组的对角线元素、下三角元素和上三角元素构成的矩阵,b代表方程组的右端向量,w代表松弛因子。
超松弛迭代法的关键在于选择合适的松弛因子w,一般情况下,可以通过试验选取一个合适的值。
在MATLAB中,可以使用sor函数来实现超松弛迭代法。
三、MATLAB中超松弛迭代法的实现方法在MATLAB中,可以通过调用sor函数来实现超松弛迭代法。
sor 函数的语法格式如下:[X,flag,relres,iter,resvec] = sor(A,b,w,tol,maxit)其中,A代表线性方程组的系数矩阵,b代表右端向量,w代表松弛因子,tol代表迭代的精度要求,maxit代表最大迭代次数,X代表迭代求解得到的解向量,flag代表迭代的结果标志,relres代表相对残差的大小,iter代表迭代次数,resvec代表迭代过程中的残差向量。
以下是一个使用sor函数求解线性方程组的示例:A = [4 -1 0 -1 0 0; -1 4 -1 0 -1 0; 0 -1 4 0 0 -1; -1 0 0 4 -1 0; 0 -1 0 -1 4 -1; 0 0 -1 0 -1 4];b = [1; 0; -1; 0; 1; 0];w = 1.25;tol = 1e-6;maxit = 100;[X,flag,relres,iter,resvec] = sor(A,b,w,tol,maxit);通过调用sor函数,可以得到方程组的解向量X,迭代的结果标志flag,相对残余resrel和迭代次数iter。
MATLAB求解非齐次线性方程组
![MATLAB求解非齐次线性方程组](https://img.taocdn.com/s3/m/8091244d326c1eb91a37f111f18583d049640fa3.png)
MATLAB 求解⾮齐次线性⽅程组根据线性代数中求解⽅程组的基本知识,⾸先应判断系数矩阵的秩是否和增⼴矩阵的秩相等,若不等,则⽆解;若有解,根据秩和未知量个数的关系,判断是唯⼀解还是⽆穷多解;若为⽆穷多解,其通解为齐次⽅程组的通解加⾮齐次⽅程组的特解。
求⾮齐次线性⽅程组Ax=b 的特解,可直接使⽤命令A\b ,求解齐次线性⽅程组的通解,可以使⽤函数null 或rref 来实现。
命令含义B = null(A,'r')求系数矩阵为A 的齐次线性⽅程组Ax=0的基础解系,结果为有理数,B 的列向量即基础解系的列向量Z = null(A)求出Ax=0的基础解系后,将基础解系的向量正交单位化,存储在Z 中C = rref(A)求出矩阵A 的⾏最简形矩阵(reduced row echelon form )function [S_H, S_P] = solveLS(A,b)% 输⼊参数A :系数矩阵% 输⼊参数b :Ax=b 的常数项列向量b% S_H :齐次线性⽅程组的基础解系% S_P :⾮齐次线性⽅程组的特解if size(A,1) ~= length(b) %size(A,1)求矩阵的⾏数error('输⼊数据错误,请重新输⼊!');return;elseB = [A,b]; %增⼴矩阵rank_A = rank(A); %求系数矩阵的秩rank_B = rank(B); %求增⼴矩阵的秩if rank_A ~= rank_B %⽆解情况disp('线性⽅程组⽆解!');S_H = [];S_P = [];else if rank_B == size(A,2) %若增⼴矩阵的秩 = 未知量个数%size(A,2)求矩阵的列数,相当于length(A)disp('线性⽅程组有唯⼀解!');S_P = A\b; %求唯⼀解S_H = [];elsedisp('线性⽅程组有⽆穷解!');S_H = null(A,'r');%求出齐次⽅程组的基础解系S_P = A\b; %求⾮齐次⽅程组的特解endendend例 使⽤Matlab 求解⽅程组A=[1 2 -2 3; 2 4 -3 4; 5 10 -8 11];b=[2 5 12]';format rat;[S_H, S_P]=solveLS(A,b)运⾏结果线性⽅程组有⽆穷解!S_H =-2 11 00 20 1S_P =7/4-1/2该线性⽅程组有⽆穷多解,通解为⎧⎩⎨+2−2+3=2x 1x 2x 3x 42+4−3+4=5x 1x 2x 3x 45+10−8+11=12x 1x 2x 3x 4x =++,,∈R k 1⎛⎝⎜⎜⎜−2100⎞⎠⎟⎟⎟k 2⎛⎝⎜⎜⎜1021⎞⎠⎟⎟⎟⎛⎝⎜⎜⎜07/40−1/2⎞⎠⎟⎟⎟k 1k 2。
线性方程组及MATLAB应用
![线性方程组及MATLAB应用](https://img.taocdn.com/s3/m/27cb7f3531b765ce04081423.png)
数值实验 线性方程组与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分解法然后调用求解。
MATLAB编程与工程应用——第7章 MATLAB解方程与函数极值
![MATLAB编程与工程应用——第7章 MATLAB解方程与函数极值](https://img.taocdn.com/s3/m/a7dfec2fb4daa58da0114af9.png)
MATLAB解方程与函数极值
7.1 线性方程组求解
二、迭代解法
迭代解法非常适合求解大型系数矩阵的方程组。在数值 分析中,迭代解法主要包括 Jacobi迭代法、GaussSerdel迭代法、超松弛迭代法和两步迭代法。 Jacobi迭代法 1.Jacobi迭代法 对于线性方程组Ax=b,如果A为非奇异方阵,即aii≠0, 则A可分解为A=D-L-U,其中D为对角阵,其元素为A的对 角元素,L与U为A的下三角阵和上三角阵,于是Ax=b化为: x=D-1(L+U)x+D-1b 与之对应的迭代公式为: x(k+1)=D-1(L+U)x(k)+D-1b 利用jacobi jacobi迭代法求方程的解 例7.4 利用jacobi迭代法求方程的解 exp7_4.m jacobi.m
如s=2时
j =1 j =1
i 1
i 1
1 ym +1 = ym + h[ f ( xm , ym ) + f ( xm + h, ym + h)] 2
matlab求线性方程组的解
![matlab求线性方程组的解](https://img.taocdn.com/s3/m/42fe13cc541810a6f524ccbff121dd36a22dc447.png)
matlab求线性方程组的解求解线性方程分为两种方法–直接法和迭代法常见的方法一共有8种直接法Gauss消去法Cholesky分解法迭代法Jacobi迭代法Gauss-Seidel迭代法超松弛迭代法共轭梯度法Bicg迭代法Bicgstab迭代法这里我从计算代码的角度来解释一下,代码按以下顺序给出。
把方程组直接带入已知条件,就可以得到答案。
适用条件Gauss消去法:求解中小规模线性方程(阶数不过1000),一般用于求系数矩阵稠密而且没有任何特殊结构的线性方程组Cholesky分解法:对称正定方程优先使用,系数矩阵A是n 阶对称正定矩阵Jacobi迭代法非奇异线性方程组,分量的计算顺序没有关系Gauss-Seidel迭代法与Jacobi迭代法相似,但计算的分量不能改变超松弛迭代法Jacobi迭代法和Gauss-Seidel迭代法的加速版,由Gauss-Seidel迭代法改进而来,速度较快共轭梯度法需要确定松弛参数w,只有系数矩阵具有较好的性质时才可以找到最佳松弛因子。
但好处是不用确定任何参数,他是对称正定线性方程组的方法也是求解大型稀疏线性方程组最热门的方法Bicg迭代法本质是用双共轭梯度求解线性方程组的方法,对求解的方程没有正定性要求Bicgstab迭代法本质是用稳定双共轭梯度求解线性方程组的方法,对求解的方程没有正定性要求Gauss消去法第一、二个函数ltri、utri是一定要掌握的,后面的几乎每个函数都要用到ltri简单来说,当Ly=bb,L(非奇异下三角矩阵)已知求yfunction y =ltri(L,b)n=size(b,1);y=zeros(n,1);for j =1:n-1y(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-y(j)*L(j+1:n,j); endy(n)=b(n)/L(n,n);utri简单来说,当Ux=yy,U(非奇异上三角矩阵)已知求xfunction x =utri(U,y)n=size(y,1);x=zeros(n,1);for j = n:-1:2x(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-x(j)*U(1:j-1,j);endx(1)=y(1)/U(1,1);gauss算法,计算时粘贴过去就好function[L,U]=gauss(A)n=size(A,1);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,-1)+eye(n);U=triu(A);使用例子已经知道一个线性方程组,这里我就不写出数学形式了,A是系数矩阵,直接把上面写好的函数复制过来在运算就可以。
MATLAB实验一 解线性方程组的直接法
![MATLAB实验一 解线性方程组的直接法](https://img.taocdn.com/s3/m/f1104ae69e31433239689374.png)
输出 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求解线性方程组、非线性方程组](https://img.taocdn.com/s3/m/84540e8eb8f67c1cfbd6b820.png)
Matlab求解线性方程组、非线性方程组姓名:罗宝晶学号:1012208015 专业:材料学院高分子系第一部分数值计算Matlab求解线性方程组AX=B或XA=B在MATLAB中,求解线性方程组时,主要采用除法运算符“/”和“\”。
如:X=A\B表示求矩阵方程AX=B的解;X=B/A表示矩阵方程XA=B的解。
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。
如果矩阵A不是方阵,其维数是m×n,则有:m=n 恰定方程,求解精确解;m>n 超定方程,寻求最小二乘解;m<n 不定方程,寻求基本解,其中至多有m个非零元素。
针对不同的情况,MATLAB将采用不同的算法来求解。
一.恰定方程组恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式:Ax=b 其中A是方阵,b是一个列向量;在线性代数中,最常用的方程组解法有:(1)利用Cramer公式来求解法;(2)利用矩阵求逆解法,即x=A-1b;(3)利用Gaussian消去法;(4)利用Lu法求解。
一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。
前三种解法的真正意义是在其理论上,而不是实际的数值计算。
MATLAB 中,出于对算法稳定性的考虑,行列式及逆的计算大都在Lu分解的基础上进行。
在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式:x=A\b。
在MATLAB的指令解释器在确认变量A非奇异后,就对它进行Lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。
如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。
此外还要注意:在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。
matlab联立方程
![matlab联立方程](https://img.taocdn.com/s3/m/14a04178b207e87101f69e3143323968011cf4ee.png)
matlab联立方程Matlab是一种常用的科学计算软件,其强大的矩阵计算功能可以用来求解线性方程组。
本文将以Matlab联立方程为主题,介绍如何使用Matlab求解线性方程组,并通过实例来展示其应用。
一、Matlab联立方程的基本原理在数学中,线性方程组是由一组线性方程组成的方程系统。
解线性方程组即找到一组满足所有方程的未知数的取值。
Matlab可以通过矩阵运算的方式来求解线性方程组。
二、Matlab求解线性方程组的步骤1. 定义系数矩阵A和常数向量b:A是由方程组的系数组成的矩阵,b是由方程组的常数项组成的向量。
2. 使用Matlab的左除运算符\求解方程组:x = A\b其中,x是方程组的解向量。
Matlab会自动使用高效的高斯消元法或LU分解法来求解方程组。
三、示例:使用Matlab求解线性方程组假设有以下线性方程组:2x + 3y - z = 73x - 2y + 4z = -2x + y + z = 41. 定义系数矩阵A和常数向量b:A = [2, 3, -1; 3, -2, 4; 1, 1, 1]b = [7; -2; 4]2. 使用Matlab求解方程组:x = A\b运行以上代码,得到方程组的解向量x:x = [3; -1; 2]即方程组的解为x=3,y=-1,z=2。
四、Matlab联立方程的应用Matlab求解线性方程组在工程、物理、经济等领域有着广泛的应用。
例如,在电路分析中,可以利用线性方程组求解电路中的电流和电压;在经济学中,可以利用线性方程组分析供求关系和市场均衡等问题。
五、总结本文介绍了Matlab联立方程的基本原理和求解步骤,并通过实例展示了如何使用Matlab求解线性方程组。
Matlab的强大矩阵计算功能使得求解线性方程组变得简单高效,广泛应用于各个领域。
使用Matlab求解线性方程组可以帮助我们解决实际问题,提高计算效率,减少人工计算的错误。
希望本文对读者能够有所帮助。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
线性方程组求解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-1for i=k+1:nif a(k,k)==0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1 %回代for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample:>> 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.00471.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-1amax=0;% 选主元for i=k:nif abs(a(i,k))>amaxamax=abs(a(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n %进行消元m=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1 %回代for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample:>> x=detGauss(A,b)x =0.9739-0.00471.0010Gauss-Jorda n 消去法:function x=GaussJacobi(a,b)% Gauss-Jacobi 消去法[n,m]=size(a);nb=length(b);x=zeros(n,1);for k=1:namax=0;% 选主元for i=k:nif abs(a(i,k))>amaxamax=abs(a(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%进行消元b(k)=b(k)/a(k,k);for j=k+1:na(k,j)=a(k,j)/a(k,k);endfor i=1:nif i~=kfor j=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendendfor i=1:nx(i)=b(i);endExample:>> x=GaussJacobi(A,b)x =0.9739-0.00471.0010LU 分解法:function [l,u]=lu(a) %LU 分解n=length(a);l=eye(n);u=zeros(n);for i=1:nu(1,i)=a(1,i);endfor i=2:nl(i,1)=a(i,1)/u(1,1);end for r=2:n%%%%for i=r:nuu=0;for k=1:r-1uu=uu+l(r,k)*u(k,i);endu(r,i)=a(r,i)-uu;end%%%%for i=r+1:nll=0;for k=1:r-1ll=ll+l(i,k)*u(k,r);endl(i,r)=(a(i,r)-ll)/u(r,r);end%%%% function x=lusolv(a,b)End%LU 分解求解线性方程组aX=bif length(a)~=length(b) error('Error in inputing!') return;endn=length(a);[l,u]=lu(a);y(1)=b(1);for i=2:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/u(n,n);for i=n-1:-1:1z=0;for k=i+1:nz=z+u(i,k)*x(k);endx(i)=(y(i)-z)/u(i,i);endExample:>> x=lusolv(A,b)x =0.9739 -0.0047 1.0010对称正定矩阵之Cholesky 分解法:function L=Cholesky(A)%对对称正定矩阵 A 进行Cholesky 分解n=length(A); L=zeros(n);for k=1:ndelta=A(k,k);for j=1:k-1delta=delta-L(k,j)A2;endif delta<1e-10return;L(k,k)=sqrt(delta);endfor i=k+1:nL(i,k)=A(i,k);for j=1:k-1L(i,k)=L(i,k)-L(i,j)*L(k,j);endL(i,k)=L(i,k)/L(k,k);endendfunction x=Chol_Solve(A,b)%利用对称正定矩阵之Cholesky 分解求解线性方程组n=length(b);l=Cholesky(A);x=ones(1,n);y=ones(1,n);for i=1:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=(b(i)-z)/l(i,i);for i=n:-1:1Ax=bendz=0;for k=i+1:nz=z+l(k,i)*x(k);endx(i)=(y(i)-z)/l(i,i);endExample:>> a=[9 -36 30 ;-36 192 -180;30 -180 180]; >> b=[1 1 1]';>> x=Chol_Solve(a,b)x =1.8333 1.0833 0.7833对称正定矩阵之LDL'分解法:function [L,D]=LDL_Factor(A)%对称正定矩阵 A 进行LDL' 分解n=length(A);L=eye(n);D=zeros(n);d=zeros(1,n);for k=1:nd(k)=A(k,k);for j=1:k-1d(k)=d(k)-L(k,j)*T(k,j);endif abs(d(k))<1e-10return;endfor i=k+1:nT(i,k)=A(i,k);for j=1:k-1T(i,k)=T(i,k)-T(i,j)*L(k,j);endL(i,k)=T(i,k)/d(k);endendD=diag(d);function x=LDL_Solve(A,b)Ax=b %利用对对称正定矩阵 A 进行LDL' 分解法求解线性方程组n=length(b);[l,d]=LDL_Factor(A);for i=2:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/d(n,n);for i=n-1:-1:1z=0;for k=i+1:nz=z+l(k,i)*x(k);endx(i)=y(i)/d(i,i)-z; endExample:>> x=LDL_Solve(a,b)1.8333 1.0833 0.78332.迭代法Richardson 迭代法:function [x,n]=richason(A,b,x0,eps,M) %Richardson 法求解线性方程组Ax=b %方程组系数矩阵:A%方程组之常数向量:b%迭代初始向量:X0%e 解的精度控制:eps%迭代步数控制:M%返回值线性方程组的解:x%返回值迭代步数:nif(nargin == 3)eps = 1.0e-6;M = 200;elseif(nargin == 4)M = 200;endI =eye(size(A));x1=x0;x=(I-A)*x0+b;n=1;while(norm(x-x1)>eps)x1=x;x=(I-A)*x1+b;n = n + 1;if(n>=M)disp('Warning: 迭代次数太多,现在退出!');return;endendExample:>> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]';x0=[0 0 0]';>> [x,n]=richason(A,b,x0)x =0.9739-0.00471.0010Jacobi 迭代法:function [x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3eps= 1.0e-6;M = 200;elseif nargin<3errorreturnelseif nargin ==5M = varargin{1};endD=diag(diag(A)); %求 A 的对角矩阵L=-tril(A,-%求 A 的下三角阵1);0.9739U=-triu(A,1); %求 A 的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多,可能不收敛!');return;endendExample:>> [x,n]=Jacobi(A,b,x0)x =0.9739-0.00471.0010Gauss-Seidel 迭代法:function [x,n]=gauseidel(A,b,x0,eps,M) if nargin==3 eps= 1.0e-6;M = 200;elseif nargin == 4M = 200;elseif nargin<3errorreturn;endD=diag(diag(A)); %求 A 的对角矩阵L=-tril(A,-1); %求 A 的下三角阵U=-triu(A,1); %求 A 的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1; % 迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多,可能不收敛!');return;endendExample:>> [x,n]=gauseidel(A,b,x0)x =0.97394超松驰迭代法:function [x,n]=SOR(A,b,x0,w,eps,M)if nargin==4-0.00471.0010D=diag(diag(A)); %求 A 的对角矩阵>> [x,n]=SOR(A,b,x0,1)x =eps= 1.0e-6;M = 200; elseif nargin<4errorreturnelseif nargin ==5M = 200;end if(w<=0 || w>=2)error;return;L=-tril(A,-1); %求 A 的下三角阵U=-triu(A,1); %求 A 的上三角阵B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1; % 迭代次数while norm(x-x0)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=M) disp('Warning: 迭代次数太多,可能不收敛!');return;endendExample:0.9739-0.00471.0010 n =4对称逐次超松驰迭代法:function [x,n]=SSOR(A,b,x0,w,eps,M) if nargin==4eps= 1.0e-6;M = 200;elseif nargin<4errorreturnelseif nargin ==5M = 200;endif(w<=0 || w>=2)error;return;endD=diag(diag(A)); %求 A 的对角矩阵L=-tril(A,-1); %求 A 的下三角阵U=-triu(A,1); %求 A 的上三角阵B1=inv(D-L*w)*((1-w)*D+w*U);B2=inv(D-U*w)*((1-w)*D+w*L);f1=w*inv((D-L*w))*b;f2=w*inv((D-U*w))*b;x12=B1*x0+f1;x =B2*x12+f2;n=1; % 迭代次数while norm(x-x0)>=epsx0=x; x12=B1*x0+f1;x =B2*x12+f2;n=n+1;if(n>=M)disp('Warning: 迭代次数太多,可能不收敛!');return;endendExample:>> [x,n]=SSOR(A,b,x0,1)x =0.9739-0.00471.0010n =3两步迭代法:function [x,n]=twostep(A,b,x0,eps,varargin)if nargin==3eps= 1.0e-6;M = 200;elseif nargin<3errorreturnelseif nargin ==5M = varargin{1};endD=diag(diag(A)); %求 A 的对角矩阵L=-tril(A,-1); %求 A 的下三角阵U=-triu(A,1); %求 A 的上三角阵B1=(D-L)\U;B2=(D-U)\L;f1=(D-L)\b;f2=(D-U)\b;x12=B1*x0+f1;x =B2*x12+f2;%迭代次数n=1;while norm(x-x0)>=epsx0 =x;x12=B1*x0+f1;x =B2*x12+f2;n=n+1;if(n>=M)disp('Warning: 迭代次数太多,可能不收敛!');return;endendExample:>> [x,n]=twostep(A,b,x0)x =0.9739-0.00471.0010最速下降法:function [x,n]=fastdown(A,b,x0,eps)if(nargin == 3)eps = 1.0e-6;n =end r = b-A*x0;d = dot(r,r)/dot(A*r,r);x = x0+d*r;n=1;while(norm(x-x0)>eps)x0 = x;r = b-A*x0;d = dot(r,r)/dot(A*r,r);x = x0+d*r;n = n + 1;endExample: >> [x,n]=fastdown(A,b,x0) x =0.9739-0.00471.0010共轭梯度法:function [x,n]=conjgrad(A,b,x0)if(nargin == 3)eps = 1.0e-6;end r1 = b-A*x0;p1 = r1;d = dot(r1,r1)/dot(p1,A*p1);x = x0+d*p1;r2 = r1-d*A*p1;f = dot(r2,r2)/dot(r1,r1);p2 = r2+f*p1;n = 1;for(i=1:(rank(A)-1))x0 = x;p1 = p2;r1 = r2;d = dot(r1,r1)/dot(p1,A*p1);x = x0+d*p1;r2 = r1-d*A*p1;f = dot(r2,r2)/dot(r1,r1);p2 = r2+f*p1;n = n + 1;end d = dot(r2,r2)/dot(p2,A*p2);x = x+d*p2;n = n + 1;Example: >> [x,n]=conjgrad(A,b,x0)0.9739-0.00471.0010 n =4预处理的共轭梯度法:当AX=B 为病态方程组时,共轭梯度法收敛很慢。