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

合集下载

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

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

MATLAB计算方法3解线性方程组计算解法线性方程组是数学中的一个重要问题,解线性方程组是计算数学中的一个基本计算,有着广泛的应用。

MATLAB是一种功能强大的数学软件,提供了多种解线性方程组的计算方法。

本文将介绍MATLAB中的三种解线性方程组的计算方法。

第一种方法是用MATLAB函数“linsolve”解线性方程组。

该函数使用高斯消元法和LU分解法求解线性方程组,可以处理单个方程组以及多个方程组的情况。

使用该函数的语法如下:X = linsolve(A, B)其中A是系数矩阵,B是常数向量,X是解向量。

该函数会根据A的形式自动选择求解方法,返回解向量X。

下面是一个使用“linsolve”函数解线性方程组的例子:A=[12;34];B=[5;6];X = linsolve(A, B);上述代码中,A是一个2×2的系数矩阵,B是一个2×1的常数向量,X是一个2×1的解向量。

运行代码后,X的值为[-4.0000;4.5000]。

第二种方法是用MATLAB函数“inv”求解逆矩阵来解线性方程组。

当系数矩阵A非奇异(可逆)时,可以使用逆矩阵求解线性方程组。

使用“inv”函数的语法如下:X = inv(A) * B其中A是系数矩阵,B是常数向量,X是解向量。

该方法先计算A的逆矩阵,然后将逆矩阵与B相乘得到解向量X。

下面是一个使用“inv”函数解线性方程组的例子:A=[12;34];B=[5;6];X = inv(A) * B;上述代码中,A是一个2×2的系数矩阵,B是一个2×1的常数向量,X是一个2×1的解向量。

运行代码后,X的值为[-4.0000;4.5000]。

第三种方法是用MATLAB函数“mldivide”(或“\”)求解线性方程组。

该函数使用最小二乘法求解非方阵的线性方程组。

使用“mldivide”函数的语法如下:X=A\B其中A是系数矩阵,B是常数向量,X是解向量。

求线性方程组AX=b通解的Matlab实现程序

