(西南交大)数值分析上机作业.docx

数值分析作业

专业班级：14级隧道1班
姓 名： 戴龙钦
学 号： 2014200095
指导老师： 赵海良

2014年12月 成都

序 言
通过数值分析的理论知识的学习，此次实验将我们学过的理论知识运用于实践之中。本次实验，我选用的计算机语言为MATLAB，其主要有一下几个特点。
1．编程效率高
MATLAB是一种面向科学与工程计算的高级语言，允许使用数学形式的语言编写程序，且比BASIC、FORTRAN和C等语言更加接近我们书写计算公式的思维方式，用MATLAB编写程序犹如在演算纸上排列出公式与求解问题。因此，MATLAB语言也可通俗地称为演算纸式科学算法语言。由于它编写简单，所以编程效率高，易学易懂
2. 用户使用方便
MATLAB语言与其他语言相比，较好的解决了上述问题，把编辑、编译、链接和执行融为一体。它能在同一画面上进行灵活操作，快速排除输入程序中的书写错误、语法错误以至语义错误，从而加快了用户编写、修改和调试程序的速度，可以说在编程和调试过程中它是一种比VB还要简单的语言。
3. 方便的绘图功能
MATLAB的绘图是十分方便的，它有一系列绘图函数（命令），例如线性坐标、对数坐标、半对数坐标及极坐标，均只需调用不同的绘图函数（命令），在图上标出图题、XY轴标注，格（栅）绘制也只需调用相应的命令，简单易行。另外，在调用绘图函数时调整自变量可绘出不变颜色的点、线、复线或多重线。这种为科学研究着想的设计是通用的编程语言所不能及的。

I

目 录
1 必做题 1
1.1 题目1 1
1.1.1 题目 1
1.1.2 计算原理 1
1.1.2.1 多项式拟合 1
1.1.2.2 插值多项式 3
1.1.3 计算过程及结果 4
1.1.3.1 多项式拟合 4
1.1.3.2 插值多项式 4
1.1.4 结果分析 4
1.2 题目3 4
1.2.1 题目 4
1.2.2 计算原理 5
1.2.2.1 Jacobi迭代法 5
1.2.2.2 Gauss-Seidel迭代法 6
1.2.3 计算过程及结果 7
1.2.4 结果分析 8
2 选做题 10
2.1 题目2 10
2.1.1 题目 10
2.1.2 计算原理 10
2.1.2.1 Euler算法 10
2.1.2.2 4阶Runge-Kutta算法 11
2.1.2.3 方法的稳定性 12
2.1.3 计算过程及结果 13
2.1.4 结果分析 14
3总结 15
4附件 16
4.1 必做第1题程序 16
附录1.1 16
附录1.2 16
4.2 必做第1题结果 16
附录2.1 16
附录2.2 16
4.3 必做第3题程序 17
附录1.3 17
附录1.4 17
附录1.5 18
附录1.6 19
附录1.7 19
附录1.8 19
附录1.9 20
附录1.10 20
附录1.11 20
附录1.12 21
附录1.13 21
附录1.14 21
4.4 必做第3题结果 22
附录2.3 22
附录2.4 22
附录2.5 22
附录2.6 23
附录2.7 23
附录2.8 23
附录2.9 24
附录2.10 24
附录2.11 24
附录2.12 25
4.5 选做第2题程序 25
附录1.15 25
附录1.16 26
附录1.17 26
附录1.18 26
附录1.19 27
附录1.20 27
4.6 选做第2题结果 27
附录2.13 27
附录2.14 27
附录2.15 27
附录2.16 27
附录2.17 28
附录2.18 28
附录2.19 28
附录2.20 28

II
2014级学硕隧道1班 戴龙钦 2014200095
1 必做题
1.1 题目1
1.1.1 题目
某过程涉及两变量x 和y，拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi与yi之间的对应数据如下，xi=1,2,…,10
yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570
24.8234 75.2795 103.5743 97.4847 78.2392
（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。
（2）请用插值多项式给出最好近似结果
下列数据为另外的对照记录，它们可以作为近似函数的评价参考数据。
xi =
Columns 1 through 7
1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000
Columns 8 through 14
4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000
Columns 15 through 17
7.1000 7.5000 7.9000
yi =
Columns 1 through 7
42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006
Columns 8 through 14
-18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840
Columns 15 through 17
79.5688 93.7700 102.3677
1.1.2 计算原理
1.1.2.1 多项式拟合
假设给定数据点 (i=0,1,…,m)，为所有次数不超过的多项式构成的函数类，现求一，使得
(1-1)
当拟合函数为多项式时，称为多项式拟合，满足式（1-1）的称为最小二乘拟合多项式。特别地，当n=1时，称为线性拟合或直线拟合。
显然

