计算方法实验.docx
《计算方法实验.docx》由会员分享，可在线阅读，更多相关《计算方法实验.docx（17页珍藏版）》请在冰豆网上搜索。
计算方法实验
计算方法实验报告
目录
计算方法实验报告-1-
实验一舍入误差与数值稳定性-2-
实验二方程求根-3-
实验三线性方程组数值解法-4-
实验四插值法-8-
实验五曲线拟合-9-
实验六数值积分与数值微分-11-
实验一舍入误差与数值稳定性
设
=
已知其精确值为
（
.
1）编制从大到小的顺序计算
的程序；
2）编制按从小到大的顺序计算
的程序；
3）按2种顺序分别计算
，
，
，并指出有效位数。
思路分析：
两种不同的累加方式。
解：
1）程序如下：
（预设N=100）
N=100;
S=0;
while（N-1）
S=S+1/（N^2-1）;
N=N-1;
end
fprintf（'%f',S）;
2）程序如下：
N=100;
S=0;
forj=2:
N
S=S+1/（j^2-1）;
end
fprintf（'%f',S）;
3）降序：
=0.749000
=0.749900
=0.749967
升序：
=0.749000
=0.749900
=0.749967
实验二方程求根
用牛顿法求下列方程的根：
思路分析：
首先应是牛顿迭代式，之后判断收敛即可。
解：
程序如下：
eps=5e-6;
delta=1e-6;
N=100;
k=0;
x0=1;
while
（1）
x1=x0-func1（x0）/func2（x0）;
k=k+1;
if（k>N|abs（x1）disp（'newtonmethodfailed'）
break;
end
ifabs（x1）<1
d=x1-x0;
else
d=（x1-x0）/x1;
end
x0=x1;
if（abs（d）break;
end
end
fprintf（'%f',x0）;
functiony=func2（x）
y=2*x-exp（x）;
end
functiony=func1（x）
y=x*x-exp（x）;
end
运行结果：
-0.703467
实验三线性方程组数值解法
1.编写用追赶法解三对角线性方程组的程序，并解下列方程：
Ax=b,其中，
=[-41
=[-27
1-41-15
1-41-15
…………
1-41-15
1-4]-15]
思路分析：
按照追赶法的递推公式求解。
解：
程序如下:
functionx=Trid（a,b,c,d）
%追赶法求解三对角的线性方程组Ax=d
%b为主对角线元素，a，c分别为次对角线元素，d为右端项
%A=[b1c1
%a2b2c2
%......
%a_（n-1）b_（n-1）c_（n-1）
%a_（n）b_（n）]
%b=[b1...b_（n）]
%a=[0a2...a_（n）]
%c=[c1...c_（n-1）]
n=length（b）;
u
（1）=b
（1）;
fori=2:
n
l（i）=a（i）/u（i-1）;
u（i）=b（i）-l（i）*c（i-1）;
end
y
（1）=d
（1）;
fori=2:
n
y（i）=d（i）-l（i）*y（i-1）;
end
x（n）=y（n）/u（n）;
fori=n-1:
-1:
1
x（i）=（y（i）-c（i）*x（i+1））/u（i）;
end
fori=1:
n
fprintf（'x[%d]=%f\n',i,x（i））;
end
运行结果：
x[1]=8.705758,x[2]=7.823032,x[3]=7.586371,x[4]=7.522453,x[5]=7.503440,x[6]=7.491306,x[7]=7.461785,x[8]=7.355835,x[9]=6.961556,x[10]=5.490389
2.分别用雅克比迭代法与高斯-赛德尔迭代法解下列方程组：
RI=V,其中
R=[31-13000-10000V=[-15
-1335-90-11000027
0-931-1000000-23
00-1079-30000-90
000-3057-70-50-20
0000747-300012
00000-304100-7
0000-50027-27
0000000-229]-10]
思路分析：
分别按照雅克比与高斯迭代法递推式编制程序。
解：
程序如下：
雅克比迭代法
functionjacobi（A,b,max,eps）
n=length（A）;x=zeros（n,1）;x1=zeros（n,1）;k=0;
while1
x1
（1）=（b
（1）-A（1,2:
n）*x（2:
n,1））/A（1,1）;
fori=2:
n-1
x1（i）=（b（i）-A（i,1:
i-1）*x（1:
i-1,1）-A（i,i+1:
n）*x（i+1:
n,1））/A（i,i）;
end
x1（n）=（b（n）-A（n,1:
n-1）*x（1:
n-1,1））/A（n,n）;
k=k+1;
ifsum（abs（x1-x））fprint（'number=%d\n',k）;
break;
end
ifk>=max
fprintf（'themethodisdisconvergent\n'）;
break;
end
x=x1;
end
ifkfori=1:
n
fprintf（'x[%d]=%f\n',i,x1（i））;
end
end
运行结果：
x[1]=-0.200550,x[2]=0.368393,x[3]=-0.731860,x[4]=-0.300318,x[5]=-0.446577,x[6]=0.399384,x[7]=0.121500,x[8]=0.151792,x[9]=-0.334359
高斯-赛德尔迭代法
functionGauseSeidel（A,b,max,eps）
n=length（A）;x=zeros（n,1）;x1=zeros（n,1）;k=0;
while1
x1
（1）=（b
（1）-A（1,2:
n）*x（2:
n,1））/A（1,1）;
fori=2:
n-1
x1（i）=（b（i）-A（i,1:
i-1）*x1（1:
i-1,1）-A（i,i+1:
n）*x（i+1:
n,1））/A（I,I）;
end
x1（n）=（b（n）-A（n,1:
n-1）*x1（1:
n-1,1））/A（n,n）;
k=k+1;
ifsum（abs（x1-x））fprintf（'number=%d\n',k）;
break;
end
ifk>=max
fprintf（'themethodisdiscomvergent\n'）;
break;
end
x=x1;
end
ifkfori=1:
n
fprint（'x[%d]=%f\n',i,x1（i））;
end
end
运行结果：
x[1]=-0.200550,x[2]=0.368393,x[3]=-0.731860,x[4]=-0.300318,x[5]=-0.446577,x[6]=0.399384,x[7]=0.121500,x[8]=0.151792,x[9]=-0.334359
实验四插值法
按下列数据做五次插值，并求x1=0.46,x2=0.55,x3=0.60时的函数近似值。
Xi
0.30
0.42
0.50
0.58
0.66
0.72
yi
1.04403
1.08462
1.11803
1.15603
1.19817
1.23223
思路分析：
用牛顿插值法多项式来解决此问题。
解：
程序如下：
functionnewtint（x,y,xhat）
n=length（y）;
c=y（:
）;
forj=2:
n
fori=n:
-1:
j
c（i）=（c（i）-c（i-1））/（x（i）-x（i-j+1））;
end
end
yhat=c（n）;
fori=n-1:
-1:
1
yhat=yhat*（xhat-x（i））+c（i）;
end
fprintf（'N（%f）=%f',xhat,yhat）;
运行结果：
x=0.46，y=1.100742;
x=0.55,y=1.141271;
x=0.60,y=1.166194;
实验五曲线拟合
试分别用抛物线y=a+bx+c
和指数曲线y=a
拟合下列数据：
X
1
1.5
2
2.5
3
3.5
4
4.5
y
33.4
79.5
122.65
159.05
189.15
214.15
238.65
252.50
x
5
5.5
6
6.5
7
7.5
8
y
267.55
280.5
296.65
301.40
310.4
318.15
325.15
比较2个拟合函数的优劣。
解：
1.抛物线拟合程序如下：
functionZXE（x,y,m）
S=zeros（1,2*m+1）;T=zeros（m+1,1）;
fork=1:
2*m+1
S（k）=sum（x.^（k-1））;
end
fork=1:
m+1
T（k）=sum（x.^（k-1）.*y）;
end
A=zeros（m+1,m+1）;a=zeros（m+1,1）;
fori=1:
m+1
forj=1:
m+1
A（i,j）=S（i+j-1）;
end
end
a=A\T;
fork=1:
m+1
fprintf（'a[%d]=%f\n',k,a（k））;
end
运行结果：
a[1]=-45.333297，a[2]=94.230200，a[3]=-6.131610；即：
y=-45.333297+94.230200x-6.131610
2.指数曲线拟合程序如下：
functionZXEI（x,y,m）
y=log（y）;
S=zeros（1,2*m-1）;T=zeros（m,1）;
fork=1:
2*m-1
S（k）=sum（x.^（k-1））;
end
fork=1:
m
T（k）=sum（x.^（k-1）.*y）;
end
A=zeros（m,m）;a=zeros（m,1）;
fori=1:
m
forj=1:
m
A（i,j）=S（i+j-1）;
end
end
a=A\T;
a
（1）=exp（a
（1））;
fork=1:
m
fprintf（'a[%d]=%f\n',k,a（k））;
end
运行结果：
a[1]=67.402,a[2]=0.238960，即：
y=67.402
实验六数值积分与数值微分
用复化梯形公式和复化辛普森公式计算积分
和
观察n为多少时所得近似值具有6位有效数字。
解：
MATLAB程序如下：
复化梯形公式
functiontrap（a,b）
n=1030;
fori=1:
3
x=a:
（b-a）/n:
b;
h=（b-a）/n;
s=（h/2）*（f（a）+2*sum（f（x（2:
1:
n-1）））+f（b））;
fprintf（'s（%d）=%f\n',n,s）;
n=n+1;
end
functiony=f（x）
y=（1+cos（x）.^2）.^0.5;
运行结果：
s（1030）=1.908574
s（1031）=1.908575
s（1032）=1.908577
n=1030时才有6位有效数字。
复化辛普森公式
functionsimpson（a,b）
n=2;
fori=1:
3
x=a:
（b-a）/（2*n）:
b;
m=2*n+1;
h=（b-a）/n;
s=（h/6）*（f（a）+2*sum（f（x（3:
2:
m-2）））+4*sum（f（x（2:
2:
m-1）））+f（b））;
fprintf（'s（%d）=%f\n',n,s）;
n=n*2;
end
functiony=f（x）
y=（1+cos（x）.^2）.^0.5;
运行结果：
s
（2）=1.910141
s（4）=1.910099
s（8）=1.910099
n=4时即有6位有效数字；
复化梯形公式
functiontrap（a,b）
n=1954;
fori=1:
3
x=a:
（b-a）/n:
b;
h=（b-a）/n;
s=（h/2）*（1+2*sum（f（x（2:
1:
n-1）））+f（b））;
fprintf（'s（%d）=%f\n',n,s）;
n=n+1;
end
functiony=f（x）
y=tan（x）/x;
运行结果：
s（1954）=0.848456
s（1955）=0.848456
s（1956）=0.848456
n=1954时才有6位有效数字。
复化辛普森公式
functionsimpson（a,b）
n=4;
fori=1:
3
x=a:
（b-a）/（2*n）:
b;
m=2*n+1;
h=（b-a）/n;
s=（h/6）*（1+2*sum（f（x（3:
2:
m-2）））+4*sum（f（x（2:
2:
m-1）））+f（b））;
fprintf（'s（%d）=%f\n',n,s）;
n=n*2;
end
functiony=f（x）
y=tan（x）./x;
运行结果：
s（4）=0.848972
s（8）=0.848967
s（16）=0.848967
n=8时即有6位有效数字；
升级会员 