求线性方程组AX=b通解的Matlab实现程序
而得 到 齐次线 性 方 程组 A X= 0 的所有 解 所构 成 的空 间 ,也就 是齐次 线性方 程组 的一 个基础解 系 。 对齐 次 线性 方 程 组Ax = 0 ,Ma t l a b 命 令 如下 :( 1 )  ̄ l f 果 A
[ n 1 】 [ 2 n 2 】
的通解加 上其 自身 的一个 特解 。在理 论基础 上 ,我们 利用下 面 的例子说 明Ma t 1 a b 实现 程序 。
解 :Ma t l a b 实现程序如下:
> > A = [ 1 1 0— 3 - 1 ; 1— 1 2— 1

7 ] ;
> > f o r ma t r a t

下面是 其简化形 式
[一 n l +7 n 2 ]
[n 1 +5 n 2 ]
C1 1 +C2 ( z 2+… +C


0 【

, ,
在这里 C , C 2 , …, C ~ 设为任 意 的常数 。 利用Ma t l a b ,我 们要 求零 空 间 ,可 以调用 函数n u l l ,从
%限定 输 出格式 为有理 式
I 2 十 X 2 一 x 3 + x 4 = 1 I 3 一 2 + 2 一 3 = 2
> >P=n u l l ( A , ' r ’ ) % 求线性方程组的解空间的有理形式
的基
> > s y ms n l ,n 2
例 2 求 方 程 组 : 1 5 + 一 X 3 + 2 X 4 : 一 1 的 通 解。
关键 词 :齐次线性 方程 组 ;非齐次线性 方程组 ;通解
引 言
求 解 线性 方 程 组 的 问题 是 数 值线 性 代 数 的三 大 问 题之

数值分析(hilbert矩阵)病态线性方程组的求解matlab程序

数值分析(hilbert矩阵)病态线性方程组的求解matlab程序

(Hilbert 矩阵)病态线性方程组的求解理论分析表明,数值求解病态线性方程组很困难。

考虑求解如下的线性方程组的求解Hx = b ,期中H 是Hilbert 矩阵,()ij n n Hh ,11ij h i j ,i ,j = 1,2,…,n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ,取(1,1,,1)nx K ?,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,比较结果。

3.结合计算结果,试讨论病态线性方程组的求解。

第1小题:condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出,1)在n小于等于13之前,图像近似直线log(cond(H))~1.519n-1.8332)在n大于13之后,图像趋于平缓,并在一定范围内上下波动,同时随着n的增加稍有上升的趋势第2小题:solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多,迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多,迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多,迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多,迭代可能不收敛');break;endend。

Matlab解方程(方程组)

Matlab解方程(方程组)

Matlab 解方程这里系统的介绍一下关于使用Matlab求解方程的一系列问题,网络上关于Matlab求解方程的文章数不胜数,但是我大体浏览了一下,感觉很多文章都只是零散的介绍了一点,都只给出了一部分Matlab函数例子,以至于刚接触的人面对不同文章中的不同函数一脸茫然,都搞不清楚这些函数各自的用途,也不知道在什么样的情况下该选择哪个函数来求解方程,在使用Matlab解方程时会很纠结。

不知道读者是否有这样的感觉,反正我刚开始接触时就是这样的感觉,面对网络搜索到一系列函数都好想知道他们之间是个什么关系。

所谓的方程就是含有未知数的等式,解方程就是找出使得等式成立时的未知数的数值。

求方程的解可以转换成不同形式,比如求函数的零点、多项式的根。

方程分类很多,按照未知数个数分为一元、二元、多元方程;按照未知数组合形式分为线性方程和非线性方程;按照非零项次数是否一致分为齐次方程和非齐次方程。

线性方程就是方程中未知数次数是一次的,未知数之间不存在指、对、2及以上幂次的关系,线性方程又分为一元线性方程,也就是一元一次方程;多元线性方程,也就是多元一次方程,多以线性方程组的形式出现(包括齐次线性方程组和非齐次线性方程组)。

在Matlab中求解方程的函数主要有roots、solve、fzero、和fsolve函数等,接下来详细的介绍一下各个Matlab函数的使用方法和使用场合。

一、直接求解法(线性方程组)直接求解法不需要借助任何的Matlab函数,主要用于求解线性方程组,也就是未知数次数是一次的方程组,包括齐次线性方程组合非齐次线性方程组。

当然既然可以求解方程组自然也就可以求解单个方程。

主要针对A x=b形式的方程,其中A是未知数系数矩阵,x是未知数列向量,b是常数列向量,当b=0时就是齐次线性方程组,b ≠0时是非齐次线性方程组。

用左除法,x=A\b例:求解线性方程组的解12341242341234251357926640x x x x x x x x x x x x x x +-+=⎧⎪-+=-⎪⎨+-=⎪⎪+--=⎩解:即直接利用b 左除以A 。

Matlab数组运算及线型方程组的求解

Matlab数组运算及线型方程组的求解

数组运算及线型方程组的求解1.“:”号的用法。

用“:”号生成行向量a=[1 2 3 4 5 6 7 89 10]、b=[5 3 1 -1 -3 -5];用线性等分命令linspace重新生成上述的a和b向量。

另,在100和10000之间用对数等分命令logspace生成10维的向量c。

linspace(1,10,10) linspace(5,-5,6) ak=logspace(2,4,10)2. 已知多项式a(x)=x2+2x+3,b(x)=4x2+5x+6,求a,b积的微分。

a=[1 2 3];b=[4 5 6];c=polyder(a,b)c=poly2str(c,'x')3.用生成下列矩阵,取出方框内的数组元素a=[1:5;10:-1:6;11:15;16:20;21:25]q=a(2,2:3)m=a(2:4,4)n=a(4:5,1:3)4. 生成一个9×9维的魔方矩阵,提取其中心的3×3维子矩阵M,利用sum函数检验其各行和各列的和是否相等。

并且实现上述中心矩阵左旋90°或右旋90°,左右翻转,上下翻转a=magic(9)a =47 58 69 80 1 12 23 34 45 57 68 79 9 11 22 33 44 46 67 78 8 10 21 32 43 54 56 77 7 18 20 31 42 53 55 66 6 17 19 30 41 52 63 65 76 16 27 29 40 51 62 64 75 5 26 28 39 50 61 72 74 4 1536 38 49 60 71 73 3 14 2537 48 59 70 81 2 13 24 35>> m=a(4:6,4:6)m =20 31 42 30 41 52 40 51 62 >> c=rot90(m)c =42 52 62 31 41 51 20 30 40 >> c=rot90(m,-1)c =40 30 20 51 41 31 62 52 42 >> s=fliplr(m)s =42 31 20 52 41 30 62 51 40 >> w=flipud(m)w =40 51 62 30 41 52 20 31 425.已知a=[1 2 3:4 5 6;7 8 0],求其特征多项式并求其根、特征值和特征多项式6. 计算二重不定积分7.求解微分方程。

matlab矩阵分解与线性方程组求解

matlab矩阵分解与线性方程组求解

格式
[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-平方根法和改进平方根法求解线性方程组例题与程序演示教学

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解线性方程组

数值计算实验——解线性方程组西南交通大学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程序(牛顿法及线形方程组)

MATLAB 程序1、图示牛顿迭代法(M 文件)文件名:newt_gfunction x = new_g(f_name,x0,xmin,xmax,n_points)clf,hold off% newton_method with graphic illustrationdel_x = 0.001;wid_x = xmax - xmin; dx = (xmax - xmin)/n_points;xp = xmin:dx:xmax; yp = feval(f_name,xp);plot(xp,yp);xlabel('x');ylabel('f(x)');title('newton iteration'),hold onymin = min(yp); ymax = max(yp); wid_y = ymax-ymin;yp = 0. * xp; plot(xp,yp)x = x0; xb = x+999; n=0;while abs(x-xb) > 0.000001if n > 300 break; endy=feval(f_name,x); plot([x,x],[y,0]);plot(x,0,'o')fprintf(' n = % 3.0f, x = % 12.5e, y = % 12.5e \ n', n, x, y);xsc = (x-xmin)/wid_x;if n < 4, text(x,wid_y/20,[num2str(n)]), endy_driv = (feval(f_name,x + del_x) - y)/del_x;xb = x;x = xb - y/y_driv; n = n+1;plot([xb,x],[y,0])endplot([x x],[0.05 * wid_y 0.2 * wid_y])text( x, 0.2 * wid_y, 'final solution')plot([ x ( x - wid_x * 0.004)], [0.01 * wid_y 0.09 * wid_y])plot([ x ( x + wid_x * 0.004)], [0.01 * wid_y 0.09 * wid_y])传热问题假设一个火炉是用厚度为0.05m 的砖单层砌成的。

matlab超松弛迭代法求方程组

matlab超松弛迭代法求方程组

一、介绍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解方程组方法

matlab解方程组方法在MATLAB中,有多种方法可以解方程组。

以下是其中几种常用的方法:1.solve函数:这是最直接的方法,适用于解线性方程组。

假设你有以下线性方程组:(Ax = b)你可以使用solve函数来求解。

例如:2.matlab复制代码A = [1, 2; 3,4];b = [5; 6];x = solve(A,b);3.\和/运算符:这两个运算符也可以用于解线性方程组。

例如:4.matlab复制代码A = [1, 2; 3, 4];b = [5; 6];x = A\b; % 使用左除运算符或者matlab复制代码x = b/A; % 使用右除运算符5.gaussj函数:这个函数使用高斯-约当消元法来解方程组。

使用方法如下:6.matlab复制代码A = [1, 2; 3,4];b = [5; 6];x = gaussj(A,b);7.mldivide函数:这个函数与\运算符相同,也是用于解线性方程组。

例如:8.matlab复制代码A = [1, 2; 3, 4];b = [5; 6];x = mldivide(A, b); % 等价于A\b9.lyap函数:对于非线性方程组,可以使用lyap函数来求解。

这个函数用于解决Lyapunov方程,通常用于控制系统和稳定性分析。

使用方法如下:10.matlab复制代码A = [1, 2; 3, 4];lyap(A); % 对于给定的A矩阵,求解Lyapunov方程。

11.fzero和root函数:这两个函数用于求解非线性方程的根。

例如,如果你有一个非线性方程(f(x) = 0),你可以使用fzero或root来找到这个方程的根。

使用方法如下:12.matlab复制代码f = @(x) x^2 - 4; % 非线性方程 f(x) = x^2 - 4x = fzero(f, [1, 2]); % 在区间[1,2]内寻找方程的根或者:matlab复制代码root(f) % 使用root函数求解非线性方程的根。

线性方程组及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分解法然后调用求解。

matlab求线性方程组的解

matlab求线性方程组的解

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计算方法3解线性方程组计算解法名师公开课获奖课件百校联赛一等奖课件

MATLAB计算方法3解线性方程组计算解法名师公开课获奖课件百校联赛一等奖课件

li1 ai1
u11
(i 2,3,, n)
k 1
ukj akj lkmumj akj
m 1
(
j
k,
k
1,,
n)
lik
aik
k 1
limumk
m 1
(i
k
1,,
n)
ukk aik
(k 2,3,, n)
例3.1
2 1 2 6 2 1 2 6
4 5 4 18 2 3 0 6
a11 a12 a1n l11
a21
a22
a2n
l21
l22
l11 l21 l n1
l22
l
n2
an1
an2
ann
l n1
l n2
l
nn
l
nn
其中aij a ji
由矩阵乘法
(1)
1)
l2 11
a11
l11
a11
(取正)
2) L第1行 LT第j列 (j 2,,n)
…….
(k)
1求u的第k行:用L的第k行 u的第j列
(j k,k 1,,n)
(lk1 , lk 2 ,, lkk,0,0) (u1 j , u2 j ,, u jj,0,0)' akj
k 1
k 1
lkmumj 1 ukj akj ukj akj lkmumj
m 1
m 1
2 求L的第k列:用L的第i行 u的第k列
利用Gauss消元法得到同解旳三角方程为
1 c1
y1
2 c2
y2
n1
ቤተ መጻሕፍቲ ባይዱ
cn1

matlab练习程序(线性常微分方程组参数拟合)

matlab练习程序(线性常微分方程组参数拟合)

matlab练习程序(线性常微分⽅程组参数拟合)⽐如我们已经有了微分⽅程模型和相关数据,如何求模型的参数。

这⾥以SEIR模型为例⼦,。

⼀般的线性⽅程我们可以⽤,⼀般的⾮线性⽅程我们可以⽤。

这⾥是线性微分⽅程组,所以我们采⽤最⼩⼆乘来解。

关键是构造出最⼩⼆乘形式,微分可以通过前后数据差分的⽅法来求。

不过这⾥还有⼀个技巧就是如果数据前后帧间隔过⼤,可以先插值,再对插值后的数据差分。

如果实际测量数据抖动过⼤导致插值后差分明显不能反映实际情况,可以先对数据平滑(拟合或是平均)再求差分。

先看SEIR微分⽅程:写成矩阵形式:到这⾥就能⽤最⼩⼆乘来求解了。

matlab代码如下:main.m:clear all;close all;clc;%%SEIR模型A = [0.50.10.050.02];[t,h] = ode45(@(t,x)SEIR(t,x,A),[0300],[0.010.980.010]); %[初始感染⼈⼝占⽐初始健康⼈⼝占⽐初始潜伏⼈⼝占⽐初始治愈⼈⼝占⽐]plot(t,h(:,1),'r');hold on;plot(t,h(:,2),'b');plot(t,h(:,3),'m');plot(t,h(:,4),'g');legend('感染⼈⼝占⽐I','健康⼈⼝占⽐S','潜伏⼈⼝占⽐E','治愈⼈⼝占⽐R');title('SEIR模型')data=[t h];data = data(1:3:80,:); %间隔取⼀部分数据⽤来拟合figure;plot(data(:,1),data(:,2),'ro');hold on;plot(data(:,1),data(:,3),'bo');plot(data(:,1),data(:,4),'mo');plot(data(:,1),data(:,5),'go');T=min(data(:,1)):0.1:max(data(:,1)); %插值处理,如果数据多,也可以不插值I=spline(data(:,1),data(:,2),T)';S=spline(data(:,1),data(:,3),T)';E=spline(data(:,1),data(:,4),T)';R=spline(data(:,1),data(:,5),T)';plot(T,I,'r.');plot(T,S,'b.');plot(T,E,'m.');plot(T,R,'g.');%求微分,如果数据帧间导数变化太⼤,可以先平均或者拟合估计⼀个导数%因为前⾯T是以0.1为步长,这⾥乘以10dI = diff(I)*10; dI=[dI;dI(end)];dS = diff(S)*10; dS=[dS;dS(end)];dE = diff(E)*10; dE=[dE;dE(end)];dR = diff(R)*10; dR=[dR;dR(end)];X = [zeros(length(I),1) -I.*S zeros(length(I),2); %构造线性最⼩⼆乘⽅程组形式-E I.*S -E zeros(length(I),1);E zeros(length(I),2) -I;zeros(length(I),2) E I];Y = [dS;dE;dI;dR];A = inv(X'*X)*X'*Y%⽤估计参数代⼊模型[t,h] = ode45(@(t,x)SEIR(t,x,A),[0300],[I(1) S(1) E(1) R(1)]); %[初始感染⼈⼝占⽐初始健康⼈⼝占⽐初始潜伏⼈⼝占⽐初始治愈⼈⼝占⽐] plot(t,h(:,1),'r');hold on;plot(t,h(:,2),'b');plot(t,h(:,3),'m');plot(t,h(:,4),'g');SEIR.m:function dy=SEIR(t,x,A)alpha = A(1); %潜伏期转阳率beta = A(2); %感染率gamma1 = A(3); %潜伏期治愈率gamma2 = A(4); %患者治愈率dy=[alpha*x(3) - gamma2*x(1);-beta*x(1)*x(2);beta*x(1)*x(2) - (alpha+gamma1)*x(3);gamma1*x(3)+gamma2*x(1)];结果:原始参数[0.5 0.1 0.05 0.02]与模型:拟合参数[0.499921929359668 0.100099242849855 0.0505821757746970 0.0199739921888752]与模型:。

线性代数方程组数值解法及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与线性方程组求解

科学从来不是以繁琐为特征的,工程问题则是很繁琐的,所以工程师需要计算机辅助,让计算机去处理大量的简单重复性工作,让人可以集中精力处理重要的关键性决策。

以现代飞行器外形设计为例,它决定了整个飞行器的空气动力学特征,因而地位十分重要。

现在流体动力学的理论和计算流体动力学软件已经很成熟,问题是如何应用到特定的外形上来。

人们采用的方法就是把飞行器的外形分成若干大的部件,每个部件沿着其表面用三维的细网格划分出许多立方体,这些立方体包括了机身表面以及此表面内外的空气。

对每个立方体列出空气动力学方程,其中包括了与它相邻的立方体的共同边界变量,这些方程通常都已经简化为线性方程。

对一个飞行器,小立方体的数目可以多达400,000个,而要解的联立方程可能多达2,000,000个。

对如此大的方程组在求解之前先要简化,简化方法包括:一是利用许多不相邻的元素之间没有关联,其交叉系数为零,在整个大联立方程组中,绝大部分的系数为零,使用稀疏矩阵的计算方法;二是把矩阵进行分解如LU分解,可以大大提高计算速度。

工程问题不断向矩阵理论提出需求,从而推动了理论的发展。

绝不是某些超人,从定义出发,冥思苦想出新理论的。

卫星遥感图象处理中,气象和地球资源卫星大约用90分钟绕地球一圈,拍摄的图像宽度约为150公里。

地球赤道长40,000公里,因此大约每16天,可以扫描地球的每个角落一遍。

卫星上用三种可见光和四种红外光进行摄像,对每一个区域,可以获得七张遥感图象。

利用多通道的遥感图可以获取尽可能多的地面信息,因为各种地貌、作物和气象特征可能对不同波段的光敏感。

而在实用上应该寻找每一个地方的主因素,再把它们合成起来,成为一张实用的图象。

每一个象素上有七个数据,形成一个多元的变量数组,在其中合成并求取主因素的问题,就与线性代数中要讨论的特征值问题有关。

国家地理信息系统在全国设立几十万个观察点,把每一点的经度、纬度和高度三个坐标建立起来。

由于地壳的变动,测量仪器的现代化和实际需求的增长,地理信息系统精细化的工程是与无止境的。

MATLAB解方程组(线性与非线性方程组)

MATLAB解方程组(线性与非线性方程组)

例7-9 求下列非线性方程组在(0.5,0.5) 附近的数值解。 (1) 建立函数文件myfun.m。 function q=myfun(p) x=p(1); y=p(2); q(1)=x-0.6*sin(x)-0.3*cos(y); q(2)=y-0.6*cos(x)+0.3*sin(y); (2) 在给定的初值x0=0.5,y0=0.5下,调用fsolve函数求方程的根。 x=fsolve('myfun',[0.5,0.5]',optimset('Display','off')) x= 0.6354 0.3734
2.Gauss-Serdel迭代法 在Jacobi迭代过程中,计算时,已经得到,不必再用,即原来的迭代
公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b, 于是得到:
x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比,Gauss-Serdel
7.1.2 迭代解法 迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中,迭代
解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法 和两步迭代法。
1.Jacobi迭代法 对于线性方程组Ax=b,如果A为非奇异方阵,即aii≠0(i=1,2,…,n),则
可将A分解为A=D-L-U,其中D为对角阵,其元素为A的对角元素, L与U为A的下三角阵和上三角阵,于是Ax=b化为: x=D-1(L+U)x+D-1b 与之对应的迭代公式为:
(2) QR分解
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 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-Jordan消去法: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);endfor 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%%%%Endfunction x=lusolv(a,b)%LU分解求解线性方程组aX=b if 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)^2;endif delta<1e-10return;endL(k,k)=sqrt(delta);for 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分解求解线性方程组Ax=b 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);endfor i=n:-1:1for 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)%利用对对称正定矩阵A进行LDL'分解法求解线性方程组Ax=b n=length(b);[l,d]=LDL_Factor(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)/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)x =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.0010n =5Jacobi迭代法: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,-1); %求A的下三角阵U=-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.97391.0010n =5Gauss-Seidel迭代法:function [x,n]=gauseidel(A,b,x0,eps,M) if nargin==3eps= 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的上三角阵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.9739-0.00471.0010n =4超松驰迭代法:function [x,n]=SOR(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的上三角阵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:>> [x,n]=SOR(A,b,x0,1)x =0.9739-0.00471.0010n =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.0010n =3最速下降法:function [x,n]=fastdown(A,b,x0,eps) if(nargin == 3)eps = 1.0e-6;endr = 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.0010n =5共轭梯度法:function [x,n]=conjgrad(A,b,x0) if(nargin == 3)eps = 1.0e-6;endr1 = 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;endd = dot(r2,r2)/dot(p2,A*p2);x = x+d*p2;n = n + 1;Example:>> [x,n]=conjgrad(A,b,x0)x =0.9739-0.00471.0010n =4预处理的共轭梯度法:当AX=B为病态方程组时,共轭梯度法收敛很慢。

相关文档
最新文档