为的多元函数，因此上述问题即为求的极值 问题。由多元函数求极值的必要条件，得
(1-2)
即
 (1-3)
（3）是关于的线性方程组，用矩阵表示为
(1-4)
式（1-3）或式（1-4）称为正规方程组或法方程组。
可以证明，方程组（1-4）的系数矩阵是一个对称正定矩阵，故存在唯一解。从式（1-4）中解出 (k=0,1,…，n)，从而可得多项式
 (1-5)
可以证明，式（1-5）中的满足式（1-1），即为所求的拟合多项式。我们把称为最小二乘拟合多项式的平方误差，记作

由式(1-2)可得
(1-6)
多项式拟合的一般方法可归纳为以下几步：
1) 由已知数据画出函数粗略的图形——散点图，确定拟合多项式的次数n；
2) 列表计算和；
3) 写出正规方程组，求出；
4) 写出拟合多项式。
在实际应用中，或；当时所得的拟合多项式就是拉格朗日或牛顿插值多项式。 
1.1.2.2 插值多项式
在众多函数中,多项式最简单、最易计算，已知函数 在n+1个互不相同的点处的函数值，为求的近似式，自然应当选n次多项式
(1-7)
使满足条件
， (1-8)
称为被插函数，城插值多项式，式（1-8）称为插值条件，称插值节点。这种求函数近似式的方法称为插值法。几何上，其实质是用通过n+1个点的多项式由曲线，当做曲线的近似曲线，如图1-1所示。

图1-1 插值多项式的几何意义
1.1.3 计算过程及结果
1.1.3.1 多项式拟合
设拟合模型为

利用matlab编程，程序见附录1.1
最终计算机运行结果见附录2.1。
由于多项式拟合中，多项式的阶数越大，拟合的精度越高，但是在超过7阶会产生龙格现象，所以，6次多项式为最好的近似结果为：

1.1.3.2 插值多项式
由题意知，x，y的十个已知点，现给出另外17个x值，利用编程求出对应这17个点的y值。本文采用三次样条插值法，Matlab编程见附录1.2
最终计算结果见附录2.2。
通过与题中所给数据对比来看，所得结果计算差别不大。
1.1.4 结果分析
通过与题中所给数据对比来看，两种方法所得结果计算差别不大。
1.2 题目3
1.2.1 题目
用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b1或Ax=b2，研究其收敛性。上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。
(1) A行分别为：
A1=[6,2,-1]，A2=[1,4,-2]，A3=[-3,1,4]；b1=[-3,2,4]T；b2=[100,-200,345]T
(2) A行分别为：
A1=[1,0,8,0.8]，A2=[0.8,1,0.8]，A3=[0.8,0.8,1]；b1=[3,2,1] T；b2=[5,0,-10]T
(3) A行分别为：
A1=[1,3]，A2=[-7,1]；b1=[4,6]T
1.2.2 计算原理
1.2.2.1 Jacobi迭代法
设有n阶方程组Ax=b，若系数矩阵非奇异，且 (i = 1, 2,…, n)，将方程组改写成同解方程组：

然后写成迭代格式：

上式也可以简单地写为：

对以上两式给定一组任意初值后，经反复迭代可得到一向量序列，如果x (k)收敛于，则就是方程组Ax=b的解，该方法称为雅克比(Jacobi)迭代法。程序实现雅克比(Jacobi)迭代法据体流程如图1-1所示。
设D = diag (a11, a22, …, ann)，将Ax=b改写为： AX = (D – (D - A)) x = b，DX = (D - A) x + b，X = (I – D-1A) x + D-1b。记B = I – D-1A，g = D-1 b。则迭代格式的向量表示为

BJ = I – D-1A称为雅克比迭代矩阵。
Jacobi迭代收敛的充要条件：迭代矩阵BJ的谱半径。特别地，若系数矩阵A满足以下特征中的任意一条，则Jacobi迭代法收敛：
A为行对角占优阵，即；
A为列对角占优阵，即；
A的元素满足。

