计算方法与实习报告.docx
《计算方法与实习报告.docx》由会员分享，可在线阅读，更多相关《计算方法与实习报告.docx（21页珍藏版）》请在冰豆网上搜索。
计算方法与实习报告
实习题1
1.用2种不同的顺序计算
，分析其误差的变化。
思路：
采用迭代法用C++来编程实现，一种用正序即n从1到10000取值，另一种用逆序
（一）正序：
采用递增的算法，将n从1到10000进行
的累加求和。
1）程序代码为：
#include
#include
usingnamespacestd;
voidmain（）
{
floatsum=0;
intn=1;
while（n<=10000）{
sum=sum+（float）1/（n*n）;
n++;
}
cout<<"Thesumthatfrommintomaxis:
"<}
（2）运行输出的结果为：
Thesumthatfrommaxtominis:
=1.64473
（二）逆序：
采用递减的算法，将n从10000到1进行
的累加求和。
（1）程序代码为：
#include
#include
usingnamespacestd;
voidmain（）
{
floatsum=0;
intn=10000;
while（n>0）{
sum=sum+（double）1/（n*n）;
n--;
}
cout<<"Thesumthatfrommaxtominis:
"<}
（2）运行输出的结果为：
Thesumthatfrommintomaxis:
1.64483
结果分析：
由运行的结果可以看出正序的结果为1.64473，逆序的结果为1.64483，而题目所给出的精确值为1.644834。
可以发现从小到大求和比从大到小求和更接近实际值，因为在从大到小的求解过程中，会出现大数“吃掉”小数的现象，从而造成较大误差。
实习题2
1：
用牛顿法求下列方程的根：
思路：
给定初始值x0，
为根的容许误差，
为|f（x）|的容许误差，N为迭代次数的容许值。
在满足根的容许误差，函数值的容许误差及迭代次数条件下，计算x1=x0－f（x0）/f'（x0）。
以此进行下次循环。
（1）程序代码为：
#include
#include
usingnamespacestd;
#defineN100
#defineeps1e-6
#defineeta1e-8
floatNewton（float（*f）（float）,float（*f1）（float）,floatx0）
{
floatx1,d;
intk=0;
do
{
x1=x0-（*f）（x0）/（*f1）（x0）;
if（k++>N||fabs（（*f1）（x1））{
printf（"\nNewton迭代发散"）;
break;
}
d=fabs（x1）<1?
x1-x0:
（x1-x0）/x1;
x0=x1;
printf（"x（%d）=%f\t",k,x0）;
}
while（fabs（d）>eps&&fabs（（*f）（x1））>eta）;
returnx1;
}
floatf（floatx）
{
returnfloat（x*x-exp（x））;
}
floatf1（floatx）
{
returnfloat（2*x-exp（x））;
}
voidmain（）{
inti=0;
floatx=0;
for（i=0;i<100;i++）{
x=x-A（x）/B（x）;
cout<<"x="<i++;
}
}
（2）运行过程：
x=-1
x=-0.733044
x=-0.703808
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
x=-0.703467
（3）运行结果为：
x=-0.703467
结果分析：
得出的结果为x=-0.703467，由运行过程可知在此值时结果已经稳定，且经检验结果精确到了小数点后第六位。
实习题3
1：
用列主元高斯消去法解方程组：
（1）x1+x2+3x4=4
2x1+x2–x3+x4=1
3x1–x2–x3+3x4=–3
–x1+2x2+3x3–x4=4
思路：
将方程用增广矩阵[A|b]=（aij）n
（n+1）表示,分别进行消元和回代过程
消元过程：
对k=1,2,.....,n-1
（1）选主元
（2）如果
=0，则矩阵A奇异，程序结束；否则执行（3）
（3）如果
，则交换第k行与第
行对应元素位置
（4）消元，对i=k+1,........,n计算
，对j=k+1,..........,n+1计算
回代过程
（1）若
，则矩阵A奇异，程序结束；否则执行
（2）
（2）
；对i=n-1,......2,1计算
（1）程序代码为：
#include
#include
#include
usingnamespacestd;
voidmain（）
{
inti;
floatx[4];
floatc[4][5]={1,1,0,3,4,
2,1,-1,1,1,
3,-1,-1,3,-3,
-1,2,3,-1,4};
voidDirect（float*,int,float[]）;
Direct（c[0],4,x）;
for（i=0;i<=3;i++）printf（"x[%d]=%f\n",i,x[i]）;
}
voidDirect（float*c,intn,floatx[]）
{
inti,j,k,t;
floatp;
for（i=0;i<=n-2;i++）
{
k=i;
for（j=i+1;j<=n-1;j++）
if（fabs（*（c+j*（n+1）+i））>（fabs（*（c+k*（n+1）+i））））k=j;
if（k!
=i）
for（j=i;j<=n;j++）
{
p=*（c+i*（n+1）+j）;
*（c+i*（n+1）+j）=*（c+k*（n+1）+j）;
*（c+k*（n+1）+j）=p;
}
for（j=i+1;j<=n-1;j++）
{
p=（*（c+j*（n+1）+i））/（*（c+i*（n+1）+i））;
for（t=i;t<=n;t++）*（c+j*（n+1）+t）-=p*（*（c+i*（n+1）+t））;
}
}
for（i=n-1;i>=0;i--）
{
for（j=n-1;j>=i+1;j--）
（*（c+i*（n+1）+n）-=x[j]*（*（c+i*（n+1）+j）））;
x[i]=*（c+i*（n+1）+n）/（*（c+i*（n+1）+i））;
}
}
（2）运行结果为：
x[0]=-1.333333
x[1]=2.333333
x[2]=-0.333333
x[3]=1.000000
4：
编写用追赶法解三对角线性方程组的程序，并解下列方程组。
（2）Ax=b其中
思路：
应用高斯消去法求解时，消元过程可以简化，得到同解的三角方程组，进行回代
消去过程为：
回代过程为
（1）程序代码为：
#include
voidmain（）
{
floatx[10];
inti;
floata[10][11]={-4,1,0,0,0,0,0,0,0,0,-27,
1,-4,1,0,0,0,0,0,0,0,-15,
0,1,-4,1,0,0,0,0,0,0,-15,
0,0,1,-4,1,0,0,0,0,0,-15,
0,0,0,1,-4,1,0,0,0,0,-15,
0,0,0,0,1,-4,1,0,0,0,-15,
0,0,0,0,0,1,-4,1,0,0,-15,
0,0,0,0,0,0,1,-4,1,0,-15,
0,0,0,0,0,0,0,1,-4,1,-15,
0,0,0,0,0,0,0,0,1,-4,-15};
voidDirect（float*,int,float[]）;
Direct（a[0],10,x）;
for（i=0;i<=9;i++）printf（"x[%d]=%f\n",i,x[i]）;
}
voidDirect（float*u,intn,floatx[]）
{
inti,r,k;
for（r=0;r<=n-1;r++）
{
for（i=r;i<=n;i++）
for（k=0;k<=r-1;k++）
*（u+r*（n+1）+i）-=*（u+r*（n+1）+k）*（*（u+k*（n+1）+i））;
for（i=r+1;i<=n-1;i++）
{
for（k=0;k<=r-1;k++）
*（u+i*（n+1）+r）-=*（u+i*（n+1）+k）*（*（u+k*（n+1）+r））;
*（u+i*（n+1）+r）/=*（u+r*（n+1）+r）;
}
}
for（i=n-1;i>=0;i--）
{
for（r=n-1;r>=i+1;r--）
*（u+i*（n+1）+n）-=*（u+i*（n+1）+r）*x[r];
x[i]=*（u+i*（n+1）+n）/（*（u+i*（n+1）+i））;
}
}
（2）运行结果为：
x[0]=8.705758
x[1]=7.823033
x[2]=7.586371
x[3]=7.522452
x[4]=7.503439
x[5]=7.491305
x[6]=7.461785
x[7]=7.355835
x[8]=6.961556
x[9]=5.490389
实习题4
2：
按下列数据：
xi0.300.420.500.580.660.72
yi1.044031.084621.118031.156031.198171.23223
作五次插值，并求x1=0.46，x2=0.55，x3=0.60时的函数近似值。
思路：
（1）输入
令
（2）对i=0,1,2,....,n计算
（1）程序代码为：
#include
#include
usingnamespacestd;
#defineN5
voidDifference（floatx[],floaty[],intn）
{
float*f=newfloat[n+1];
intk,i;
for（k=1;k<=n;k++）
{
f[0]=y[k];
for（i=0;if[i+1]=（f[i]-y[i]）/（x[k]-x[i]）;
y[k]=f[k];
}
deletef;
return;
}
voidmain（）
{
inti;
floatb,varx;
floatx[N+1]={0.30,0.42,0.50,0.58,0.66,0.72};
floaty[N+1]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};
Difference（x,y,N）;
b=y[N];
for（i=N-1;i>=0;i--）
b=b*（varx-x[i]）+y[i];
printf（"Nn（%f）=%f",varx,b）;
}
在main函数中：
分别另varx为0.46，0.55，0.60可得：
（2）运行结果为：
=0.46时，
=1.100724；
=0.55时，
=1.141271；
=0.60时，
=1.166194；
实习题5
1：
试分别用抛物线
和指数曲线
拟合下列数据：
xi11.522.533.544.5
yi33.479.50122.65159.05189.15214.15238.65252.50
xi55.566.577.58
yi267.55280.50296.65301.40310.40318.15325.15
比较2个拟合函数的优劣。
（1）用抛物线
拟合：
a.程序代码为：
matlab代码为：
>>A=[11.522.533.544.555.566.577.58];
B=[33.479.5122.65159.05189.15214.15238.65252.5267.55280.5296.65301.4310.4318.15325.15];
S=zeros（1,5）;T=zeros（3,1）;
>>fork=1:
5
S（k）=sum（A.^（k-1））;
end
>>fork=1:
3
T（k）=sum（A.^（k-1）.*B）;
end
>>C=zeros（3,3）;
>>a=zeros（3,1）;
>>fori=1:
3
forj=1:
3
C（i,j）=S（i+j-1）;
end
end
>>a=C\T;
>>fork=1:
3
disp（a（k））;
end
b.运行结果为：
a[0]=-45.3333
a[1]=94.2302
a[2]=-6.1316
c.拟合曲线：
（2）用指数曲线
拟合
将等式两边取自然对数得：
lny=lna+bx
记y=lnya=lnay=a+bx
得到新数据：
x11.522.533.544.555.566.577.58
y3.514.384.815.075.245.375.475.535.595.645.695.715.745.765.78
a.程序代码为：
#include
#include
#include
voidmain（）
{
voidColPivot（float*,int,float[]）;
inti;
floatx[2];
floatc[2][3]={15,67.5,79.29,
67.5,373.5,373.5};
ColPivot（c[0],2,x）;
for（i=0;i<=1;i++）printf（"x[%d]=%f\n",i,x[i]）;
}
voidColPivot（float*c,intn,floatx[]）
{
inti,j,k,t;
floatp;
for（i=0;i<=n-2;i++）
{
k=i;
for（j=i+1;j<=n-1;j++）
if（fabs（*（c+j*（n+1）+i））>（fabs（*（c+k*（n+1）+i））））k=j;
if（k!
=i）
for（j=i;j<=n;j++）
{
p=*（c+i*（n+1）+j）;
*（c+i*（n+1）+j）=*（c+k*（n+1）+j）;
*（c+k*（n+1）+j）=p;
}
for（j=i+1;j<=n-1;j++）
{
p=（*（c+j*（n+1）+i））/（*（c+i*（n+1）+i））;
for（t=i;t<=n;t++）*（c+j*（n+1）+t）-=p*（*（c+i*（n+1）+t））;
}
}
for（i=n-1;i>=0;i--）
{
for（j=n-1;j>=i+1;j--）
（*（c+i*（n+1）+n）-=x[j]*（*（c+i*（n+1）+j）））;
x[i]=*（c+i*（n+1）+n）/（*（c+i*（n+1）+i））;
}
}
b.运行结果为：
x[0]=4.208903
x[1]=0.239355
c.拟合曲线为：
分析：
带入数据进行检验，我们可以发现抛物线
与指数曲线
拟合相比，抛物线更加精确，而指数曲线具有较大的误差，因为在指数曲线拟合时，将指数函数按照一次函数来进行处理的，而拟合时，拟合函数的复杂度越高越准确，所以抛物线的拟合方法更准确。
实习题6
用复化梯形公式和复化Simpson公式计算积分
和
观察n为多少时所得近似值具有6位有效数字。
（1）用复化梯形公式计算积分
a.程序代码为
#include
#include
intn;
floatAutoTrap（float（*f）（float）,floata,floatb）
{
inti;
floatx,s,h=b-a;
floatt1,t2=h/2*（（*f）（a）+（*f）（b））;
n=1;
do
{
s=0;
t1=t2;
for（i=0;i<=n-1;i++）
{
x=a+i*h+h/2;
s+=（*f）（x）;
}
t2=（t1+s*h）/2;
n*=2;
h/=2;
}
while（fabs（t2-t1）>0.5e-6）;
returnt2;
}
floatf（floatx）
{
returnsqrt（1+cos（x）*cos（x））;
}
voidmain（）
{
floats;
doublepi=3.141592653;
s=AutoTrap（f,0,pi/2）;
printf（"T（%d）=%f\n",n,s）;
}
b.运行结果为：
T[8]=1.910099
结果分析：
运行结果表示n为8时所得近似值具有6位有效数字。
（2）复化Simpson公式计算积分
a.程序代码为
#include
#include
floatSimpson（float（*f）（float）,floata,floatb,intn）
{
intk;
floats,s1,s2=0;
floath=（b-a）/n;
s1=（*f）（a+h/2）;
for（k=1;k<=n-1;k++）
{
s1+=（*f）（a+k*h+h/2）;
for（k=1;k<=n-1;k++）
{
s1+=（*f）（a+k*h+h/2）;
s2+=（*f）（a+k*h）;
}
s=h/6*（（*f）（a）+4*s1+2*s2+（*f）（b））;
returns;
}
}
floatf（floatx）
{
if（x==0）return1;
elsereturnsin（x）/x;
}
voidmain（）
{
doublepi=3.141592653;
inti,n=2;
floats;
for（i=0;i<=2;i++）
{
s=Simpson（f,0,pi/4,n）;
printf（"s（%d）=%f\n",n,s）;
n*=2;
}
}
b.运行结果为：
s
（2）=1.005897
s（4）=0.887991
s（8）=0.824189
计算方法上机实验心得
通过用计算方法这门课的上机实验学习我对于c++和matlab程序的拥有一定掌握，并能完成一些在课本上的题目，如牛顿迭代法，列主元高斯消去法，雅可比迭代法，高斯赛德尔迭代法，插值法，复化梯形公式，复化Simpson公式等等。
在求解题目的过程中发现自己的编程能力逐步得到了提高，并且能够发现并解决一些简单的程序问题，同时加深了我们对课本上的知识的理解和认识，可以说是使得我们能够学以致用，加强了理论与实际的联系，我们从中获益良多。
升级会员 