计算方法(B)实验报告
姓名:
学号:
学院:
专业:
实验一 三对角方程组Tx f =的求解
一、 实验目的
掌握三对角方程组Tx f =求解的方法。 二、 实验内容
求三对角方程组Tx f =的解,其中:
4 -1 -1 4 -1 -1 4 1 -1 4T ????????=??
??
???? , 3223f ?? ? ? ?= ? ? ???
三、 算法组织
设系数矩阵为三对角矩阵
11222333111
b c a b c a b c a b c b n n n n T ---????????
=??????
??????
则方程组Tx f =称为三对角方程组。
设矩阵T 非奇异,T 可分解为T=LU ,其中L 为下三角矩阵,U 为单位上三角矩阵,记
1
1
212
313
1
1
1111
,11n n n n n r l r l r L U l r l μμμμμ---????
? ? ? ? ? ?==
?
? ? ? ? ? ? ? ? ??
?
?
?
可先依次求出,L U 中的元素后,令Ux y =,先求解下三角方程组Ly f =得出
y ,再求解上三角方程组Ux y =。 追赶法的算法组织如下:
1.输入三对角矩阵T 和右端向量f ;
2.将Tx f =压缩为四个一维数组{}{}{}{}i i i i a b c d 、、、,{}{}{}i i i a b c 、、是T 的三对角线性方程组的三个对角,{}i d 是右端向量。将分解矩阵压缩为三个一维数组
{}{}{}i i i l r μ、、。
3.对T 做Crout 分解(也可以用Doolittle 分解)导出追赶法的计算步骤如下:
1111,b r c μ==
for 2i n =
111, ,
,i i i i i i i i i i i i i l a b a r r c y d l y μμ---==-==-
end
4.回代求解x
/n n n x y μ=
for 11i n =-
1()/i i i i i x y c x μ+=-
end
5. 停止,输出结果。 四、 MATLAB 程序 MATLAB 程序见附件1. 五、 结果及分析 实验结果为:
(1.0000 1.0000
1.0000 1.0000)T x =
实验二 Jacobi 迭代和Gauss-Seidel 迭代解线性方程组
一、 实验目的
掌握Jacobi 迭代和Gauss-Seidel 迭代解线性方程组的方法。 二、 实验内容
用Jacobi 迭代和Gauss-Seidel 迭代解电路电流方程组,使各部分电流误差均小于310-。
2
12132343452328103831001025150154503050i i i i i i i i i i i i -=??-+-=??
-+-=??-+=?
-+=?? 三、 算法组织
形如Ax b =的方程组,用Jacobi 迭代求解x ,算法组织如下:
1. 将系数矩阵A 分解成对角元D 、下三角部分元E -和上三角部分元F -,于是
A D E F =--.
2. 由11()()()Ax b D E F x b Dx E F x b x D E F x D b --=?--=?=++?=++。
3. 从而构成形如(1)()k k x Gx d +=+迭代格式:
(1)1()1()k k x D E F x D b +--=++
其中
11
()
G D E F d D b
--?=+?=? 4. 选取初始向量(0)x 进行迭代计算。
5. 当迭代后的解满足题中的约束条件(1)()1max k k i i i n
x x ε+≤≤-<时迭代停止。
形如Ax b =的方程组,用Gauss —Seidel 迭代求解x ,算法组织如下: 1. 将系数矩阵A 分解成对角元D 、下三角部分元E -和上三角部分元F -,于是
A D E F =--.
2. 由
11()()()()Ax b D E F x b D E x Fx b x D E Fx D E b --=?--=?-=+?=-+-。 3. 从而构成形如(1)()k k x Gx d +=+迭代格式:
(1)1()1()()k k x D E Fx D E b +--=-+- 其中
11
()()G D E F
d D E b
--?=-?=-? 4. 选取初始向量(0)x 进行迭代计算。
5. 当迭代后的解满足题中的约束条件(1)()1max k k i i i n
x x ε+≤≤-<时迭代停止。
四、 MATLAB 程序
MATLAB 程序见附件2,其中1为Jacobi 迭代,2为Gauss —Seidel 迭代。 五、 结果及分析 Jacobi 迭代结果:
方程组的解为(0.36070.03350.01630.00540.0055)T x = 迭代次数8i =
Gauss —Seidel 迭代结果:
方程组的解为(0.36070.03350.01660.00550.0056)T x = 迭代次数4i =
由以上结果可知,达到相同的计算精度,Gauss —Seidel 迭代比Jacobi 迭代的速度快,Gauss —Seidel 迭代比Jacobi 迭代次数少。
实验三 多项式插值及误差计算
一、实验目的
掌握多项式插值的原理和基本方法。 二、实验内容
已知2
1
()(11)125f x x x
=-≤≤+,对5,10,20n = a. 计算函数()f x 在点2
1,(0,1,2,,)i x i i n n
=-+=处的值()i f x ;
b. 求插值数据点(){}
,(0,1,2,,)i i x y i n =的Newton 插值多项式()n N x 和三
次样条插值多项式()n S x ; c. 对5,20n =,计算2
1,(110,9099)100
k x k k =-+
=和相应的函数值()
,(),()k n k n k k y f x N x S x =;
d. 计算()()
max n n k k k
E N y N x =-,()()
max n n k k k
E S y S x =-,解释所得到
结果。
三、 算法组织
(一)本题第一问是简单的用matlab 程序可以计算,算法很简单。
(二) 本题在算法上第二问中的Newton 插值多项式)(x N n 和三次样条插值多项式()n S x 。计算两种插值多项式的算法如下: 1. 求Newton 插值多项式
)
(x N n ,算法组织如下:
Newton 插值多项式的表达式如下:
010011()()()()()n n n N x c c x x c x x x x x x -=+-+???+--???-
其中每一项的系数c i 的表达式如下:
12011010
(,,,)(,,,)
(,,,)i i i i i f x x x f x x x c f x x x x x -???-???=???=
-
根据上述公式,为了得到系数需计算:
1) 一阶差商01[],[],,[]n f x f x f x ??? 2) 二阶差商01121[,],[,],[,]n n f x x f x x f x x -??? … … … …
3) n 阶差商01111[,,,],[,,,]n n n f x x x f x x x --?????? 4) n+1阶差商011[,,,,]n n f x x x x -???
2. 求三次样条插值多项式,算法组织如下:
所谓三次样条插值多项式()n S x 是一种对区间进行分段的分段函数,然后在每一段上进行分析,即它在节点i x 011()n n a x x x x b -=<??<<=分成的每个小区间
1[,]i i x x -上是3次多项式,其在此区间上的表达式如下:
2233
1111111()[()()]()()666[,]1,2,,.
i i i i i i i i i i i i i i i
i i h x x h x x S x x x M x x M y M y M h h h x x x i n --------=-+-+-+-∈=???,
, 因此,只要确定了i M 的值,就确定了整个表达式,i M 的计算方法如下: 令:
111
11111116()6(,,)i i i i i i i i i i i i i i
i i i i i i i h h h h h h y y y y d f x x x h h h h μλμ++++--+++?
===-?++??
--?=-=?+?
, 则i M 满足如下n-1个方程:
1121,2,,1i i i i i i M M M d i n μλ-+++==???-,
方程中有n+1个未知量,则令0M 和n M 分别为零,则由上面的方程组可得到(11)i M i n ≤≤-的值,可得到整个区间上的三次样条插值多项式()n S x 。 (三) 第三问和第四问的算法与第二问的算法类似,不再赘述。
四、 MATLAB 程序 MATLAB 程序见附件3。 五、 结果及分析 第一问
当n=5时,各节点及f(x)值为: x(0)=-1,y(0)=3.846154e-02
x(1)=-6.000000e-01,y(1)=1.000000e-01 x(2)=-2.000000e-01,y(2)=5.000000e-01 x(3)=2.000000e-01,y(3)=5.000000e-01 x(4)=6.000000e-01,y(4)=1.000000e-01 x(5)=1,y(5)=3.846154e-02
当n=10时,各节点及f(x)值为:
x(0)=-1,y(0)=3.846154e-02
x(1)=-8.000000e-01,y(1)=5.882353e-02 x(2)=-6.000000e-01,y(2)=1.000000e-01 x(3)=-4.000000e-01,y(3)=2.000000e-01 x(4)=-2.000000e-01,y(4)=5.000000e-01 x(5)=0,y(5)=1
x(6)=2.000000e-01,y(6)=5.000000e-01 x(7)=4.000000e-01,y(7)=2.000000e-01 x(8)=6.000000e-01,y(8)=1.000000e-01 x(9)=8.000000e-01,y(9)=5.882353e-02 x(10)=1,y(10)=3.846154e-02
当n=20时,各节点及f(x)值为:
x(0)=-1,y(0)=3.846154e-02
x(1)=-9.000000e-01,y(1)=4.705882e-02 x(2)=-8.000000e-01,y(2)=5.882353e-02 x(3)=-7.000000e-01,y(3)=7.547170e-02 x(4)=-6.000000e-01,y(4)=1.000000e-01 x(5)=-5.000000e-01,y(5)=1.379310e-01 x(6)=-4.000000e-01,y(6)=2.000000e-01 x(7)=-3.000000e-01,y(7)=3.076923e-01