图1-1 Jacobi迭代法流程图 图1-2 Gauss-Seidel迭代法流程图
1.2.2.2 Gauss-Seidel迭代法
在雅可比迭代中，每次迭代时只用到前一次的迭代值，而在高斯-塞德尔迭代中，每次迭代时充分利用最新的迭代值。一旦一个分量的新值被求出，就立即用于后续分量的迭代计算，而不必等到所有分量的新值被求出以后。如果迭代收敛，应该比更接近于原方程的解(i = 1, 2,…, n)，因此在迭代过程中及时地以代替（i = 1, 2,…, n-1），可得到更快的收敛效果。这样可将迭代格式写成：

上式可简写为：
(i = 1, 2,…, n)
该式称为Gauss-Seidel迭代格式。程序实现Gauss-Seidel迭代法据体流程如图1-2所示。
对于上述Gauss-Seidel迭代式，如写成矩阵形式为：
= D-1+ D-1b, =。
其中，L和U分别为：

则Gauss-Seidel迭代法的迭代矩阵为。
其收敛的充要条件为谱半径。特别地，若系数矩阵A满足Jacobi迭代法三条特征中的任意一条，则Gauss-Seidel迭代法收敛。
1.2.3 计算过程及结果
当b1=时，用雅阁比法结合matlab编程计算。首先要编写雅阁比与高斯-赛德尔迭代法的自定义函数，然后通过调用此函数来解题。
雅阁比自定义函数程序见附录1.3。高斯-赛德尔迭代法自定义函数见附录1.4。完成后，把以上程序保存成m文件，以备调用。
现对本题中的三个小题分别利用雅阁比法进行计算, 选取x(0) =，迭代10 次。精度选。
(1) 首先用雅阁比法进行计算，当b1=时，程序见附录1.5。当b2=时，程序见附录1.6。通过计算机计算，当为b1时，见附录2.3，当为b2时，见附录2.4。
现在用高斯-赛德尔法进行计算：当b1=时，程序见附录1.7，当b2=时，程序见附录1.8。通过计算，当b1=时，计算结果见附录2.5，当为b2时，所得结果见附录2.6。
通过以上对比，为表达清楚，现将总结如下：
用雅阁比法计算b1、b2结果如下：
=[-0.726719 0.809460 0.251958]
=[36.426331 -2.015201 114.049096]
用高斯-赛德尔法计算结果如下：
=[-0.727247 0.808074 0.252546]
=[36.363673 -2.075136 114.041539]
(2) 仿照（1）中的计算步骤，先对其进行雅阁比法计算，b1=时，程序见附录1.9，b2=时，程序见附录1.10。运算结果，当为b1时，见附录2.7，当为b2时，见附录2.8。
现对其进行高斯-赛德尔迭代法计算，当b1=[3,2,1]时，程序见附录1.11，当b2=[5,0,-10]时，程序见附录1.12。运算结果，当为b1时，见附录2.9，当为b2时，见附录2.10。
(3) 同样，由（1）题的计算方法，先用雅阁比进行计算，程序见附录1.13，运行结果见附录2.11。再用高斯-赛德尔迭代法进行计算，程序见附录1.14，运行结果见附录2.12。
1.2.4 结果分析
(1) 通过以上计算分析，由结果来看，他们的收敛性都比较好，理论分析正确，但是从过程来看，通过对b1、b2本身的对比，高斯-赛德尔法的收敛速度要快于雅阁比法，基本上在第五六次迭代过程中，高斯-赛德尔法就已经很接近真实值，而雅阁比法在第七八次迭代过程中才基本接近真实值。将b1与b2进行对比后得出，在相同计算方法的前提下，右端项对迭代收敛有一定的影响，b2的收敛性要差于b1的收敛性。
(2) 由以上方法得出，当用雅阁比法进行计算时，无论取b1还是b2，都没有收敛性，所得数据过于离散。雅阁比法在本题失效（原因是本题给出的行列式无法满足雅阁比迭代收敛的条件：A既不是行对角占优矩阵也不是列对角占优矩阵，同时）。但是用高斯-赛德尔法可以取得较好的收敛性。等式右端项对迭代收敛性没有什么明显的影响。
(3) 由此可得，对于本小题而言，无论是雅阁比法还是高斯-赛德尔法都已经失效，无法得到收敛的数列。因为此题给出的行列式无法满足雅阁比法与高斯-赛德尔迭代法所给出的收敛条件。A既不是行对角占优矩阵，也不是列对角占优矩阵，且。不满足其充要条件。所以不能用这两种方法计算。

