数值计算方法上机实验报告Word格式文档下载.docx
- 文档编号:3872207
- 上传时间:2023-05-02
- 格式:DOCX
- 页数:30
- 大小:121.49KB
数值计算方法上机实验报告Word格式文档下载.docx
《数值计算方法上机实验报告Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《数值计算方法上机实验报告Word格式文档下载.docx(30页珍藏版)》请在冰点文库上搜索。
4.输出结果
丈锹F}精⑹«
5C(OJ查看㈣WW{H)
——kkkkk
**ajc*怖]数表**址*
K(0)=0.10000
3t{l)=o.06923
X
(2)=0.05781
x{3)=0.05564
x{4)=0.05564
计篦结果=
1/18.000000^0.056
5.结果分析
当k3时,得5位有效数字0.05564此时,X3X40.000000.0005,
故取X*X30.055640.056。
此种迭代格式仍存在一定的缺陷,经实验后发现当初值X0X时必收敛,但
是当X0X*(0)时迭代结果发散,较小尚不确定。
6.心得体会
起初对题目的理解并不是很透彻,另外对构建牛顿迭代公式理论依据不是特别充分,比如说为什么在原有直接得到的式子两边各乘一个X,只是试出来的。
在编程方面不够成熟。
当然也加深了对牛顿迭代法的理解和应用的具体实现。
用列主元消去法求解方程组
12x1
3x2
3x3
15
18x1
x3
x1
x2
6
并求出系数矩阵A的行列式的值detA。
2.1顺序高斯消去法
顺序高斯消去法是利用线性方程组初等变换中的一种变换,即用一个不为零
的数乘一个方程后加至另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
这样,顺序高斯消去法可分成“消去”和“回代”两个过程。
在用顺序高斯消去法时,在消元之前检查方程组的系数矩阵的顺序主子式,当阶数较高时是很难做到的。
若线性方程组的系数具有某种性质时,如常遇到的对角占优方程组,自然能够用高斯消去法求解。
2.2列选主元消去法
线性方程组只要系数矩阵非奇异,就存在惟一解,但是按顺序消元过程中可
能出现主元素a(kkk)0,这时尽管系数矩阵非奇异,消元过程无法再进行,或者即使akk(k)0,但如果其绝对值很小,用它作除数也会导致其他元素的数量级急剧
增大和使舍入误差扩大,将严重影响计算的精度。
为避免在校园过程确定乘数时的所用除数是零或绝对值小的数,即零主元或
小主元,在每一次消元之前,要增加一个选主元的过程,将绝对值大的元素交换到主对角线的位置上来。
列选主元是当高斯消元到第k步时,从k列的akk以下(包括akk)的各元素中选出绝对值最大的,然后通过行交换将其交换到akk的位置上。
交换系数矩阵中的两行(包括常数项),只相当于两个方程的位置交换了,因此,列选主元不影响求解的结果。
列选主元消去法常用来求行列式。
设有矩阵
a1n
M
ann
a11L
AM
n1L
a
用列主元消去法将其化为上三角形矩阵,对角线上元素为a11
(1),a22
(2),L,a(333),于是行
列式
detA
(1)m事a22)Lann)
其中m为所进行的行交换次数。
这是实际中求行列式值的可靠方法。
i+l=i
、
ajt,akja」,takj
jk,k1,...,n
blt,bkbi,tQ
1
I结束;
■
返回至程床图3伽选主元屆
d?
d
-bibL,bn丿
开始J
从主程序来
』-J3事本
立电F)精〔町痕直看㈣蛰賦巴xCl)=kOO()Gx
(2)=S.0000x[:
3)=3.0000
detA=-GG-0000
采用计算机运算在计算大数据时有明显的优点,另外也需要考虑到存储。
高斯消去法的使用条件是a:
:
)0,k1,2,L,n,而列选主元法可以保证这一
条件。
并且可以避免在消元过程确定乘数时所用除数是绝对值小的数,相对全选主元的运算量小,一般也可以满足精度要求。
此次上机不仅需要对原理了解透彻,而且要求的编程能力较强。
在定义和思路上没问题,只是在编程软件的使用上遇到些不稳定的问题,如头文件的使用。
在存储空间上得到了新的认识,另外发现了当代码多时流程框图的好处。
编程是一件很需要耐心的事,自己还有很大进步空间。
实验三例3-10
2.1杜里特尔分解法
设线性方程组Axb,对系数矩阵A进行除不交换两行位置得初等行变换相
当于用初等矩阵Mi左乘A,在对方程组第一次消元后,(A和(b分别化为(A和
b
(2),即
其中
消元过程是对k1~n1进行的,因此有
Mn1LM2M1A
(1)A(n)
Mn1LM2M1(b1)b(n)
将上三角形矩阵A(n)记为U,于是有
m43
为单位下三角形矩形。
这样高斯消去法的实质是将系数矩阵A分解为两个三角形矩阵L和U相乘,即ALU
在上述矩阵描述中遇到了下三角形矩阵运算。
主对角线以上元素全为零的方阵称为下三角形矩阵。
下三角形矩阵的乘积仍是下三角形矩阵。
若下三角形矩阵可逆,其逆矩阵仍是下三角形矩阵,而且下三角形矩阵的乘积和逆矩阵很容易求得。
把A分解成一个单位下三角阵和一个上三角阵U的乘积成为杜里特尔分解。
这种分解是惟一的。
2.2高斯-约当法
高斯消去法有消元和回代两个过程,当对消元过程稍加改变便可以使方程组化为对角形方程组
Dx的形式,其中矩阵D为对角形矩阵,即a1(11)
(2)
a22
O
a(nnn)nn的各元素(包括常数项),这个个过程成为归一化,这时方程组的系数阵转化为单位阵。
为减小误差,高斯-约当消去法还常用列选主元技术。
-CUtpLrt^.DCt-记亭萍
A的逆矩阵为・
2.0000-1.00000-0000
1.5000-0.50000-EOOO
2.5000-1.5000fkSOCX)
采用杜里特尔分解法求解方程组时,由于把对系数矩阵的计算和对右端项的计算分开,这使计算线性方程组系非常方便。
只需进行一次矩形三角分解,然后再解多个三个方程组,且多解一个方程组仅需要增加大约n2次乘除法运算。
采用高斯约当法仅需要进行消元归一,而不需要回代,为编程实现提供了便利。
步骤很重要,审题--确定算法--解题步骤--流程图--程序--代入简单值进行验证。
在编程时先在代码输入区打好框架,并且尽量在每一命令后注释,方便检查错误和日后复习。
定义和变量存储很灵活,如我把单位向量直接赋给了A矩阵
变量中,还有根据最终的目的直接简化计算。
另外赋值前,确定存储空间并且要定义初值为零。
实验四例4-6
已知f(X)的观测数据
X
3
4
f(X)
-5
-6
构造插值多项式。
首先构造基函数
lk(X)
Xi
,可以证明基函数满足下列条件:
对于给定n1个节点,
i0Xk
ik
l(X)
ki
1ik
n次拉格朗日插值多项式由下式给出:
n
L(x)lk(x)yk
k0
XX
i0XkX
由于Ik(X)是一个关于X的n次多项式,所以L(X)为关于X的不高于n次的代数多
项式。
当XXi时,L(Xi)yi,满足插值条件。
V
_'
Output4.tKt-i3事本
交申F)騒㈢帘式心1童看MS期(H)
疋二久hOOOOOO处的函数值为:
尸-&
3馬0000
由于所知的拉格朗日计算机算法只能实现计算某一特定值的近似函数值,而不知如何导出表达式,故例求x=2.5处的函数值以说明表达式以得出,只是在计算机程序中。
并且也能达到拉格朗日插值法使用的目的。
编程不够细心,程序没问题,却因为不知道是输入文件错了而检查了好长时间。
但同时也加深了对拉格朗日基函数性质的认识和理解。
实验五习题5-2
给出平面函数z(x,y)axbyc的数据
i
5
0.1
0.2
0.4
0.6
0.9
yi
0.3
0.5
0.7
0.8
z
0.58
0.63
0.73
0.83
0.92
按最小二乘原理确定a,b,c。
2.1最小二乘原理
设已知某物理过程yf(X)的一组观测数据
(Xi,f(Xi)),i1,2L,m
要求在某特定函数类(X)中寻找一个函数(X)作为yf(X)的近似函数,使得
二者在Xi上的误差或称残差
i(Xi)f(x),i1,2L,m
按某种度量标准为最小,这就是拟合问题。
要求残差
i按某种度量标准为最小,即要求由残差
i构成的残差向量
01,L,mT的某种范数II
II
为最小,这本来都是很自然的,
为最小。
例如,要求
max
I4,或III即
max(X)f(Xi)
可是计算不太方便。
所以通常要求:
23
(Xi)f(Xi)
i0
m
21
或者
2m
m(X)f(X)
这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法。
2.2多变量数据拟合
对于给定的一组数据(Xi,yj,i=1,2…,
m,寻求做n次多项式
akX
使性能指标
1z011n
J(a,a丄,a)
(y
Xk)2为最小。
k=0,
由于性能指标J可以被看做关于ak,
合多项式的构造问题可转化为多元函数的极值问题。
令
1,…
n的多元函数,故上述拟
Jak
从而有正则方程组
ao
对多变量(或称多元)线性模型
进行了m次观测
这个称为回归方程组,写成矩阵形式
y=
Aa
1*y
yy2*,A
*
ym
x11
x
12
1m
x21
22
2m
xn1
xn2
nm
a0
当mn时,要确定一组ai,i
0,1,L,n,使之精确地满足
个方程,这是超定方
程组的问题,只能在最小平方误差的基础上确定i。
定义残差向量[1,2,L,m]T,则
S=y-Aa
其中y[y1,y2,L,ym]T代替输出向量。
取性能指标
使之最小,以此确定出a。
由
JSTS
(yAa)T(yAa)aTAyyTAa+aTATAa
a。
J=§
TS
T
=yy
利用向量和矩阵的运算公式,有
AAa=Ay
此即为正则方程组,当ATA非奇异时,可求得
a=(AtA)1Aty
9]OLttpLfT5.CKt-爭本
文件旧舷目褂式◎S1(V)
所求参数为:
沪a2000b^O.3000c=0.5000
曲线拟合的最小二乘法是反映所给数据点的总的趋势,并不是严格的通过每个数据点,这样就避免了大量数据插值时需要高次多项式,同时又去掉了数据所含的测量误差。
整理思路越来越熟练了,所以执行各个步骤也相对简单了很多。
另外对原理也加深了认识。
#inc_ude=ssio.h=
整nc_ude=3afh.h=
#defineN30
void3ain()
壬二
floafXLNLG
F_LEAPUPN
fplHfopenunpun-Xr-Jr)
fp2Hfopen(=oufpun.S-JW-)
fscanf(fp二-%f"
cc=
fscanf(fp」r-cxs辽裆血
fprinff(fp2J
for(naiAW++)/、卡宜旃并
x=+」帀x=rx=rQ(2*c*x宇」=
fprinff(fp2JkHwdux(%d)H%.5f\rr」HX=2
if(rabs(x=+=-x=DAH0.0§
5)
break-
e-se
continue-
fprinff(fp2>
n云w%fHW.3f\n?
-GX=i2
fc-ose(fpD
fc_ose(fp2=
2>
输入文件:
"
input1.txt"
180.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<
n;
i++)for(j=0;
j<
j++)fscanf(fp1,"
%f"
a[i][j]);
i++)fscanf(fp1,"
b[i]);
//***数据输入*///*****************************************************************************
*************
高斯消去*/
消元*/
列选主元函数*/
//*****************************************//****************************for(k=0;
k<
n-1;
k++)//从第一次消元到第N-1次消元
d=a[k][k];
l=k;
for(i=k+1;
i++)//找出绝对值最大的a[i][k]和i行
if(fabs(a[i][k])>
fabs(d)){d=a[i][k];
i=l;
}if(i==n)//判断是否奇异,不奇异进行行交换
if(d==0)fprintf(fp2,"
奇异"
//如果所有行的“首列”都为0,为奇异elseif(l!
=k)//如果第k行的“首列”并不是最大
det=det*(-1);
for(j=k;
=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++)//第(k+1)次消元要得到N-(k+1)个乘数
a[i][k]=a[i][k]/a[k][k];
for(j=k+1;
j++)//
第(i+1)行各列向量对应该行与第(k+1)行各列
向量的减法
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;
i++)//输出方程组的解
fprintf(fp2,"
x(%d)=%.4f\n"
i+1,b[i]);
i<n;
i++)det=det*a[i][i];
detA=%.4f\n"
det);
//输出系数矩阵行列式的值
fclose(fp1);
fclose(fp2);
//************************************
2>输入文件:
“input2.txt”
12-33-183-1
15-156
3>输出文件:
“output2.txt”
x
(1)=1.0000x
(2)=2.0000x(3)=3.0000detA=-66.0000
3.例3-10
1>
源程序:
#defineN30voidmain()inti,j,r,k,n;
floata[N][N]={0},s;
input3.txt"
output3.txt"
i++)//输入单位阵
j=i+n;
a[i][j]=1;
*************LU
分解*/
for(i=1;
i++)//r=O时:
"
1"
区间不变化,"
2"
区间变化。
a[i][0]=a[i][0]/a[0][0];
i++)//第二行列至末行列进行变化时
for(j=i;
2*n;
j++)//第"
2(r+1)-1"
区间的变化行
{s=0.
for(k=0;
i;
k++)//对应行列元素的乘积进行求和
s+=a[i][k]*a[k][j];
a[i][j]=a[i][j]-s;
j++)//"
2(r+1)"
区间的变化列
s+=a[k][i]*a[j][k];
a[j][i]=(a[j][i]-s)/a[i][i];
}//现在已将AI--->
LUY
//*****************************************************************************
a[0][0]=1;
//第一列其余行已为零
*********高斯约当法解Ux=Y*/
i++)
for(j=0;
j++)a[i][j]=0;
j++)//首行归一化
a[0][j]=a[0][j]/a[0][0];
i++)//(n-1)次归一消元
j++)//第i+1行的各列进行归一化
a[i][j]=a[i][j]/a[i][i];
a[i][i]=1;
for(r=0;
r<
r++)//对第一行至第i-1行的i列进行置零
for(j=r+2;
j++)//r行的各列与第r+1行的对应列进行减法运算
a[r][j]=a[r][j]-a[i][j]*a[r][r+1];
//第r+1行归一后,乘数即为要置零处的值
a[r][r+1]=0;
//乘数用完之后置零即可
数据输出*/
//************************************fprintf(fp2,"
A的逆矩阵为\n\n"
i++)fprintf(fp2,"
|"
j++)fprintf(fp2,"
%
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 上机 实验 报告