追赶法求解三对角线性方程组
- 格式:doc
- 大小:150.00 KB
- 文档页数:9
Matlab追赶法和迭代法解线性⽅程组实验⽬的:1)追赶法解三对⾓阵;2)掌握解线性⽅程组的迭代法;3)⽤Matlab实现Jacobi及超松弛迭代法实验要求:1)掌握追赶法解三对⾓阵2)掌握解线性⽅程组的迭代法3)提交追赶法、Jacobi及超松弛迭代法的m⽂件实验内容:1)追赶法解三对⾓矩阵⽅程(m⽂件)习题1. ⽤追赶法的m⽂件求解2)Jacobi迭代法解线性⽅程组(m⽂件)对不同初值⽤Jacobi迭代法解习题1并⽐较结果。
3)超松弛迭代法解线性⽅程组(m⽂件)对不同松弛因⼦解习题1并⽐较结果。
实验步骤: 代码:1 %追赶法2 %输⼊:系数矩阵A和因变量d;3 %输出:⾃变量x4 function z=zuigan(A,d)5 n=length(d);6 %取三对⾓元素a,b,c7for i=1:n-18 a(i)=A(i,i);9 b(i)=A(i+1,i);10 c(i)=A(i,i+1);11 end12 a(n)=A(n,n);13 %分解系数矩阵14 u(1)=a(1);15 l(1)=c(1)/a(1);16for i=2:n-117 u(i)=a(i)-b(i-1)*l(i-1);18 l(i)=c(i)/u(i);19 end20 u(n)=a(n)-c(n-1)*l(n-1);21 %解y22 y(1)=d(1)/u(1);23for k=2:n24 y(k)=d(k)-c(k-1)*y(k-1)/u(k);25 end26 %解x27 x(n)=y(n);28for k=n-1:-1:129 x(k)=y(k)-l(k)*x(k+1);30 end31 z=x;32 endzuigan 运⾏: 所得结果,较为粗糙。
代码:1 %雅克⽐迭代法2 %输⼊系数矩阵A,因变量b,初始向量x0,容许误差eps,最⼤迭代次数t3 %输出⾃变量x和迭代数n4 function [z,k]=jacobi(A,b,x0,e,t)5 %默认eps和最⼤迭代次数m6if nargin==37 e=1e-6;8 m=200;9 elseif nargin<310 error('输⼊的参数不⾜');11return;12 elseif nargin==513 m=t;14 end15 n=length(b);16 x(1,:)=x0;17 z(1,:)=x0;18for k=2:m19 sum=0;20for i=1:n21 w=0;22 u=0;23for j=i+1:n24 w=w+A(i,j)*x(k-1,j);25 end26for j=1:i-127 u=u+A(i,j)*x(k-1,j);28 end29 x(k,i)=(-1/A(i,i))*(u+w-b(i));30if sum<abs(x(k,i)-x(k-1,i))31 sum=abs(x(k,i)-x(k-1,i));32 end33 end34if sum<e35 z(k,:)=x(k,:);36return;37 end38 z(k,:)=x(k,:);39 end40 endjacobi 运⾏⽰例,初始向量x0=[0 0 0 0 0 0];和初始向量x0=[1 1 1 1 1 1]; 初始值不同,迭代次数可能不同。
题目:三对角矩阵追赶法在Matlab中的应用一、概述三对角矩阵是一类特殊的矩阵,在科学计算领域中有着广泛的应用。
而追赶法(也称为托马斯算法)是解三对角线性方程组的一种高效方法。
在Matlab中,可以通过内置函数或者编写自定义代码来实现三对角矩阵追赶法,本文将介绍该方法的原理及在Matlab中的应用。
二、三对角矩阵及其特点三对角矩阵是指除了主对角线外,只有上下相邻两条对角线上有非零元素的矩阵。
具体来说,三对角矩阵的形式如下:```a1 b1 0 0 0c1 a2 b2 0 00 c2 a3 b3 00 0 c3 a4 b40 0 0 c4 a5```其中a1, a2, ..., an为主对角线上的元素,b1, b2, ..., bn-1为上对角线上的元素,c2, c3, ...,为下对角线上的元素。
三对角矩阵具有以下特点:1. 矩阵的维数较大时,存储和计算均占据较少的空间和时间;2. 存在高效的解法,例如追赶法;3. 在有限差分法、对流扩散方程等实际问题中有广泛应用。
三、追赶法原理及算法流程追赶法是一种用于解三对角线性方程组的常用方法,其原理基于LU 分解和回代求解。
算法流程如下:1. 首先对三对角矩阵进行LU分解,将其分解为一个单位下三角矩阵L和一个上三角矩阵U;2. 然后进行前代和回代,分别求解Ly=d和Ux=y。
这两个步骤可以通过追赶法来高效求解。
四、Matlab中的三对角矩阵追赶法实现在Matlab中,可以通过编写自定义代码来实现三对角矩阵追赶法,也可以利用内置函数进行求解。
下面分别介绍这两种方法:1. 自定义代码实现可以编写自定义函数来实现三对角矩阵追赶法。
下面以一个简单的示例代码来说明实现过程:```matlabfunction x = tridiagonal_solver(a, b, c, d)n = length(d);for i = 2:nm = a(i) / b(i-1);b(i) = b(i) - m * c(i-1);d(i) = d(i) - m * d(i-1);endx(n) = d(n) / b(n);for i = n-1:-1:1x(i) = (d(i) - c(i) * x(i+1)) / b(i);endend```通过上述代码,可以实现对三对角矩阵的追赶法求解。
在探讨MATLAB追赶法解101阶三对角方程组之前,我们首先需要了解什么是追赶法和什么是三对角方程组。
追赶法又称托马斯算法,是一种用于求解带状矩阵(即只有主对角线和两条相邻的对角线上有非零元素的矩阵)的线性方程组的方法。
而三对角矩阵就是只有主对角线和两条相邻的对角线上有非零元素的矩阵。
在实际应用中,求解带状矩阵的线性方程组是非常常见的,特别是在数值计算和科学工程领域。
现在,让我们深入探讨MATLAB追赶法解101阶三对角方程组的方法和具体步骤。
一、MATLAB追赶法解101阶三对角方程组1. 概念介绍101阶三对角方程组是一个非常大的线性方程组,通常使用传统的高斯消元法来求解会耗费大量的时间和计算资源。
而MATLAB追赶法通过利用三对角矩阵的特殊性质,可以有效地简化计算过程,并且节省大量的内存和计算资源。
2. 追赶法步骤(1)将原方程组化为追赶法所需的形式;(2)利用追赶法求解三对角线性方程组。
二、追赶法求解101阶三对角方程组的实现过程1. 将原方程组化为追赶法所需的形式对于101阶三对角方程组,我们首先需要将其化为追赶法所需的形式。
这个过程涉及到选取合适的追赶元和追赶子以及对原方程组的变形,将其化为追赶法能够直接处理的形式。
2. 利用追赶法求解线性方程组一旦将原方程组化为追赶法所需的形式,我们就可以利用追赶法对其进行求解。
追赶法的核心是通过追赶子的迭代计算,逐步求得线性方程组的解。
在MATLAB中,可以使用内置的追赶法求解函数,也可以编写自定义的追赶法算法来实现对101阶三对角方程组的求解。
三、个人观点和理解在实际工程和科学计算中,追赶法是一种非常有效的求解带状矩阵线性方程组的方法。
对于大规模的三对角方程组,特别是高阶的情况,传统的直接求解方法往往会遇到内存和计算资源的限制,而追赶法能够通过精巧的迭代计算,在保证解的精度的显著提高计算效率。
在MATLAB中,通过调用内置的追赶法函数,可以快速地求解大规模的三对角方程组,极大地方便了工程实践中的数值计算工作。
追赶法,高斯消元法,逆矩阵法,迭代法 —— 解线性方程组精仪学院 马金玉 1012202030本文主要详细介绍了追赶法,高斯法,逆矩阵法的方法原理,运用这三种方法分别进行线性方程的求解举例,给出MATLAB 相应程序,最后做结果分析,比较说明追赶法和高斯法的特点。
最后对三种典型迭代方法Jacobi 迭代,Gauss-Seidel 迭代,SOR 迭代进行简单的分析比较。
1. 追赶法1.1).追赶法方法介绍追赶法用于求解以下形式的方程组(三对角方程组)d Ax =其中 1[,,]T n d d =d ,系数矩阵(三对角矩阵)11222111n n n n n b c a bc a b c a b ---⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦A系数矩阵A 的元素满足1100 0 (2,,1)0i i i i i n n b c b a c a c i n b a ⎧>>⎪≥+>≠=-⎨⎪>>⎩第一步:实现A=LU 的分解,按照递推公式1111//()i i i i i c b c b a βββ-=⎧⎨=-⎩ 计算 123,,...........βββ:第二步:求解方程组LY=f,相应的递推公式 11111/()/()i i i i i i i y f b y f a y b a β--=⎧⎨=--⎩ 第三部:求解方程组UX=Y ,相应的递推公式1()n nii i i x y x y x β-=⎧⎨=-⎩ 求得x因为计算1231......n ββββ-→→→→ 及 1231......n y y y y -→→→→的过程是追赶的过程,结出结果X 。
1.2).追赶法解线性方程组的matlab实例解线性方程组第一步:编写M文件如下:function [x,y,beta]=zhuiganfa(a,b,c,f)%a,b,c是三对角阵的对角线上的元素,f是自由项.n=length(b);beta(1)=c(1)/b(1);for i=2:nbeta(i)=c(i)/(b(i)-a(i)*beta(i-1));endy(1)=f(1)/b(1);for i=2:ny(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endx(n)=y(n);for i=n-1:-1:1x(i)=y(i)-beta(i)*x(i+1);enddisp(sprintf('k x(k) y(k) beta(k)')); for i=0:n-1disp(sprintf('%d %15.4f %15.4f %15.4f',i,x(i+1),y(i+1),beta(i+1))); end追赶法M文件程序截图如图1所示图1 追赶法M文件程序截图第二步:根据所求方程,在命令窗口中输入如下命令,并按ENTER 键确认。
3)三对角形线性方程组123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014x x x x x x x x x x -⎡⎤⎡⎤⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎢⎥--⎢⎢⎥⎢⎢⎥-⎣⎦⎣⎦7513261214455⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥-⎣⎦*(2,1,3,0,1,2,3,0,1,1)T x =--- 二、数学原理设系数矩阵为三对角矩阵1122233111000000000000000n n n nn b c a b c a b A a b c a b ---⎛⎫ ⎪ ⎪ ⎪=⎪ ⎪ ⎪⎪ ⎪⎝⎭则方程组Ax=f 称为三对角方程组。
设矩阵A 非奇异,A 有Crout 分解A=LU ,其中L 为下三角矩阵,U 为单位上三角矩阵,记1122233110000100000001000000100,00000000000001n n nn b L U γαβγββγβ--⎛⎫⎛⎫ ⎪⎪ ⎪ ⎪ ⎪ ⎪∂==⎪⎪ ⎪⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪∂⎝⎭⎝⎭可先依次求出L ,U 中的元素后,令Ux=y ,先求解下三角方程组Ly=f 得出y ,再求解上三角方程组Ux=y 。
事实上,求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑。
其计算公式为:1111,1111,111,2,3,,,1,2,,1ii i i i i i i ii i i i i n ni i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+⎧===⎪⎪=⎪⎪⎪==-=⎪⎪⎨-⎪=⎪⎪=⎪⎪=--⎪=-⎪⎩对对(*)三、程序设计function x=chase(a,b,c,f)%求解线性方程组Ax=f,其中A 是三对角阵 %a 是矩阵A 的下对角线元素a(1)=0 %b 是矩阵A 的对角线元素%c 是矩阵A 的上对角线元素c(n)=0 %f 是方程组的右端向量 n=length(f);x=zeros(1,n);y=zeros(1,n); d=zeros(1,n);u= zeros(1,n); %预处理 d(1)=b(1); for i=1:n-1 u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i); end%追的过程y(1)=f(1)/d(1); for i=2:ny(i)=(f(i)-a(i)*y(i-1))/d(i); end%赶的过程 x(n)=y(n); for i=n-1:-1:1x(i)=y(i)-u(i)*x(i+1); end>> a=[0,-1,-1,-1,-1,-1,-1,-1,-1,-1];>> b=[4,4,4,4,4,4,4,4,4,4];>> c=[-1,-1,-1,-1,-1,-1,-1,-1,-1,0];>> f=[7,5,-13,2,6,-12,14,-4,5,-5];>> x=chase(a,b,c,f)x =2.00001.0000-3.00000.00001.0000-2.00003.0000-0.00001.0000-1.0000四、结果分析和讨论追赶法求解的结果为x=(2,1,-3,0,1,-2,3,0,1,-1)T。
追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容1、 学习和理解追赶法求解三对角线性方程组的原理及方法;2、 利用MA TLAB 编程实现追赶法;3、 举例进行求解,并对结果进行分。
三 实验原理设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦………… 这种方程组称为三对角线性方程组。
显然,A 是上下半宽带都是1的带状矩阵。
设A 的前n-1个顺序主子式都不为零,根据定理2.5的推论,A 有唯一的Crout 分解,并且是保留带宽的。
其中L 是下三角矩阵,U 是单位上三角矩阵。
利用矩阵相乘法,可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………得到:由上列各式可以得到L 和U 。
引入中间量y ,令yUx =,则有:已知L 和d ,可求得y 。
则可得到y 的求解表达式:11/12,3,,()(1)*y()=()[()(1)]/y d l i nm i y i li i di y i di m i y i li==-+=--…1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui=*===≤≤+=+++≤≤-=∙=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦……………由y Ux =得:111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦………… 可得到X 的求解表达式:()()1,2,,1()()u()(1)x n y n i n n x i y i i x i ==--=-+… 从而得到Ax=d 的解x 。
托马斯方法方程组
托马斯方法(Thomas decomposition)或追赶法是一种求解三对角线性
方程组的快速算法。
这种算法在特殊情况下,即对角占优的情况下,其复杂度可以与N同阶,为O(3N-2)。
三对角线性方程组是一类特殊的线性方程组,其系数矩阵是一个三对角矩阵。
这种方程组在物理学、工程学和经济学等领域有广泛应用。
托马斯方法的基本思想是利用三对角矩阵的特性,通过一系列迭代步骤求解方程组。
在使用托马斯方法时,需要先输入方程的个数,然后根据算法步骤逐步求解。
具体实现时,需要定义一个类,包含输入函数和求解函数等成员函数,以及相应的私有变量。
在输入函数中,需要读取方程的系数和常数项,并存储在相应的数组中。
在求解函数中,需要实现托马斯方法的算法步骤,包括初始化、迭代计算和输出结果等。
托马斯方法是一种高效的求解三对角线性方程组的算法,其优点在于具有较低的时间复杂度和空间复杂度。
然而,它也有一定的局限性,只能适用于对角占优的三对角线性方程组。
对于其他类型的线性方程组,可能需要使用其他算法进行求解。
3)三对角形线性方程组123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014x x x x x x x x x x -⎡⎤⎡⎤⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎢⎥--⎢⎢⎥⎢⎢⎥-⎣⎦⎣⎦7513261214455⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥-⎣⎦*(2,1,3,0,1,2,3,0,1,1)T x =--- 二、数学原理设系数矩阵为三对角矩阵112223311100000000000000000n n n n n b c a b c a b A a b c a b ---⎛⎫ ⎪⎪ ⎪=⎪ ⎪ ⎪⎪⎪⎝⎭L L L M M MM M M L L则方程组Ax=f 称为三对角方程组。
设矩阵A 非奇异,A 有Crout 分解A=LU ,其中L 为下三角矩阵,U 为单位上三角矩阵,记11222331100001000000010000000100,00000000000001n n nn b L U γαβγββγβ--⎛⎫⎛⎫⎪ ⎪⎪ ⎪ ⎪ ⎪∂==⎪⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎪ ⎪∂⎝⎭⎝⎭L L L LL L M M M MM M M MM M L L LLL可先依次求出L ,U 中的元素后,令Ux=y ,先求解下三角方程组Ly=f 得出y ,再求解上三角方程组Ux=y 。
事实上,求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑。
其计算公式为:1111,1111,111,2,3,,,1,2,,1ii i i i i i i ii i i i i n ni i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+⎧===⎪⎪=⎪⎪⎪==-=⎪⎪⎨-⎪=⎪⎪=⎪⎪=--⎪=-⎪⎩L L 对对(*)三、程序设计function x=chase(a,b,c,f)%求解线性方程组Ax=f,其中A 是三对角阵 %a 是矩阵A 的下对角线元素a(1)=0 %b 是矩阵A 的对角线元素%c 是矩阵A 的上对角线元素c(n)=0 %f 是方程组的右端向量 n=length(f);x=zeros(1,n);y=zeros(1,n); d=zeros(1,n);u= zeros(1,n); %预处理 d(1)=b(1); for i=1:n-1 u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i); end%追的过程y(1)=f(1)/d(1); for i=2:ny(i)=(f(i)-a(i)*y(i-1))/d(i); end%赶的过程 x(n)=y(n); for i=n-1:-1:1x(i)=y(i)-u(i)*x(i+1); end>> a=[0,-1,-1,-1,-1,-1,-1,-1,-1,-1];>> b=[4,4,4,4,4,4,4,4,4,4];>> c=[-1,-1,-1,-1,-1,-1,-1,-1,-1,0];>> f=[7,5,-13,2,6,-12,14,-4,5,-5];>> x=chase(a,b,c,f)x =2.00001.0000-3.00000.00001.0000-2.00003.0000-0.00001.0000-1.0000四、结果分析和讨论追赶法求解的结果为x=(2,1,-3,0,1,-2,3,0,1,-1)T。
matlab追赶法求解三对角方程组追赶法是一种求解三对角方程组的有效方法,可以通过简化矩阵的求解过程来提高计算效率。
以下是使用追赶法求解三对角方程组的具体步骤:1. 将三对角矩阵表示为下三角矩阵L、上三角矩阵U和对角矩阵D的乘积形式:A=LDU,其中L是下三角矩阵,U是上三角矩阵,D是对角矩阵。
2. 将方程组Ax=b转化为LDUx=b。
3. 首先使用前向代入法(forward substitution)解下三角方程Ly=b,其中y是临时向量。
4. 然后使用对角方程Dz=y解决z=D^-1y,其中z是临时向量。
5. 最后使用后向代入法(backward substitution)解上三角方程Ux=z,得到方程组的解x。
追赶法的时间复杂度为O(n),相比于高斯消元法等其他方法,追赶法在求解三对角方程组时具有更快的计算速度。
在MATLAB中可以使用专门的函数tridiag来实现追赶法,示例如下:```matlabfunction x = tridiag_solver(a, b, c, d)n = length(b);% 前向代入for i = 2:nm = a(i) / b(i-1);b(i) = b(i) - m * c(i-1);d(i) = d(i) - m * d(i-1);end% 后向代入x = zeros(n, 1);x(n) = d(n) / b(n);for i = n-1:-1:1x(i) = (d(i) - c(i) * x(i+1)) / b(i);endend```其中a、b和c分别是三对角方程组的次对角线、主对角线和超对角线上的元素,d是方程组的右侧常数向量。
函数返回方程组的解x。
追赶法求解三对角方程组要求:对于给定的三对角系数矩阵和右端项,可以求解线性代数方程组一、追赶法的数学理论设系数矩阵为三对角矩阵则方程组Ax=f称为三对角方程组。
设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记可先依次求出L,U中的元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y。
事实上,求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑。
其计算公式为:(*)二、追赶法的算法和流程图1.预处理生成方程组的系数及其除数,事实上,按式(*)可交替生成与:→→→…→→其计算公式为2.追的过程顺序生成方程组右端:→→…→据式(*)的计算公式为3.赶的过程逆序得出方程组的解:→→…→其计算公式按式为三、追赶法的Matlab实现function x=chase(a,b,c,f)%求解线性方程组Ax=f,其中A是三对角阵%a是矩阵A的下对角线元素a(1)=0%b是矩阵A的对角线元素%c是矩阵A的上对角线元素c(N)=0%f是方程组的右端向量N=length(f);x=zeros(1,N);y=zeros(1,N);d=zeros(1,N);u= zeros(1,N);%预处理d(1)=b(1);for i=1:N-1u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i);end%追的过程y(1)=f(1)/d(1);for i=2:Ny(i)=(f(i)-a(i)*y(i-1))/d(i);end%赶的过程x(N)=y(N);for i=N-1:-1:1x(i)=y(i)-u(i)*x(i+1);end四、追赶法的算例实现算例用追赶法求解方程组解答令a=[0,-1,-1,-3]; b=[2,3,2,5]; c=[-1,-2,-1,0]; f=[6,1,0,1];在命令窗口运行语句 x=chase(a,b,c,f) 得结果为x=5 4 3 2。
数值计算方法实验报告实验序号:实验二实验名称:追赶法算法设计及MATLAB实现实验人:专业年级:教学班:学号:实验时间:实验二追赶法算法设计及MATLAB实现一、实验目的1.初步掌握算法设计规则;2.初步掌握MATLAB程序设计规则.二、实验内容1.构造利用追赶法求解三对角线性方程组的算法;2.在MATLAB环境下编写追赶法的程序(函数);3.自由选择若干个三对角线性方程组求解。
三、实验步骤1.追赶法算法:算法名称:thomas输入参数:向量a,b,c,f输出参数:输出解信息x算法的自然语言:Step1:u1=b1,y1=b1;Step2:对于i=2,3,….n;Step2.1:当u1-i≠,否则转step5l i =ai/u1-i;ui =bi-li*c1-i;yi =fi-li*y1-i;Step3:当un≠时,xn=yn/un,否则转step5Step4:对于:i=n-1,n-2,…..,2,1,转step6xi =(yi-ci*x1+i)/uiStep5:无解信息,转step7Step6:输出xStep7:关机2.MATLAB程序function [x,L,U]=thomas(a,b,c,f) n=length(b);% 对A进行分解u(1)=b(1);for i=2:nif(u(i-1)~=0)l(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);elsebreak;endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;% 求解Ly=by(1)=f(1);for i=2:ny(i)=f(i)-l(i-1)*y(i-1);end% 求解Ux=yif(u(n)~=0)x(n)=y(n)/u(n);endfor i=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1))/u(i);end3.求解实例例1.方程组例2.方程组例3.方程组⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎭⎫⎝⎛⎪⎪⎪⎪⎪⎭⎫⎝⎛1131132132134321xxxx四、实验结论对于追赶法我最先写的是如下的程序:但是出现了如上截图中的错误,后来与同学讨论还是没能解决我的问题,最后借鉴了她的算法得到了正确的结果。
追赶法求解三对角线性方程组
一 实验目的
利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容
1、 学习和理解追赶法求解三对角线性方程组的原理及方法;
2、 利用MATLAB 编程实现追赶法;
3、 举例进行求解,并对结果进行分。
三 实验原理
设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵
11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦
………… 这种方程组称为三对角线性方程组。
显然,A 是上下半宽带都是1的带状矩阵。
设A 的前n-1个顺序主子式都不为零,根据定理2.5的推论,A 有唯一的Crout 分解,并且是保留带宽的。
其中L 是下三角矩阵,U 是单位上三角矩阵。
利用矩阵相乘法,可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………
得到:
由上列各式可以得到L 和U 。
引入中间量y ,令
y Ux =,则有:
已知
L 和d ,可求得y 。
则可得到y 的求解表达式:
11/1
2,3,,()(1)*y()=()[()(1)]/y d l i n
m i y i li i di y i di m i y i li
==-+=--…
1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui
=*===≤≤+=+++≤≤-=•=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
……………
由y Ux =得:
111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
………… 可得到X 的求解表达式:
()()
1,2,,1
()()u()(1)
x n y n i n n x i y i i x i ==--=-+…
从而得到Ax=d 的解x 。
四 Matlab 编程
根据以上的实验原理,在Matlab 中编程如下函数x=Trid (A,d )
在Matlab中新建Trid.m,其中程序为如上虚线框内代码,放在工
作目录。
在Command Window输入以下语句:
clear all;clc;
fprintf('输入非奇异三对角系数矩阵A\n');
A=input('A matrix='); % 输入系数矩阵
fprintf('系数矩阵');A
if det(A)==0 % 判断系数矩阵是否奇异
fprintf('系数矩阵A奇异!!!请重新输入!\n');
fprintf('重新输入非奇异三对角系数矩阵A\n');
A=input('A matrix=');
fprintf('系数矩阵');A
end
fprintf('输入矩阵d\n');
d=input('d matrix=');
fprintf('矩阵d');d
x=Trid(A,d) % 调用Trid.m中的Trid函数进行求解
fid = fopen('Ax=d.txt', 'wt'); % 生成Ax=d.txt文件
fprintf(fid,'%s\r\n','利用三对角线追赶法求解Ax=d');
fprintf(fid,'%s\r\n','==================================');
for i=1:size(A,1)
fprintf(fid, '%.1f\t', A(i,:)); % 输出Ax=d,以上=为分隔符fprintf(fid, '%s', 'x');
fprintf(fid, '%u\t', i);
fprintf(fid, '%.1f\n', d(i));
end
fprintf(fid,'%s\r\n','==================================');
五 举例计算及分析
以课本(《数值分析》第4版,颜庆津,北京航空航天大学出版社)27页例3为例进行计算,输入系数矩阵A 和d :
41141=4114A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
1411,10.5=132d ⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦ 调用x=Trid(A,d)后,并生成Ax=d.txt 文件,Command Window
同时也会输出解为0.20.2=0.50.80.3x ⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦
,与课本答案一致; Ax=d.txt 文件内容
如下图:
fprintf(fid,'%s\r\n','求解得到结果如下:'); for i=1:size(A,1)
fprintf(fid,'%s','x'); % 输出解向量x (i )
fprintf(fid, '%u', i);
fprintf(fid, '%s\t', '='); fprintf(fid,'%.5f\r\n',x(i)); end fclose(fid)
再次对程序进行验证,输入矩阵如下:
12213=120.531A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦,21=13d ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦
利用matlab 计算Ax 得到Ax=d ,验证所得到的x 即为方程组的解。
由以上两组计算可表明,该程序能满足追赶法解三对角线性方程
组,其中要求系数矩阵满足要求,即非奇异,三对角且前(n-1)个顺序主子式都不为零。
在又换了一组系数矩阵后,系数矩阵如下:
A=[1 1 0 0 0;2 2 1 0 0;0 2 3 1 0;0 0 2 4 1;0 0 0 2 5]
即11221=23124125A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
,11=111d ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ 运行程序后出错,查找问题发现A 进行Crout 分解时,
l1=1
m2=2
u1=1
l2=a2-m2*u1=2-2×1=0
再计算u2=c2/l2时无法计算,分析发现A 的2阶顺序主子式为零,无法进行Crout 分解,所以该程序还无法判断输入的系数矩阵的(n-1)阶顺序主子式是否为零,这是不足的地方,也无法判断系数矩阵是否是三对角矩阵,还需要在程序上进行修改完善。
六 实验总结
1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………
通过本次的课程实验设计,对追赶法求解三对角线性方程组有了深刻的认识和了解,用matlab编程实现了该解法;同时在实验过程中明白了系数矩阵需要满足一定的要求才能使用该种追赶法方法。