河北工业大学数值分析实验报告.docx
《河北工业大学数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《河北工业大学数值分析实验报告.docx（23页珍藏版）》请在冰豆网上搜索。
河北工业大学数值分析实验报告
河北工业大学
数值分析
2014
实验报告
学院:
计算机科学与软件学院
班级:
姓名:
学号:
实验一
1.实验名称
用两种不同的顺序计算
，分析其误差的变化。
2.实验目的
（1）深刻理解舍入误差所引起的数值的不稳定性。
（2）理解利用计算机进行数值计算中舍入误差所引起的数值不稳定性。
（3）理解初始小的舍入误差可能造成误差累积，从而对计算结果造成巨大影响。
（4）理解对同一问题采用不同的算法，由于舍入误差的影响，也可能得到截然不同的结果。
3.算法描述
（1）用float精度正序相加得到10000项的和
（2）用float精度逆序相加得到10000项的和
（3）用double精度正序相加得到10000项的和
（4）用double精度逆序相加得到10000项的和
（5）根据结果分析并比较问题
4.源程序
#include
main（）
{
floati,s1,s2,n;
doubles3,s4;
s1=0;
s2=0;
s3=0;
s4=0;
for（i=1;i<=10000;i++）
{
n=i*i;n=1/n;s1=s1+n;
}
printf（"用float精度正序相加结果为:
%f\n",s1）;
n=0;
for（i=10000;i>=1;i--）
{
n=i*i;n=1/n;s2=s2+n;
}
printf（"用float精度倒序相加结果为:
%f\n",s2）;
n=0;
for（i=1;i<=10000;i++）
{
n=i*i;n=1/n;s3=s3+n;
}
printf（"用double精度正序相加结果为:
%f\n",s3）;
n=0;
for（i=10000;i>=1;i--）
{
n=i*i;n=1/n;s4=s4+n;
}
printf（"用double精度倒序相加结果为:
%f\n",s4）;
}
5.运行结果
6.编程体会
舍入误差会引起数值计算的不稳定性。
舍入误差对于计算结果精确度影响小的算法具有较好的稳定性；然而对于计算结果精确度影响大的算法，则会产生误差的积累。
就本实验来说，分别使用float型还有double型算法，由于小数点后保留的位数比较多，并没有出现太大的误差，也避免了大数吃小数的现象。
我在计算机上运用int类型也计算了一下，发现正序倒序结果都是1，这就说明由于舍入误差的累积，带来了计算的不稳定性。
因此运用计算机解决问题的时候要采用合适的算法，才能得到更期望的结果。
实验二
1.实验名称
1拉格朗日插值
按下列数据
x
-3.0
-1.0
1.0
2.0
3.0
y
1.0
1.5
2.0
2.0
1.0
作二次插值，并求x
=-2，x
=0，x
=2.75时的函数近似值
2牛顿插值
按下列数据
x
0.30
0.42
0.50
0.58
0.66
0.72
y
1.04403
1.08462
1.11803
1.15603
1.19817
1.23223
作五次插值，并求x
=0.46，x
=0.55，x
=0.60时的函数近似值。
2.实验目的
（1）了解拉格朗日插值和牛顿插值的计算方法。
（2）比较两种算法并分析优点。
3.算法描述
1．拉格朗日插值
1）从五个点中选取三个点作为已知点
2）由选出的三个点根据拉格朗日公式求出对应的系数L0,L1,L2；
3）将所得系数代入公式，计算结果。
2．牛顿插值
1）根据所给的六组数据代入，分别求其对应的差商；
2）将对应求出的系数a0,a1,a2,a3,a4,a5代入公式；
3）将所求数据代入，计算结果。
4.源程序
1.拉格朗日插值
#include
main（）
{
charb;
floatx[50],y[50],t[50],q[50],s,p;
inti,j,k,n,m;
printf（"您选择的是拉格朗日插值：
\n"）;
printf（"题目：
按下列数据\n"）;
printf（"x：
-3.0-1.01.02.03.0\n"）;
printf（"y：
1.01.52.02.01.0\n"）;
printf（"作二次插值，并求x=-2，x=0，x=2.75时的函数近似值。
\n\n"）;
printf（"请输入插值次数：
"）;
scanf（"%d",&n）;
for（i=0;i<=n;i++）
{printf（"请输入第%d个节点：
",i+1）;
scanf（"%f,%f",&x[i],&y[i]）;}
te1:
printf（"请输入要求的X的个数：
"）;
scanf（"%d",&k）;
for（m=0;m{printf（"请输入x%d：
",m+1）;
scanf（"%f",&t[m]）;}
for（m=0;m{s=0;
for（i=0;i<=n;i++）
{p=1;
for（j=0;j<=n;j++）
{if（i!
=j）p=p*（（t[m]-x[j]）/（x[i]-x[j]））;}
s=s+p*y[i];}
q[m]=s;}
for（m=0;mprintf（"Y%d的值为：
%f\n",m+1,q[m]）;
}
2．牛顿插值
#include
main（）
{
floatx[50],y[50][50],q[50],t[50],a,s,p;
inti,j,n,m,k;
printf（"您选择的是牛顿插值：
\n"）;
printf（"题目：
按下列数据\n"）;
printf（"x：
0.300.420.500.580.660.72\n"）;
printf（"y：
1.044031.084621.118031.156031.198171.23223\n"）;
printf（"作五次插值，并求x=0.46，x=0.55，x=0.60时的函数近似值。
\n\n"）;
printf（"请输入插值的结点个数：
"）;
scanf（"%d",&k）;
for（i=0;i{printf（"请输入第%d个节点：
",i+1）;
scanf（"%f,%f",&x[i],&y[0][i]）;}
for（n=1;n{for（i=n;iy[n][i]=（y[n-1][i]-y[n-1][i-1]）/（x[i]-x[i-n]）;}
printf（"请输入要求的X的个数：
"）;
scanf（"%d",&j）;
for（i=0;i{printf（"请输入x%d：
",i+1）;
scanf（"%f",&t[i]）;}
for（m=0;m{s=0;
for（n=1;n{p=1;
for（i=0;i{a=t[m]-x[i];
p=p*a;}
p=p*y[n][n];
s=s+p;}
q[m]=s+y[0][0];}
for（i=0;iprintf（"Y%d的值为：
%f\n",i+1,q[i]）;
}
5.运行结果
6.编程体会
两种方法都是通过先构造插值函数然后再计算求解。
通过实验发现牛顿插值比拉格朗日插值的精度要高。
拉格朗日插值是通过对n+1个n次基本多项式线性组合得到的，而牛顿插值法是通过求各阶差商，递推得到的。
用拉格朗日插值多项式计算函数近似值时，如需增加插值节点，那么原来算出的数据均不能利用，必须重新计算。
牛顿插值法比拉格朗日插值法节省计算量，且便于程序设计，计算增加节点时，计算只需要增加一项，并且牛顿插值更容易计算高次插值。
实验三
1.实验名称
分别用复化梯形公式和复化辛卜生公式计算f（x）=sin（x）/x的积分，并与准确值比较判断精度。
2.实验目的
（1）通过实际计算掌握复化梯形和复化辛卜生求积公式。
（2）比较复化梯形和复化辛卜生求积公式的计算精确。
3.算法描述
（1）输入积分的上下限
（2）输入划分区间的个数
（3）让计算机计算步长h=（b-a）/n
（4）带入复化公式计算结果
4.源程序
#include
#include
#include
#defineI0.9460831
main（）{
intn1,n2,k;
floata,b,h1,h2,d1,d2;
doublet1=0,t2=0,p=0,T=0,S=0;
printf（"请输入积分下限a="）;
scanf（"%f",&a）;
printf（"请输入积分上限b="）;
scanf（"%f",&b）;
if（a==b）
{printf（"积分上下限相等，定积分的值为0\n"）;
exit（0）;}
printf（"请输入梯形公式划分区间的个数n1="）;
scanf（"%d",&n1）;
printf（"请输入辛卜生公式划分区间的个数n2="）;
scanf（"%d",&n2）;
h1=（b-a）/n1;
h2=（b-a）/n2;
printf（"梯形公式步长h1=%f\n",h1）;
printf（"辛卜生公式步长h2=%f\n",h2）;
for（k=1;k{t1+=sin（a+k*h1）/（a+k*h1）;
}
for（k=1;k{
t2+=sin（a+k*h2）/（a+k*h2）;
}
for（k=0;k{
p+=sin（a+（2*k+1）*h2/2）/（a+（2*k+1）*h2/2）;
}
if（a==0）{
T=（1+2*t1+sin（b）/b）/2*h1;
S=（1+2*t2+4*p+sin（b）/b）/6*h2;}
elseif（b==0）{
T=（sin（a）/a+2*t1+1）/2*h1;
S=（sin（a）/a+2*t2+4*p+1）/6*h2;}
else{
T=（sin（a）/a+2*t1+sin（b）/b）/2*h1;
S=（sin（a）/a+2*t2+4*p+sin（b）/b）/6*h2;
}
printf（"复化梯形公式T=%f\n",T）;
printf（"复化辛普森公式S=%f\n",S）;
printf（"准确值I=%f\n",I）;
d1=fabs（I-T）;
d2=fabs（I-S）;
printf（"复化梯形公式的误差d=%f\n",d1）;
printf（"复化辛卜生公式的误差d=%f\n",d2）;
if（d1elseif（d1==d2）printf（"复化辛卜生公式和复化梯形公式的精度相等\n"）;
elseprintf（"复化辛卜生公式的精度高于复化梯形公式的精度\n"）;
}
5.运行结果
6.编程体会
高次插值会产生龙格现象，故采用分段低次插值，然后利用积分的区间可加性得积分，即所谓的分段低次合成的复化求积公式。
复化求积法能改善求积精度。
通过编程，我体会到复化辛卜生公式的代数精度要高于复化梯形公式，复化辛卜生公式的计算结果与准确值更接近。
因此，实际应用中应采用复化辛卜生公式计算积分，可以得到与准确值更接近的近似值。
编程中要注意的问题：
（1）这里要注意讨论a=0或者b=0或者a=b的情况，因为公式中有sinx/x，不讨论a、b为0的情况可能会出现越界。
（2）复化公式中有求和符号，写for循环就可以解决，要特别注意循环的次数。
实验四
1.实验名称
用改进欧拉方法解初值问题y’=x+y;y（0）=1。
02.实验目的
（1）掌握改进欧拉公式，并熟练应用。
（2）利用改进欧拉公式求解常微分方程的初值问题，并与准确值比较结果。
3.算法描述
（1）yp=yi+h*f（xi,yi）
（2）yc=yi+h*f（xi+1,yp）
（3）yi+1=（yp+yc）/2
4.源程序
#include
#include
voidmain（）
{
inti;
floatyp,yc,y=1,h,x=0;
printf（"请输入步长h="）;
scanf（"%f",&h）;
for（i=0;i<10;i++）
{x=h*i;
yp=y+h*（x+y）;
yc=y+h*（x+h+yp）;
y=（yp+yc）/2;
printf（"实验值y（%0.1f）=%f\n",h*（i+1）,y）;
printf（"
精确值y（%0.1f）=%f\n",h*（i+1）,-h*（i+1）-1+2*exp（h*（i+1）））;
}
}
5.运行结果
6.编程体会
改进欧拉法首先运用欧拉法求出预测值然后再利用梯形公式进行校正，每迭代一次就校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式，有效的提高了计算精度。
设计计算方案主要包括六个环节：
预测、改进、计算、校正、改进、计算。
编程的时候要注意用平均化形式来表示改进欧拉公式，更易于机器操作。
实验五
1.实验名称
分别用下列方法求f（x）=x3-3x-1=0在x0=2附近的根。
根的准确值为x*=1.87938524…，要求准确到四位有效数字，并对比各种算法的计算量。
（1）二分法
（2）简单迭代法
（3）牛顿迭代法
2.实验目的
（1）掌握二分法、简单迭代法和牛顿迭代法。
（2）比较三种方法的特点、计算精度及计算量。
3.算法描述
1.二分法
（1）准备计算f（x）在有根区间[a,b]端点处的值f（a）,f（b）.
（2）二分计算f（a）在区间中点（a+b）/2处的值f[（a+b）/2]。
（3）判断若f[（a+b）/2]=0,则（a+b）/2为根，计算结果结束，否则检验；若f[（a+b）/2]<0,则以（a+b）/2代替b,否则替代a.
（4）反复执行步骤二和步骤三，直到区间的长度小于允许误差，此时的中点即为所求。
2.简单迭代法
（1）将方程改写成等价的形式，即变为迭代函数。
（2）代入一个初始近似值X0,反复进行多次迭代求出其根。
3.牛顿迭代法
（1）准备选定初始近似值X0,并计算其函数值，及其对应的导数；
（2）迭代按公式X1=X0-f0/f0ˊ迭代一次，得到新的近似值X1,并求其函数值，及其相对应的导数值；
（3）控制如果X1满足所需的误差，则终止迭代，即X1为所求的根，否则转步骤四；
（4）修改如果迭代次数预先达到指定的次数N,则此方法失败，否则转步骤二继续迭代。
4.源程序
#include
#include
voidfun1（）
{
doublex0,x1=1.0,x2=2.0,fx0,fx1,fx2;
intcount=0;
fx1=x1*x1*x1-3*x1-1;
fx2=x2*x2*x2-3*x2-1;
do{
x0=（x1+x2）/2;
printf（"第%d个二分点为%f\n",count+1,x0）;
fx0=x0*x0*x0-3*x0-1;
if（fx0==0）{
printf（"失败"）;break;}
if（fx0*fx1<0）{
x2=x0;fx2=fx0;}
else{x1=x0;fx1=fx0;}
count++;
}
while（fabs（x1-x2）>5e-4）;
printf（"共迭代了%d次\n",count）;
}
voidfun2（）
{
doublex0,x1=2.0;
intcount=0;
do
{
x0=x1;
x1=exp（log（3*x0+1）/3）;
count++;
}
while（fabs（x1-x0）>5e-4）;
printf（"简单迭代法x=%f\n共迭代了%d次\n",x1,count）;
}
voidfun3（）
{
doublex1,x2=2.0,f1,f2;
intcount=0;
do
{x1=x2;
f1=x1*x1*x1-3*x1-1;
f2=3*x1*x1-3;
x2=x1-f1/f2;
count++;
}
while（fabs（x1-x2）>5e-4）;
printf（"牛顿迭代法x=%f\n共迭代了%d次\n",x2,count）;
}
voidmain（）
{
fun1（）;
fun2（）;
fun3（）;
}
5.运行结果
6.编程体会
（1）二分法是逐步搜索方法的一种改进，在有根区间[a，b]上无限地进行二分过程，这些区间最终必将收缩于一点，该点就是所求的根。
（2）简单迭代法是一种逐次逼近法，其基本思想是将隐式方程归结为一组显式的计算公式，迭代过程实质上就是一个逐步显式化的过程，并且不好的迭代公式可能造成最后得不到计算结果。
（3）牛顿迭代法是将非线性方程逐步化为线性方程来求解，可以加速迭代速度。
但是牛顿迭代法依赖于初始值x0的选取，如果x0距离准确值接近，则迭代次数少，如果距离准确值远，则迭代次数很多。
并且，如果方程有两个跟，不同的初始值可能最终会收敛到别的根。
（4）总的来说，二分法算法简单，易于操作，但是精度较低，计算量大;简单迭代法精度较二分法有所提高但在迭代时需要选择合适的迭代公式，以得到较快的收敛速度；牛顿迭代法精度最高，计算量小，但是初始值的选取很重要。
实验六
1.实验名称
分别用高斯列主元消去法和直接三角分解法（LU分解）求方程组的解
系数矩阵：
10787常向量：
10
75658
861096
759107
精确解为：
（-60，102，-27，16）
2.实验目的
（1）掌握高斯列主元消去法和直接三角分解法（LU分解）求方程组解的方法
（2）比较高斯列主元消去法和直接三角分解法的特点
3.算法描述
1．高斯列主元消去法（以4阶为例）
第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式
第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为
第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为
按x4→x3→x2→x1的顺序回代求解出方程组的解。
2.直接三角分解法
将高斯消去法改写成紧凑的形式，可以直接从矩阵A的元素得到计算L,U元素的递推公式，而不需要任何中间步骤，这就是所谓的直接三角分解法。
一旦实现了矩阵A的LU分解，那么求解AX=b的问题就是等价于求解两个三角形方程组。
（1）Ly=b,求y
（2）Ux=y,求x
4.源程序
1.高斯消元法
#include
voidmain（）
{
floata[4][4]={{10,7,8,7},{7,5,6,5},{8,6,10,9},{7,5,9,10}},
y[4],c[4][4],x[4],d[4],m,b;
inti,n,j,f;
printf（"请输入右端顶:
\n"）;
for（i=0;i<=3;i++）
scanf（"%f",&y[i]）;
for（n=0;n<=2;n++）
{
m=a[n][n];
f=n;
for（i=（n+1）;i<=3;i++）
{
if（m{
m=a[i][n];
f=i;
}
}
if（f!
=n）
{
for（j=0;j<=3;j++）
{
c[n][j]=a[n][j];
a[n][j]=a[f][j];
a[f][j]=c[n][j];
}
d[n]=y[n];
y[n]=y[f];
y[f]=d[n];
}
for（i=（n+1）;i<=3;i++）
{
b=-a[i][j]/a[n][n];
for（j=0;j<=3;j++）
a[i][j]=a[n][j]*b+a[i][j];
y[i]=y[n]*b+y[i];
}
}
x[3]=y[3]/a[3][3];
x[2]=（y[2]-a[2][3]*x[3]）/a[2][2];
x[1]=（y[1]-a[1][3]*x[3]-a[1][2]*x[2]）/a[1][1];
x[0]=（y[0]-a[0][3]*x[3]-a[0][2]*x[2]-a[0][1]*x[1]）/a[0][0];
printf（"x1的值为%f\nx2的值为%f\nx3的值为%f\nx4%f\n",x[0],x[1],x[2],x[3]）;
｝
2.直接三角分解法
#defineN5
#include
voidmain（）
{
inti,j,r,k;
floata[N][N],l[N][N],u[N][N],b[N],x[N],y[N];
floatde=0;
printf（"pleaseimputthe%d*%dmatrixA:
\n",N-1,N-1）;
for（i=1;ifor（j=1;jscanf（"%f",&a[i][j]）;
printf（"pleaseimputthel*%dmatrixb:
\n",N-1）;
for（j=1;jscanf（"%f",&b[j]）;
for（j=1;j{
u[1][j]=a[1][j];
l[j][1]=a[j][1]/u[1][1];
}
for（r=2;r{
for（j=r;j{
for（k=1;k{
de+=l[r][k]*u[k][j];
}
u[r][j]=a[r][j]-de;
de=0;
}
for（i=r+1;r{
for（k=1;k
升级会员 