2 选做题
2.1 题目2
2.1.1 题目
给定初值问题y/=-1000(y-x2)+2x , 0£x£1。y(0)=1，请分别用Euler和4阶Runge-Kutta算法按步长h=0.1,0.01,0.001,0.0001计算其数值解，分析其中遇到的现象及问题。
2.1.2 计算原理
2.1.2.1 Euler算法
一阶常微分方程初值问题
（2-1）
第一步，对式2-1的求解区域作一致网格剖分，得到一致网格节点。

第二步，对连续方程作节点离散，得到节点离散方程。
（2-2）
第三步，由泰勒公式

改写上式为

同样地，有

由离散节点方程得
（2-3）
忽略高阶项, 并用作为真解的近似值，得差分方程：

即
（2-4）
总的，
（2-5）
开始
输入x0,y0,h,n
X1=x0+h
y1=y0+h*f(x0,y0);
i=1
输出x1,y1
i=n?
结束
i=i+1
x0=x1,y0=y1
否
图2-1 Euler算法流程图
这就是著名的欧拉(Euler)公式。
以上方法称为欧拉方法或欧拉折线法。程序实现雅克比(Jacobi)迭代法据体流程如图1-1所示。
2.1.2.2 4阶Runge-Kutta算法
一阶常微分方程初值问题如式2-1所示，常微分方程初值问题的数值解法是求方程（2-1）的解在点列上的近似值，这里是到的步长，一般略去下标记为。
常微分方程初值问题的数值解法一般分为两大类：
（1）单步法：这类方法在计算时，只用到、和，即前一步的值。因此，在有了初值以后就可以逐步往下计算。典型方法如龙格–库塔方法。
（2）多步法：这类方法在计算时，除用到、和以外，还要用，即前面步的值。典型方法如Adams方法。
经典的方法是一个四阶的方法，它的计算公式是：
（2.6）
开始
输入x0,y0,h,N
x1=x0+h
k1=f(x0,y0),k2=f(x0+h/2,y0+hk1/2)
k3=f(x0+h/2,y0+hk2/2),k4=f(x1,y0+hk3)
y1=y0+h(k1+2k2+2k3+k4)/6
i=1
输出x1,y1
i=n?
结束
i=i+1
x0=x1,y0=y1
否
是
图2-2 4阶Runge-Kutta算法流程图
方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值。程序实现Runge-Kutta算法据体流程如图2-2所示。
2.1.2.3 方法的稳定性
粗略的说：如果是某方法第k步的精确值，为其近似值，其绝对误差为，即。假设第k步后的计算中不再有舍入误差，只是由引起的扰动（，），都有，则称此方法是绝对稳定的。
因此，对于此题，根据绝对稳定的定义，对于Euler算法设有一扰动，此时

故。
如果，则必有，即时，Euler法是绝对稳定的，称为Euler法的绝对稳定区间。同理，可得4阶Runge-Kutta算法的绝对稳定区间为。
2.1.3 计算过程及结果
用matlab编程计算方程初值Euler和4阶Runge-Kutta算法时，首先要编写Euler和4阶Runge-Kutta算法的自定义函数，然后通过调用此函数来解题。
Euler算法自定义函数程序见附录1.15。4阶Runge-Kutta算法自定义函数见附录1.16。完成后，把以上程序保存成m文件，以备调用。
现对本题中按步长h=0.1,0.01,0.001,0.0001计算其数值解分别利用Euler和4阶Runge-Kutta算法进行计算。
根据式y/=-1000(y-x2)+2x，以及初始条件y(0)=1，可以得到实际的方程，取绝对误差小于最为计算停止条件。
(1) 当h=0.1时，a=0，b=1，y(0)=1，即分段数n=10时，Euler算法程序见附录1.17。4阶Runge-Kutta算法程序见附录1.17。通过计算机计算，Euler算法结果见附录2.13，4阶Runge-Kutta算法结果见附录2.14。
通过对两种算法结果对比，为表达清楚，现将结果列出如下：
用Euler算法计算结果如下：
i
x(i)
y(i)
yx(i)
wucha(i)
0
0.0000
1.0000
1.0000
0.0000
1
0.1000
-99.0000
0.0100
99.0100
2
0.2000
9802.0200
0.0400
9801.9800
3
0.3000
-970395.9400
0.0900
970396.0300
用4阶Runge-Kutta算法计算结果如下：
i
x(i)
y(i)
yx(i)
wucha(i)
0
0.0000
1.0000
1.0000
0.0000
1
0.1000
16174067.2633
0.0100
16174067.2533
(2) 当h=0.01时，即分段数n=100时，Euler算法程序见附录1.18。4阶Runge-Kutta算法程序见附录1.17。通过计算机计算，Euler算法结果见附录2.15，4阶Runge-Kutta算法结果见附录2.16。
(3) 当h=0.001时，即分段数n=1000时，Euler算法程序见附录1.19。4阶Runge-Kutta算法程序见附录1.17。通过计算机计算，Euler算法结果见附录2.17，4阶Runge-Kutta算法结果见附录2.18。
(4) 当h=0.0001时，即分段数n=10000时，Euler算法程序见附录1.20。4阶Runge-Kutta算法程序见附录1.19。通过计算机计算，Euler算法结果见附录2.13，4阶Runge-Kutta算法结果见附录2.20。
2.1.4 结果分析
(1) 当h=0.1时，h在绝对稳定区间之外，故结果应该是不稳定的，从结果可以看见，不论Euler还是4阶Runge-Kutta算法，结果失去稳定性，而且还可以看出，收敛阶数越大，结果在是稳定情况下，精度越差。
(2) 当h=0.01时，h在绝对稳定区间之外，结果失稳，从附录2.15和2.16可以看出，结果失稳。
(3) 当h=0.001与0.0001时，h在绝对稳定区间之内，两种算法均收敛，而且从附录2.17与2.18的对比，附录2.19与2.20的对比可以看出，4阶Runge-Kutta算法精度高于Euler算法。但是并不是步长h越大越好，从结果可以看出，h=0.001与h=0.0001在达到相同的精度时，h=0.0001的计算步远大于h=0.001的计算步。

