追赶法解三对角矩阵
TDMA追赶法
做三次样条曲线时,需要解三对角矩阵(Tridiagonal Matrices)。
常用解法为Thomas Algorithm,又叫The tridiagonal matrix algorithm (TDMA)。
它是一种基于高斯消元法的算法,分为两个阶段:向前消元forward elimination和回代backward substitution。
本文以一个6乘6矩阵为例,介绍一下使用TDMA的求解过程。
1.范例求解步骤1:将矩阵变为上三角矩阵首先要把上面公式中的系数矩阵变为一个上三角矩阵。
第一行:将上式除以b1:可写作:所以矩阵方程可写为:第二行:将变换后的第一行乘以a2,再与第二行相减,即可消去x1,得:所以新的矩阵方程为:同理可推,第三行:第四行:第五行:第六行:最后得到新的上三角矩阵公式为:步骤2:求解x逆序可以求出,如下:2. 一般性公式:注意:使用TDMA求解,系数矩阵需时diagonally dominant,即:3. 实现代码(C语言)void tdma(float x[], const size_t N, constfloat a[], constfloat b[], float c[]){size_t n;c[0] = c[0] / b[0];x[0] = x[0] / b[0];for (n = 1; n < N; n++) {float m = 1.0f / (b[n] - a[n] * c[n - 1]);c[n] = c[n] * m;x[n] = (x[n] - a[n] * x[n - 1]) * m;}for (n = N - 1; n-- >0; )x[n] = x[n] - c[n] * x[n + 1];}。
三对角方程组的追赶法
2013-2014(1)专业课程实践论文题目:三对角方程组的追赶法一、算法理论在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对角线方程组11112222211111n n n n n n n n n b c x f a b c x f a b c x f a b x f -----⎛⎫⎛⎫⎛⎫⎪⎪ ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪= ⎪⎪⎪ ⎪⎪⎪ ⎪⎪ ⎪⎝⎭⎝⎭⎝⎭, 简记为Ax f =。
求解Ax f =:等价于解两个三角形方程组,Ly f y =求;,Ux y x =求.从而得到解三对角线方程组的追赶法公式:(1)计算{}i β的递推公式()111/,/,2,3,,1;i i i i i c b c b a i n βββ==-=- (2) 解Ly f =()()11111/,/,2,3,,;i i i i i i i y f b y f a y b a i n β--==--=(3)解Ux y =1,,1,2,2,1.n n i i i i x y x y x i n n β+==-=--我们将计算系数12112n n y y y βββ-→→→→→→ 及的过程称为追的过程, 将计算方程组的解11n n x x x -→→→ 的过程称为赶的过程。
#include <stdio.h>#include <math.h>#include<stdlib.h>#define N 20double a[N], b[N], c[N-1], f[N], r[N];int n;int i;void LUDecompose(); // LU分解void backSubs(); // 回代void main(){printf("请输入方程的维数n=");scanf("%d",&n);getchar();if(n>N||n<=0){printf("由于该维数过于犀利, 导致程序退出!");return;}printf("\n输入下三角元素\n");printf("输入%d个a值: ", n-1);for (i=1; i<n; i++)scanf("%lf", &a[i]);getchar();printf("\n输入主对角线元素\n");printf("输入%d个b值: ", n);for (i=0; i<n; i++)scanf("%lf", &b[i]);getchar();printf("\n输入上三角元素\n");printf("输入%d个c值: ", n-1);for (i=0; i<n-1; i++)scanf("%lf", &c[i]);getchar();printf("\n输入%d个方程组右端项: \n", n);for (i=0; i<n; i++)scanf("%lf", &f[i]);getchar();LUDecompose();backSubs();printf("\n线性方程组的解为: \n");for (i=0; i<n; i++)printf("x%d=%lf\n", i+1, f[i]);}void LUDecompose(){ //α被b取代, β被c取代, 以节省存储空间c[0]=c[0]/b[0];for(i=1;i<n-1;i++){r[i]=a[i];b[i]=b[i]-r[i]*c[i-1];c[i]=c[i]/b[i];}r[i]=a[i];b[i]=b[i]-r[i]*c[i-1];}void backSubs(){ // y被f取代, x也被f取代, 以节省存储空间f[0]=f[0]/b[0];for(i=1; i<n; i++)f[i]=(f[i]-r[i]*f[i-1])/b[i];f[n-1]=f[n-1];for(i=n-2;i>=0;i--)f[i]=f[i]-c[i]*f[i+1];}四、 算法实现例1.用该程序计算三对角线方程组2100012100A 012100012100012⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭---=-----, 10000b ⎛⎫⎪ ⎪⎪= ⎪ ⎪ ⎪⎝⎭计算其方程组的解。
数值分析追赶法实验报告(3篇)
第1篇一、实验目的通过本次实验,掌握追赶法的基本原理和计算步骤,了解追赶法在解三对角线性方程组中的应用,并学会利用C++编程实现追赶法,提高编程能力。
二、实验原理追赶法是一种解三对角线性方程组的迭代方法,其基本原理是利用递推公式逐步求解未知数。
对于形如Ax=b的三对角线性方程组,其中系数矩阵A具有如下形式:A = [a_00, a_01, 0, ..., 0;a_10, a_11, a_12, ..., 0;0, a_21, a_22, ..., a_2n-1;...;0, ..., 0, a_n2]追赶法将系数矩阵A分解为两个因子L和U,其中L为下三角矩阵,U为上三角矩阵,即:A = LU其中:L = [1, 0, ..., 0;a_10/a_11, 1, 0, ..., 0;..., ..., ..., ...;a_n1/a_n2, ..., ..., ..., 1]U = [a_11, a_12, ..., a_1n;0, a_22, ..., a_2n-1;..., ..., ..., ...;0, ..., ..., a_n2]通过递推公式求解L和U中的元素,进而得到解向量x:x_1 = b_1 / a_11x_i = (b_i - ∑(j=1 to i-1) l_ij x_j) / u_ij, i = 2, ..., n三、实验步骤1. 编写C++程序,实现追赶法的基本算法。
2. 生成三对角线性方程组的系数矩阵A和解向量b。
3. 调用C++程序,计算追赶法的结果,并输出解向量x。
4. 分析追赶法的计算过程,验证结果是否正确。
四、实验数据及结果1. 生成三对角线性方程组的系数矩阵A和解向量b。
假设A为:A = [2, -1, 0, 0, 0;-1, 2, -1, 0, 0;0, -1, 2, -1, 0;0, 0, -1, 2, -1;0, 0, 0, -1, 2]b = [1, 1, 1, 1, 1]2. 追赶法计算结果。
追赶法(Thomas算法)
二、解三对角线性方程组的追赶法 定理1:满足引理1条件的三对角方阵A有如下形式的 唯一的克劳特分解。
p1 a2 A= pn
p2 a3 pn 1 an
1 q1 1 q2 =PQ 1 qn 1 1
其中
p1 = b1 (i = 1,2,, n 1) qi = ci pi p = b a q (i = 2,3,, n ) i i i 1 i
解三对角线方程组Ax = f可化为求解两个三角形 方程组
Py = f
Qx = y
(1) 解 Py = f
p1 a2 ( P, f ) = p2 a3 pn 1 an f1 f2 f3 pn f n
得
{
y1 = f1 / p1
yi = ( f i ai yi 1 ) / piຫໍສະໝຸດ (i = 2,3,, n )
( 2) 解 Qx = y
1 q1 1 q2 1 qn 1 1
x1 y1 x2 = y2 x y n n
得
xn = y n
xi = yi qi xi +1
i = n 1 , , 2 ,1
作业: P50 习题11
§2-4
追赶法(Thomas算法 算法) 追赶法 算法
一、对角占优矩阵
若矩阵A = ( aij )n× n 满足
|aii |> ∑|aij |
j =1 j ≠i ≠i
n
i = 1 , 2 , , n
则称A为严格对角占优矩阵.
若矩阵A = ( aij )n× n 满足
|aii | ∑|aij | ≥
j =1 j ≠i
n
i = 1 , 2 , , n
则称A为弱对角占优矩阵.
matlab追赶法解101阶三对角方程组
在探讨MATLAB追赶法解101阶三对角方程组之前,我们首先需要了解什么是追赶法和什么是三对角方程组。
追赶法又称托马斯算法,是一种用于求解带状矩阵(即只有主对角线和两条相邻的对角线上有非零元素的矩阵)的线性方程组的方法。
而三对角矩阵就是只有主对角线和两条相邻的对角线上有非零元素的矩阵。
在实际应用中,求解带状矩阵的线性方程组是非常常见的,特别是在数值计算和科学工程领域。
现在,让我们深入探讨MATLAB追赶法解101阶三对角方程组的方法和具体步骤。
一、MATLAB追赶法解101阶三对角方程组1. 概念介绍101阶三对角方程组是一个非常大的线性方程组,通常使用传统的高斯消元法来求解会耗费大量的时间和计算资源。
而MATLAB追赶法通过利用三对角矩阵的特殊性质,可以有效地简化计算过程,并且节省大量的内存和计算资源。
2. 追赶法步骤(1)将原方程组化为追赶法所需的形式;(2)利用追赶法求解三对角线性方程组。
二、追赶法求解101阶三对角方程组的实现过程1. 将原方程组化为追赶法所需的形式对于101阶三对角方程组,我们首先需要将其化为追赶法所需的形式。
这个过程涉及到选取合适的追赶元和追赶子以及对原方程组的变形,将其化为追赶法能够直接处理的形式。
2. 利用追赶法求解线性方程组一旦将原方程组化为追赶法所需的形式,我们就可以利用追赶法对其进行求解。
追赶法的核心是通过追赶子的迭代计算,逐步求得线性方程组的解。
在MATLAB中,可以使用内置的追赶法求解函数,也可以编写自定义的追赶法算法来实现对101阶三对角方程组的求解。
三、个人观点和理解在实际工程和科学计算中,追赶法是一种非常有效的求解带状矩阵线性方程组的方法。
对于大规模的三对角方程组,特别是高阶的情况,传统的直接求解方法往往会遇到内存和计算资源的限制,而追赶法能够通过精巧的迭代计算,在保证解的精度的显著提高计算效率。
在MATLAB中,通过调用内置的追赶法函数,可以快速地求解大规模的三对角方程组,极大地方便了工程实践中的数值计算工作。
基于高斯消去法的追赶法
基于高斯消去法的追赶法
基于高斯消去法的追赶法是一种求解三对角矩阵线性方程组的方法。
三对角矩阵的特点是除了主对角线上的元素外,只有两个相邻的副对角线上有非零元素。
追赶法的思想是通过对三对角矩阵进行分解,将其转化为两个上、下三角矩阵的乘积形式,从而简化方程组的求解过程。
具体步骤如下:
1. 首先,对于三对角矩阵A的第一行进行变换,使得A的第一个元素为1。
2. 然后,对于i = 2, 3, ..., n,进行追赶过程,将A的第i 行的第i-1个元素(副对角线上)消去,将A的第i-1行的第i个元素(副对角线上)消去,并更新A的第i行的第i个元素(主对角线上)。
3. 最后,通过回代法求解得到方程组的解。
追赶法的优点是简单、高效,适用于解决大规模的三对角线性方程组。
它在科学计算和工程领域中有广泛的应用,例如求解抛物型偏微分方程、求解材料传输过程等。
追赶法求解三对角线性方程组
追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容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 ,令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 。
追赶法求解三对角方程及其算例
追赶法求解三对角方程组要求:对于给定的三对角系数矩阵和右端项,可以求解线性代数方程组一、 追赶法的数学理论设系数矩阵为三对角矩阵112223311100000000000000n 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 βγββαβγγβαβγ--+⎧===⎪⎪=⎪⎪⎪==-=⎪⎪⎨-⎪=⎪⎪=⎪⎪=--⎪=-⎪⎩对对(*)二、 追赶法的算法和流程图1.预处理生成方程组的系数i u 及其除数i d ,事实上,按式(*)可交替生成i d 与i u :1d →1u →2d →…→1-n u →n d其计算公式为⎪⎩⎪⎨⎧-=-===+++1,...,2,1,,/c u b 111i i 11n i u a b d d d i i i i i 2.追的过程顺序生成方程组右端:i y →2y →…→n y据式(*)的计算公式为n i d y a f y d f y i i i i i ,...,3,2,/)(/1111=⎩⎨⎧-==-逆序得出方程组的解i x :n x →1-n x →…→1x其计算公式按式为1,2,1,1,⋯--=⎩⎨⎧-==+n n i x u y x y x i i i in n 三、 追赶法的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); endy(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四、 追赶法的算例实现算例 用追赶法求解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡1016432153-1-21-2-31-1-2x x x x 解答令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。
解三对角方程组的追赶法
~ A LU
~ 角矩阵,且 A L DP
, unn ), P D1U , 则 P 为单位上三
A
对称,故有
A P DL A LDP
T T T
15
第二章 解线性方程组的直接方法
AT PT DLT A LDP
由Doolitle分解的唯一性,得
P L, 即
T
A P DP
T
对任意非零向量 是正定的
T
x,
T
y P1x ,
T 1 1
于是由 A
x Dx x ( P ) AP x T 1 T 1 ( P x) AP x y Ay 0
即 D 正定,所以 D 的对角元素均为正数。令
~ D diag( u11 ,, u nn )
16
第二章 解线性方程组的直接方法
(li )
(ui )
di ai di 1 di
( yi )
12
第二章 解线性方程组的直接方法
3. dn bn dn
4. 对
( xn )
i n 1, n 2,,1 (d i ci d i 1 ) bi d i
5.输出
dx
( xi )
,停机.
追赶法的基本思想与Gauss消去法及三角分解法相同, 只是由于系数中出现了大量的零,计算中可将它们撇开, 从而使得计算公式简化,也大大地减少计算量.其乘除运 算量为 5n 4
ck bk 1 c k 1
0 c n 1 a n bn
8
第二章 解线性方程组的直接方法
其中 uk bk ck 1 ak uk 1
(k 2,, n)
数值分析实验报告之追赶法求三对角矩阵
xi yi ui xi1 , i n 1, ,1 ,显示求解结果,过程结束。
实 验 总 结
本次试验也是针对求解线性方程组解的问题。当一个矩阵的结构比较特殊时,利用 特殊的方法进行对待,这样可以提高解题的速度,本次实验所针对的一类矩阵是比较特 殊的, 因此用了追赶法进行求解。 因为对实验的核心算法有所了解, 在编程实现的时候, 基本上没有遇到算法的逻辑错误。在网上找到了几组实验数据,我分别对它们进行了测 试,在确定结果无误之后,然后求解老师实验所给数据。本次实验使我对线性方程组的 求解有了进一步的认识,增加了一种方程组的求解方法,对更加深入的学习线性方程组 有很大的帮助。
u1 ci y yi 1 c1 y , y1 1 , ui , yi i , i 2, , n 1, b1 b1 bi ui 1ai bi ui 1ai
yn
yn yn1an ,然后求解 xi yi ui xi 1 , i n 1, ,1 。 bn un1an
实验原理
b1 c1 x1 y1 x1 u1 x2 y1 x u x y a2 b2 c2 x2 y 2 2 2 3 2 化为 x an1 bn1 cn 1 xn1 y n 1 u n 1 xn y n 1 n 1 an bn xn y n xn y n
实验步骤
Step5:
方程组,过程结束,否则转到 Step5:; 2i;
Step6: 计算; ui ci / bi ui1ai , yi ( yi yi1 ) / bi ui1ai Step7: 判断 i 是否大于 n ,若是,转到 Step8,否则,令 r 1 r ,返回 Step6。 Step8:
用追赶法求解三对角方程组
用追赶法求解三对角方程组1. 三对角方程组的背景大家好,今天咱们来聊聊一个有点学术味儿的话题——三对角方程组。
不过别担心,我会尽量让这件事情变得轻松有趣,就像跟朋友聊天一样。
三对角方程组呢,其实就是那些系数在对角线附近的线性方程组,听起来是不是有点复杂?别急,咱们慢慢来,打个比方,它就像是一个田字格,只在主要的对角线上有数字,其他地方都是零。
哎,生活中有很多时候我们会遇到这样的方程,比如在物理、工程或者计算机科学里。
这时候,咱们就得想办法求解它们。
2. 追赶法的简介2.1 追赶法是什么好啦,接下来咱们来介绍一下追赶法。
这法子听上去是不是有点像小时候玩捉迷藏的感觉?其实它就是一种巧妙的迭代算法,专门用来解决那些三对角的线性方程组。
为什么叫追赶法呢?因为它能快速“追赶”到正确的解,就像小兔子在草地上跑得飞快一样。
它的基本思路就是把这个三对角方程组转化为一个更简单的形式,从而一步一步找到答案。
2.2 为什么用追赶法那为什么不直接用其他的方法呢?哦,朋友们,真相是,追赶法在处理这类方程的时候特别高效,速度快得像闪电!想象一下,如果你在一场马拉松里,你会选择走路还是飞奔?当然是飞奔啦!同样的道理,追赶法能节省大量的计算资源和时间,让我们轻松愉快地拿到想要的解。
3. 追赶法的步骤3.1 初始准备咱们要开始追赶了,首先得准备一下。
你需要把方程组写成标准的形式,通常我们可以把它表示成一个矩阵。
这样一来,咱们就能更清晰地看到那些三角形的结构。
接着,得设定好初始条件,这就好比你出发前检查好背包里有没有水、食物和地图。
没有这些东西,你可不敢贸然出门啊!3.2 逐步追赶准备好之后,追赶法就开始工作了。
第一步,咱们需要对三角形的每一行进行“消元”,也就是让下面的元素逐渐变为零。
听起来是不是有点复杂?其实就像在厨房里切菜,先把最上面的部分处理掉,然后逐步往下进行。
一步一步来,绝对不能急,这样才能确保每一刀都精准无误。
接着,咱们要开始反向代入,也就是从最后一行开始,逐行算出未知数。
MATLAB-追赶法求解三对角方程组的算法原理例题与程序
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。
追赶法解三角矩阵
cout<<"请输入三对角阵的规模:"<<endl;
cin>>n;
cout<<"请输入三对角阵:"<<endl;
for(i=0;i<n;i++){ //输入a[][]
for(s=0;s<n;s++){
cin>>a[i][s];
}
cout<<endl;
}
cout<<"请输入b的值:"<<endl;
u[i][i+1]=a[i][i+1];
l[i][i]=1.0;
l[i+1][i]=a[i+1][i]/u[i][i];
u[i+1][i+1]=a[i+1][i+1]-l[i+1][i]*a[i][i+1];
}
cout<<"根据Doolittle分解:"<<endl;
cout<<"L矩阵为:"<<endl;
*/
cout<<u[i][s]<<" \t ";
}
cout<<endl;
}
y[0]=b[0]; //解Ly=b(追过程)
for(i=1;i<n;i++){
y[i]=b[i]-l[i][i-1]*y[i-1];
}
x[n-1]=y[n-1]/u[n-1][n-1]; //解Ux=y(赶过程)
matlab追赶法求解三对角方程组
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。
追赶法
解:对于对角占优的三对角线性矩阵可用追赶法求解 采用 Fortran 编程计算如下:
program main implicit none integer::i,j,in,nb real::x real::a(1:100),b(1:100),c(1:100),f(1:100) real::d(1:100),y(1:100),v(1:100) open(unit=100,file='a.dat') open(unit=101,file='b.dat') open(unit=102,file='c.dat') open(unit=103,file='f.dat') open(unit=104,file='v.dat') do while(.not.eof(100)) read(unit=100,fmt=*)in,x a(in)=x end do nb=0 do while(.not.eof(101)) read(unit=101,fmt=*)in,x b(in)=x nb=nb+1 end do do while(.not.eof(102)) read(unit=102,fmt=*)in,x c(in)=x end do do while(.not.eof(103))
read(unit=103,fmt=*)in,x f(in)=x end do d(1)=c(1)/b(1) do i=2,nb-1 d(i)=c(i)/(b(i)-a(i)*d(i-1)) end do y(1)=f(1)/b(1) do i=2,nb y(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*d(i-1)) end do v(nb)=y(nb) do i=nb-1,1,-1 v(i)=y(i)-d(i)*v(i+1) end do write(unit=104,fmt=201)'y(i)',',','x(i)' do i=1,nb write(unit=104,fmt=200)y(i),',',v(i) end do 200 format(f7.4,a1,f7.4) 201 format(a7,a1、c.dat 和 f.dat 的内容分别如图 1、图 2、图 3 和图 4 所示;计算结果保存在 v.dat 中,如图 5 所示:
追赶法(Thomas算法)
定理1:满足引理1条件的三对角方阵A有如下形式的 唯一的克劳特分解。
p1 a2 p2
1 q1
1 q2A Nhomakorabea
a3 pn1
an
pn
=PQ
1 qn1
1
其中
qpi1
b1 ci
pi
i 1,2,, n 1
得
xi yi qi xi1
i n 1,,2,1
作业:
P50 习题11
感谢您的下载让小编的努力能帮助到您, 最后一页是小编对你的谢谢哦,提醒一下, 下载好了几个全部自己看一遍,把用不上 的删除哦!包括最后一页!
yi fi ai yi1 / pi i 2,3,, n
(2) 解Qx y
1
q1 1
q2
1
qn1
1
x1 x2
y1 y2
xn yn
xn yn
pi
bi
aiqi1
i 2,3,, n
解三对角线方程组 Ax f可化为求解两个三角形 方程组
Py f Qx y
(1) 解Py f
p1
f1
a2 p2
f2
(P, f )
a3
f3
pn1 an
pn
f
n
得
y1 f1 / p1
追赶法解三对角矩阵
实验追赶法解三对角方程组一、实验目的学会用追赶法解三对角方程组,并应用该算法于实际问题.二、实验要求给定三对角方程组,应用追赶法解得方程组的解。
三、实验内容1、追赶法2、以课本数值试验2为实例3、如果有错,修改直至运行成功,查看运行结果;四、实验环境matlab五、实验步骤和方法1、程序设计2、带入实例3、撰写实验报告。
六、实验预习要求得到实例的解一、[源程序]function x = my_zgf2(A,d,flag)%MY_ZGF2 Summary of this function goes here[m,n]=size(A); %计算矩阵的大小if nargin==2; %输入变量等于2的时候,A中储存所有元素的值for i=1:na(i)=A(i+1,i);b(i)=A(i,i);c(i)=A(i,i+1);enda(1)=0; %补充不足的值b(n)=A(n,n);c(n)=0;elsec=[A(1,:) 0]; %flag==1时b=A(2,:);a=[0 A(3,:)];endu(1)=b(1);for i=2:n %第一次追赶,得到上、下三角矩阵l(i)=a(i)/u(i-1);u(i)=b(i)-c(i-1)*l(i);endy(1)=d(1); %解Ly=dfor i=2:ny(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n); %解Ux=yfor i=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1))/u(i);end二、带入实例A =-2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 05.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000-2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 0d= 8.1400 0 0 0 0 0 0 0>> d=A(4,:);my_zgf2(A,d,1)ans =2.0350 1.0174 0.5086 0.2541 0.1267 0.0626 0.0298 0.0119 >>。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验追赶法解三对角方程组
一、实验目的
学会用追赶法解三对角方程组,并应用该算法于实际问题.
二、实验要求
给定三对角方程组,应用追赶法解得方程组的解。
三、实验内容
1、追赶法
2、以课本数值试验2为实例
3、如果有错,修改直至运行成功,查看运行结果;
四、实验环境
matlab
五、实验步骤和方法
1、程序设计
2、带入实例
3、撰写实验报告。
六、实验预习要求
得到实例的解
一、[源程序]
function x = my_zgf2(A,d,flag)
%MY_ZGF2 Summary of this function goes here
[m,n]=size(A); %计算矩阵的大小
if nargin==2; %输入变量等于2的时候,A中储存所有元素的值for i=1:n
a(i)=A(i+1,i);
b(i)=A(i,i);
c(i)=A(i,i+1);
end
a(1)=0; %补充不足的值
b(n)=A(n,n);
c(n)=0;
else
c=[A(1,:) 0]; %flag==1时
b=A(2,:);
a=[0 A(3,:)];
end
u(1)=b(1);
for i=2:n %第一次追赶,得到上、下三角矩阵
l(i)=a(i)/u(i-1);
u(i)=b(i)-c(i-1)*l(i);
end
y(1)=d(1); %解Ly=d
for i=2:n
y(i)=d(i)-l(i)*y(i-1);
end
x(n)=y(n)/u(n); %解Ux=y
for i=n-1:-1:1
x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
二、带入实例
A =
-2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 0
5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000
-2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 0
d= 8.1400 0 0 0 0 0 0 0
>> d=A(4,:);
my_zgf2(A,d,1)
ans =
2.0350 1.0174 0.5086 0.2541 0.1267 0.0626 0.0298 0.0119 >>。