数值计算方法上机实验报告.docx
《数值计算方法上机实验报告.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实验报告.docx（40页珍藏版）》请在冰点文库上搜索。
数值计算方法上机实验报告
华
北
电
力
大
学
科
技
学
院
实
验报告
课程名称：
数值计算方法上机实验报告
姓名：
牛玺童班级：
电气11k6
学号：
111904010415
1.题目
实验一造倒数表
造倒数表，并例求18的倒数。
（精度为0.0005）
2.算法原理
kk
2.1牛顿迭代法牛顿迭代法是通过非线性方程线性化得到迭代序列的一种方法。
对于非线性方程f（x）=0，若已知根x*的一个近似值x，将f（x）在x处展
成一阶泰勒公式后忽略高次项可得：
f（x）≈
f（xk）+f'（xk）（x-xk）
右端是直线方程，用这个直线方程来近似非线性方程f（x）。
将非线性方程
f（x）=0的根x*代入f（x*）=0，即
*
f（xk）+f'（xk）（x-xk）≈0
k
解出x*≈x-
f（xk）
f'（xk）
*
将右端取为xk+1，则xk+1是比xk更接近于x的近似值，即
xk+1
≈xk-
f（xk）
f'（xk）
这就是牛顿迭代公式，相应的迭代函数是
ϕ（x）=x-
f（x）f'（x）
2.2牛顿迭代法的应用
计算1是求cx-1=0的解，解出x，即得到1。
取
f（x）=cx-1，f'（x）=c，
c
有牛顿迭代公式
xk+1
=xk
c
-cxk-1=1
cc
这样就失去了迭代的意义，达不到迭代的效果。
故重新构造方程：
cx2-x=0，1也是该式的解。
故取f（x）=cx2-x，
c
f'（x）=2cx-1，则有牛顿迭代公式
kkk
cx2-xcx2
xk+1=xk-
2cx
=,
-12c-1
k=0,1,...
1的值在1
kk
~1之间，取初值x=0.1。
18
3.流程图
20100
4.输出结果
5.结果分析
当k=3时，得5位有效数字0.05564。
此时，x3-x4=0.00000<0.0005，
3
故取x*=x=0.05564≈0.056。
此种迭代格式仍存在一定的缺陷，经实验后发现当初值x0
>x*时必收敛，但
是当x0
0）时迭代结果发散，ς较小尚不确定。
6.心得体会
起初对题目的理解并不是很透彻，另外对构建牛顿迭代公式理论依据不是特别充分，比如说为什么在原有直接得到的式子两边各乘一个x，只是试出来的。
在编程方面不够成熟。
当然也加深了对牛顿迭代法的理解和应用的具体实现。
1.题目
实验二例3-4
用列主元消去法求解方程组
⎧
12x1-3x2+3x3=
15⎫
⎪⎪
⎨-18x1-3x2-x3=-15⎬
⎪⎪
⎩x1+
x2+
x3=6⎭
并求出系数矩阵A的行列式的值detA。
2.算法原理
2.1顺序高斯消去法
顺序高斯消去法是利用线性方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加至另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。
这样，顺序高斯消去法可分成“消去”和“回代”两个过程。
在用顺序高斯消去法时，在消元之前检查方程组的系数矩阵的顺序主子式，当阶数较高时是很难做到的。
若线性方程组的系数具有某种性质时，如常遇到的对角占优方程组，自然能够用高斯消去法求解。
2.2列选主元消去法
kk
线性方程组只要系数矩阵非奇异，就存在惟一解，但是按顺序消元过程中可能出现主元素a（k）=0，这时尽管系数矩阵非奇异，消元过程无法再进行，或者
即使（k）0
akk≠，但如果其绝对值很小，用它作除数也会导致其他元素的数量级急剧
增大和使舍入误差扩大，将严重影响计算的精度。
为避免在校园过程确定乘数时的所用除数是零或绝对值小的数，即零主元或
小主元，在每一次消元之前，要增加一个选主元的过程，将绝对值大的元素交换到主对角线的位置上来。
列选主元是当高斯消元到第k步时，从k列的akk以下（包括akk）的各元素中选出绝对值最大的，然后通过行交换将其交换到akk的位置上。
交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。
列选主元消去法常用来求行列式。
设有矩阵
⎛a11
La1n⎞
A=⎜MM⎟
⎜⎟
a
⎜n1Lann⎟
⎝⎠
112233
用列主元消去法将其化为上三角形矩阵，对角线上元素为a
（1）,a
（2）,L,a（3）,于是行
列式
detA=（-1）m
aaLa
（1）
（2）（n）
1122nn
其中m为所进行的行交换次数。
这是实际中求行列式值的可靠方法。
3.流程图
d?
d
≠
4.输出结果
5.结果分析
采用计算机运算在计算大数据时有明显的优点，另外也需要考虑到存储。
kk
高斯消去法的使用条件是a（k）≠0,
k=1,2,L,n，而列选主元法可以保证这一
条件。
并且可以避免在消元过程确定乘数时所用除数是绝对值小的数，相对全选主元的运算量小，一般也可以满足精度要求。
6.心得体会
此次上机不仅需要对原理了解透彻，而且要求的编程能力较强。
在定义和思路上没问题，只是在编程软件的使用上遇到些不稳定的问题，如头文件的使用。
在存储空间上得到了新的认识，另外发现了当代码多时流程框图的好处。
编程是一件很需要耐心的事，自己还有很大进步空间。
1.题目
实验三例3-10
用杜里特尔分解法求矩阵A的逆矩阵A-1。
⎛11
A=⎜12
-1⎞
-2⎟
2.算法原理
⎜⎟
⎝⎠
⎜-211⎟
2.1杜里特尔分解法
（1）
（1）
（2）
设线性方程组Ax=b，对系数矩阵A进行除不交换两行位置得初等行变换相当于用初等矩阵M1左乘A，在对方程组第一次消元后，A和b分别化为A和
b
（2），即
⎧MA
（1）=A
（2）
⎪1
⎨
Mb
（1）
=b
（2）
⎪⎩1
⎧1⎫
⎪⎪
⎪⎪-m211⎪⎪
其中M1=⎨-m311⎬
⎪⎪
⎪MO⎪
⎪⎩-mn11⎭⎪
第k次消元时，A（k）和b（k）分别化为A（k+1）和b（k+1），即
⎧MA（k）=A（k+1）
⎪k
⎨
Mb（k）=
b（k+1）
⎩⎪k
⎧1⎫
⎪⎪
⎪
M=
⎪
其中1⎨
⎪
O⎪
1⎪
⎬
-mk+1,kO⎪
⎪MO⎪
⎪⎪
⎩⎪-mnk
1⎭⎪
消元过程是对k=1~n-1进行的，因此有
⎧⎪M
LMMA
（1）=A（n）
⎨
⎪⎩M
（1）
n-121
b（n）
将上三角形矩阵A（n）记为U，于是有
A=M-1-1
LMn-1U=LU
⎛1⎞
⎜⎟
其中L=M-1M-1
M-1
⎜m21
⎜m
=31
1⎟
m321⎟
12L
n-1
⎜
⎜MM
⎟
m43O⎟
⎜MMMO⎟
⎜⎜⎟⎟
为单位下三角形矩形。
⎝mn1
mn2
mn3L1⎠
这样高斯消去法的实质是将系数矩阵A分解为两个三角形矩阵L和U相乘，即A=LU
在上述矩阵描述中遇到了下三角形矩阵运算。
主对角线以上元素全为零的方阵称为下三角形矩阵。
下三角形矩阵的乘积仍是下三角形矩阵。
若下三角形矩阵可逆，其逆矩阵仍是下三角形矩阵，而且下三角形矩阵的乘积和逆矩阵很容易求得。
把A分解成一个单位下三角阵和一个上三角阵U的乘积成为杜里特尔分解。
这种分解是惟一的。
2.2高斯-约当法
高斯消去法有消元和回代两个过程，当对消元过程稍加改变便可以使方程组化为对角形方程组
Dx=b
的形式，其中矩阵D为对角形矩阵，即
（1）
11
⎜
D=⎜
⎜
（2）⎟
22⎟
O⎟
⎜⎟
⎜a（n）⎟
⎝nn⎠
当高斯-约当消去法消元的每一步都先用主元去除其所在行的各元素（包括常数项）时，方程组便可化成
（n）
⎛1⎞⎛x1⎞
⎛b1⎞
⎜⎟⎜⎟⎜⎟
1xb（n）
⎜⎟⎜
2⎟=⎜2⎟
⎜O⎟⎜M⎟
⎜M⎟
⎜⎟⎜⎟⎜⎟
⎝1⎠⎝xn⎠
⎜b（n）⎟
⎝n⎠
这是等号右端即为方程组的解。
高斯-约当消去法每一步都用主元去除其所在行
的各元素（包括常数项），这个个过程成为归一化，这时方程组的系数阵转化为单位阵。
为减小误差，高斯-约当消去法还常用列选主元技术。
3.流程图
4.输出结果
5.结果分析
采用杜里特尔分解法求解方程组时，由于把对系数矩阵的计算和对右端项的计算分开，这使计算线性方程组系非常方便。
只需进行一次矩形三角分解，然后再解多个三个方程组，且多解一个方程组仅需要增加大约n2次乘除法运算。
采用高斯约当法仅需要进行消元归一，而不需要回代，为编程实现提供了便利。
6.心得体会
步骤很重要，审题--确定算法--解题步骤--流程图--程序--代入简单值进行验证。
在编程时先在代码输入区打好框架，并且尽量在每一命令后注释，方便检查错误和日后复习。
定义和变量存储很灵活，如我把单位向量直接赋给了A矩阵变量中，还有根据最终的目的直接简化计算。
另外赋值前，确定存储空间并且要定义初值为零。
实验四例4-6
1.题目
已知f（x）的观测数据
x
1
2
3
4
f（x）
0
-5
-6
3
构造插值多项式。
2.算法原理
首先构造基函数lk
（x）=
n
Õ
i=0
i≠k
x-xixk-xi
，可以证明基函数满足下列条件：
⎧0
l（x）=⎨
i≠k
，
ki⎩1
i=k
对于给定n+1个节点，n次拉格朗日插值多项式由下式给出：
n
L（x）=∑lk（x）yk
k=0
其中lk
（x）=
n
Õ
i=0
i≠k
x-xixk-xi
由于lk（x）是一个关于x的n次多项式，所以L（x）为关于x的不高于n次的代数多项式。
当x=xi时，L（xi）=yi，满足插值条件。
3.流程图
4.输出结果
5.结果分析
由于所知的拉格朗日计算机算法只能实现计算某一特定值的近似函数值，而不知如何导出表达式，故例求x=2.5处的函数值以说明表达式以得出，只是在计算机程序中。
并且也能达到拉格朗日插值法使用的目的。
6.心得体会
编程不够细心，程序没问题，却因为不知道是输入文件错了而检查了好长时间。
但同时也加深了对拉格朗日基函数性质的认识和理解。
1.题目
实验五习题5-2
给出平面函数z（x,y）=ax+by+c的数据
i
1
2
3
4
5
xi
0.1
0.2
0.4
0.6
0.9
yi
0.2
0.3
0.5
0.7
0.8
zi
0.58
0.63
0.73
0.83
0.92
按最小二乘原理确定a，b，c。
2.算法原理
2.1最小二乘原理
设已知某物理过程y=f（x）的一组观测数据
（xi,f（xi））,i=1,2,L,m
要求在某特定函数类Φ（x）中寻找一个函数ϕ（x）作为y=f（x）的近似函数，使得
二者在xi上的误差或称残差
δi=ϕ（xi）-f（xi）,i=1,2,L,m
按某种度量标准为最小，这就是拟合问题。
要求残差δi按某种度量标准为最小，即要求由残差δi构成的残差向量
δ=[δ,δ,L,δ]T
的某种范数δ为最小。
例如，要求δ
，或δ即
1∞
mm
δ1=∑δi
i=0
=∑
i=0
ϕ（xi）-f（xi）
δ=maxδ
∞
=max
ϕ（xi）-f（xi）
ii
为最小，这本来都是很自然的，可是计算不太方便。
所以通常要求：
11
⎛m2⎞2⎧m2⎫2
δ2=⎜∑δi⎟
=⎨∑[ϕ（xi）-f（xi）]⎬
⎝i=0⎠⎩i=0⎭
或者
δ2=
mm
δ2=
[ϕ（x）-f（x）]2
2∑i
i=0
∑ii
i=0
为最小。
这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法。
2.2多变量数据拟合
对于给定的一组数据（xi,yi），i=1，2…，m，寻求做n次多项式
使性能指标
n
k
y=∑
k=0
akx
mn
01n
∑i∑ki
J（a,a,L,a）=
（y
i=1
-
k=0
axk）2为最小。
由于性能指标J可以被看做关于ak，k=0，1，…，n的多元函数，故上述拟合多项式的构造问题可转化为多元函数的极值问题。
令
从而有正则方程组
∂J=0
∂ak
∑i∑i
⎢
∑i⎥a0⎢∑i⎥
⎡mxx2L
xn⎤⎡⎤⎡y⎤
⎢xx2
3n+1⎢⎥
xLx⎥a⎢
xy⎥
∑i∑i∑i
⎢
∑i⎥⎢
1⎥=⎢∑ii⎥
⎢MMMLM
⎢M⎥M
⎢⎥
⎢xn
xn+1
xn+2
x2n
⎥⎣a⎦⎢
xny⎥
⎣∑i∑i
∑iL
∑i⎦n
⎣∑ii⎦
对多变量（或称多元）线性模型
01122nn
y*=a+ax+ax+L+ax
进行了m次观测
⎧y*=a
+ax
+ax
+L+ax
10111221
⎪
n1n
⎪y*=a
+
ax
+
ax
+L+ax
⎨20112222
n2n
⎪M
⎪y*=a
+ax
+ax
+L+ax
⎩n011n
22nnnn
这个称为回归方程组，写成矩阵形式
y=Aα
1
⎡y*⎤
⎡1x11
x21L
xn1⎤
⎡a0⎤
⎢⎥⎢⎥⎢⎥
y*1xxLxa
其中y=⎢
2⎥,A=⎢
1222
n2⎥,α=⎢
1⎥。
⎢M⎥
⎢MMMMM⎥
⎢M⎥
⎢⎥⎢⎥⎢⎥
*1xx
Lxa
⎢⎣ym⎥⎦
⎣1m
2mnm⎦⎣n⎦
当m>n时，要确定一组ai,i=0,1,L,n,使之精确地满足m个方程，这是超定方程组的问题，只能在最小平方误差的基础上确定αi。
12m
定义残差向量δ=[δ,δ,L,δ]T,则
δ=y-Aα
12m
其中y=[y,y,L,y]T
代替输出向量。
取性能指标
J=δTδ
使之最小，以此确定出α。
由
J=δTδ=（y-Aα）T（y-Aα）
=yTy-αTAy-yTAα+αTATAα
利用向量和矩阵的运算公式，有
ATAα=ATy
此即为正则方程组，当ATA非奇异时，可求得
α=（ATA）-1ATy
3.流程图
4.输出结果
5.结果分析
曲线拟合的最小二乘法是反映所给数据点的总的趋势，并不是严格的通过每个数据点，这样就避免了大量数据插值时需要高次多项式，同时又去掉了数据所含的测量误差。
6.心得体会
整理思路越来越熟练了，所以执行各个步骤也相对简单了很多。
另外对原理也加深了认识。
附录：
1.造倒数表
1>源程序：
#include"stdio.h"
#include"math.h"
#defineN30voidmain（）
{
inti;
floatx[N],c;FILE*fp1,*fp2;
fp1=fopen（"input1.txt","r"）;
fp2=fopen（"output1.txt","w"）;
fscanf（fp1,"%f",&c）;
fscanf（fp1,"%f",&x[0]）;//初值fprintf（fp2,"****倒数表****\n"）;
for（i=0;i{
x[i+1]=x[i]*x[i]*c/（2*c*x[i]-1）;
fprintf（fp2,"k=%d\tx（%d）=%.5f\n",i,i,x[i]）;
if（fabs（x[i+1]-x[i]）<=0.0005）break;
elsecontinue;
}
fprintf（fp2,"k=%d\tx（%d）=%.5f\n",i+1,i+1,x[i]）;fprintf（fp2,"\n计算结果：
\n1/%f=%.3f\n\n",c,x[i+1]）;
fclose（fp1）;fclose（fp2）;
}
2>输入文件：
"input1.txt"18
0.1
3>输出文件：
“output1.txt”
****倒数表****k=0x（0）=0.10000k=1x
（1）=0.06923k=2x
（2）=0.05781k=3x（3）=0.05564k=4x（4）=0.05564
计算结果：
1/18.000000=0.056
2.例3-4
1>源程序：
#include"iostream"
#include"cmath"usingnamespacestd;
#defineN10voidmain（）
{
inti,j,k,l,n;
floatb[N],a[N][N],t,d,det=1.0;
//***数据输入*/
FILE*fp1,*fp2;fp1=fopen（"input2.txt","r"）;
fp2=fopen（"output2.txt","w"）;
fscanf（fp1,"%d",&n）;for（i=0;i//***数据输入*/
//*****************************************************************************
*************高斯消去*/
//*****************************************消元*/
//****************************列选主元函数*/for（k=0;k{
d=a[k][k];l=k;for（i=k+1;i{
if（fabs（a[i][k]）>fabs（d））
{d=a[i][k];i=l;}
}
if（i==n）//判断是否奇异，不奇异进行行交换
{
if（d==0）
fprintf（fp2,"奇异"）;//如果所有行的“首列”都为0，为奇异else
{
if（l!
=k）//如果第k行的“首列”并不是最大
{
det=det*（-1）;for（j=k;j<=n;j++）//交换系数矩阵中的两行
{t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}t=b[l];b[l]=b[k];b[k]=t;//交换右端常向量中的两行
}
}
}
//****************************列选主元函数*/
for（i=k+1;i{
a[i][k]=a[i][k]/a[k][k];for（j=k+1;j{
a[i][j]=a[i][j]-a[i][k]*a[k][j];
}
b[i]=b[i]-a[i][k]*b[k];//右端常向量
}
}
//*****************************************消元*/
//*回代*/
b[n-1]=b[n-1]/a[n-1][n-1];//计算x（N）的解
for（i=n-2;i>=0;i--）//从倒数第二项开始依次回代N-1次
{t=0;
for（j=i+1;j{t=t+a[i][j]*b[j];}
b[i]=（b[i]-t）/a[i][i];
}
//*****************************************************************************
*************高斯消去*/
//************************************数据输出*/for（i=0;ifor（i=0;ifprintf（fp2,"
升级会员 