3 总结
本报告系西南交通大学2014级硕士研究生《数值分析》课程的上机实习报告，由本人严格按照实习要求独立完成。完成必做第一题（插值多项式和多项式拟合）、必做第三题（雅格比迭代法、高斯—赛德尔迭代法求解方程组的问题）和选做第二题（利用四阶Runge-Kutta算法求解微分方程的初值问题）。我分别利用C语言和Matlab完成了程序的编写，并从中感受到了两种编程语言的特点。
本次实习总共经历了六个步骤，分别是：选题；学习相关知识（尤其是C语言基础知识的回顾）；分析解题思路；编写程序代码；验证实验数据；编写实习报告。通过使用C语言和Matlab编写数值计算的程序，让我加深了对数学编程软件的使用，同时也使我对各种算法的原理有了更深一层的理解。
这是我第一次独立使用数学编程工具，在编程之前首先务必要详细了解其基本原理，然后设计好算法流程，最终完成源程序的编写。由于平时接触的很少，在编写程序的时候往往显得力不从心，很多东西都需要重新学习并向师兄师姐请教，在编写过程中也遇到了许多困难。通过本次实习，算是对C语言和MATLAB进行了一次全面而系统的复习，也对其有了进一步的深入了解，尤其是对其基本结构、优特点和使用方法。这其中的苦与乐是值得回味的，这个过程也是让我受益匪浅的。
通过此次数值分析上机实习，认识到了计算机编程在数值分析课程中的重要性。也让我意识到了学习编程语言的重要性。要想要将数学应用于实际工程中，特别是对于工科的学生来讲，这是一门非常重要的实践课程。
在此次报告中，首次接触了Matlab这门软件，之所以选这门，是因为很多工程中对数据的处理需要使用，这对我本身的专业是非常重要的。
本学期学习的数值分析方法，让我在对数学的理解上又有了新的认识。在一连串看似杂乱无章的数据中，用数值分析进行处理，有了很多的收获。在学习的过程中，意识到了这门课对我专业的重要性，让我对自己的专业有了更好的发挥。对我写论文以及解决实际工程问题有很大的帮助。
最后，感谢老师给我这样的机会去接触这门语言，虽然只了解了皮毛，可是仍然收获颇多。由于初次接触这门软件，在报告中仍难免会有不完善甚至错误的地方，望谅解！
4 附件
4.1 必做第1题程序
附录1.1
X=[1 2 3 4 5 6 7 8 9 10]
Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392]
P3=polyfit(X, Y,3)
P4=polyfit(X, Y,4)
P5=polyfit(X, Y,5)
P6=polyfit(X, Y,6)
附录1.2
X=[1 2 3 4 5 6 7 8 9 10]
Y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392]
Xi=[1.5 1.9 2.3 2.7 3.1 3.5 3.9 4.3 4.7 5.1 5.5 5.9 6.3 6.7 7.1 7.5 7.9]
Yi= interp1(X, Y, Xi,'spline')
4.2 必做第1题结果
附录2.1
P3 = -1.0326 19.3339 -94.4787 131.7944
P4 = -0.3818 7.3680 -42.1433 73.5334 0.7450
P5 = 0.0981 -3.0789 34.5020 -163.5107 304.7282 -139.5019
P6 = 0.0194 -0.5408 5.1137 -16.8973 -0.8670 66.3750 -18.6991
附录2.2
Yi = Columns 1 through 11
43.0510 41.6393 34.7861 24.1356 11.3345 -1.6605 -12.2916 -18.0792 -17.8286 -11.0129 2.1296
Columns 12 through 17
19.8846 40.3192 61.0680 79.5479 93.6798 102.3680
4.3 必做第3题程序
附录1.3
function tx=jacobi(A,b,imax,x0,tol)
del=10^-6;
tx=[x0] ; n=length(x0);
for i=1:n
dg=A(i,i);
if abs(dg)< del
disp('diagonal element is too small');
return
end
end
for k = 1:imax
for i = 1:n
sm=b(i) ;
for j = 1:n
if j~=i
sm = sm -A(i,j)*x0(j) ;
end
end %for j
x(i)=sm/A(i,i) ;
end
tx=[tx ;x] ;
if norm(x-x0)<tol
return
else
x0=x ;
end
end
附录1.4
function tx= gseidel( A,b,imax,x0,tol)
del=10^-6;
tx=[x0]; n=length(x0);
for i=1:n
dg=A(i,i);
if abs(dg)< del
disp('diagonal element is too small');
return
end
end
for k = 1:imax
x=x0;
for i = 1:n
sm=b(i);
for j = 1:n
if j~=i
sm = sm -A(i,j)*x(j);
end
end
x(i)=sm/A(i,i);
end
tx=[tx;x];
if norm(x-x0)<tol
return
else
x0=x;
end
end
附录1.5
A=[6 2 -1;1 4 -2;-3 1 4];
b=[-3 2 4];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx=jacobi(A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.6
A=[6 2 -1;1 4 -2;-3 1 4];
b=[100,-200,345];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx=jacobi(A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.7
A=[6 2 -1;1 4 -2;-3 1 4];
b=[-3 2 4];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx= gseidel (A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.8
A=[6 2 -1;1 4 -2;-3 1 4];
b=[100,-200,345];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx= gseidel (A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.9
A=[1 0.8 0.8;0.8 1 0.8;0.8 0.8 1];
b=[3 2 1];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx=jacobi(A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.10
A=[1 0.8 0.8;0.8 1 0.8;0.8 0.8 1];
b=[5 0 -10];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx=jacobi(A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.11
A=[1 0.8 0.8;0.8 1 0.8;0.8 0.8 1];
b=[3 2 1];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx= gseidel (A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.12
A=[1 0.8 0.8;0.8 1 0.8;0.8 0.8 1];
b=[5 0 -10];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,3);
tx= gseidel (A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f %f\n', j, tx(j,1),tx(j,2),tx(j,3))
end
附录1.13
A=[1 3;-7 1];
b=[4 6];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,2);
tx=jacobi(A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f\n', j, tx(j,1),tx(j,2))
end
附录1.14
A=[1 3;-7 1];
b=[4 6];
tol=1.0*10^-6 ;
imax =10;
x0= zeros(1,2);
tx= gseidel (A,b,imax,x0,tol) ;
for j=1:size(tx,1)
fprintf('%d %f %f\n', j, tx(j,1),tx(j,2))
end
4.4 必做第3题结果
附录2.3
迭代次数i
x1
x2
x3
1
0.000000
0.000000
0.000000
2
-0.500000
0.500000
1.000000
3
-0.500000
1.125000
0.500000
4
-0.791667
0.875000
0.343750
5
-0.734375
0.869792
0.187500
6
-0.758681
0.777344
0.231771
7
-0.720486
0.805556
0.236654
8
-0.729076
0.798448
0.258247
9
-0.723108
0.811392
0.253581
10
-0.728201
0.807567
0.254821
11
-0.726719
0.809460
0.251958
附录2.4
迭代次数i
x1
x2
x3
1
0.000000
0.000000
0.000000
2
16.666667
-50.000000
86.250000
3
47.708333
-11.041667
111.250000
4
38.888889
-6.302083
124.791667
5
39.565972
2.673611
116.992187
6
35.274161
-1.395399
115.256076
7
36.341146
-1.190502
113.054470
8
35.905912
-2.558051
113.803485
9
36.486598
-2.074736
113.818947
10