数值分析上机题概况.docx
《数值分析上机题概况.docx》由会员分享，可在线阅读，更多相关《数值分析上机题概况.docx（26页珍藏版）》请在冰豆网上搜索。
数值分析上机题概况
数值分析上机题
习题1
17．（上机题）舍入误差与有效数
设
，其精确值为
。
（1）编制按从大到小的顺序
，计算
的通用程序。
（2）编制按从小到大的顺序
，计算
的通用程序。
（3）按两种顺序分别计算
，
，
，并指出有效位数。
（编制程序时用单精度）
（4）通过本上机题你明白了什么？
按从大到小的顺序计算
的通用程序为：
#include
floatsum（floatN）
{
floatj,s,sum=0;
for（j=2;j<=N;j++）
{
s=1/（j*j-1）;
sum+=s;
}
returnsum;
}
按从小到大的顺序计算
的通用程序为：
#include
floatsum（floatN）
{
floatj,s,sum=0;
for（j=N;j>=2;j--）
{
s=1/（j*j-1）;
sum+=s;
}
returnsum;
}
从大到小的顺序的值
从小到大的顺序的值
精确值
有效位数
从大到小
从小到大
0.740049
0.74005
0.740049
6
5
0.749852
0.7499
0.7499
4
4
0.749852
0.749999
0.749999
3
6
通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。
习题2
20．（上机题）Newton迭代法
（1）给定初值
及容许误差
，编制Newton法解方程
根的通用程序。
（2）给定方程
，易知其有三个根
，
，
。
1．由Newton方法的局部收敛性可知存在
，当
时，Newton迭代序列收敛于根
。
试确定尽可能大的
。
2．试取若干初始值，观察当
，
，
，
，
时Newton序列是否收敛以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？
解：
（1）编制的通用程序：
#include
#include
#defineeps0.000001/给定容许误差
floatf（floatx）//定义函数f（x）
{
floatf;
f=x*x*x/3-x;//f（x）的表达式;
return（f）;
}
floatdf（floatx）//定义函数df（x），计算f（x）的导函数
{
floatdf;
df=x*x-1;//f（x）导函数的表达式;
return（df）;
}
voidmain（void）
{
floatx0,x1,a;
intk=0;
cout<<"请输入初值x0:
";
cin>>x0;
do
{
a=-f（x0）/df（x0）;
x1=x0+a;
k++;
x0=x1;
}
while（fabs（a）>eps）;
cout<//输出迭代的次数和根值
}
（2）计算迭代序列收敛于根
的尽可能大的
的函数为：
#include
#include
voiddelay（intn）//定义延时函数
{for（n=10000;n>0;n--）;}
#defineeps0.000001
floatf（floatx）//定义函数f（x）
{
floatf;
f=x*x*x/3-x;//f（x）的表达式;
return（f）;
}
floatdf（floatx）//定义函数df（x），计算f（x）的导函数
{
floatdf;
df=x*x-1;//f（x）导函数的表达式;
return（df）;
}
intjudgement（floatz）
{
intcount=5;
floatx0,x1,type,type1;
x0=z;
while（count-->0）
{
x1=x0-f（x0）/df（x0）;
type=fabs（x1）;
type1=fabs（x1-x0）;//调试值用
cout<<"count="<if（fabs（x1-x0）return1;
x0=x1;
delay（30000）;//调试值用
}
return0;
}
voidmain（void）
{
floatdelta=0;
intflag=1;
while（flag==1）
{
cout<<"方程的根为:
"<<'\n';
delta+=eps;
flag=judgement（delta）;
}
cout<<"输出方程根收敛的区间值:
\n";
cout<}
当步长为0.001时，程序计算出的δ的为δ=0.774，即在区间（-0.774，0.774）内迭代序列收敛于0。
对于不同得初始值收敛于不同的根,
在（-∞，-1）内收敛于
，在（-0.774，0.774）内收敛于
，在（1，+∞）内收敛于
，但在内（0.774，1）和（－1，0.774）均可能收敛于
和
。
，
，
分别为方程的精确解。
分析：
对于不同的初值，迭代序列会收敛于不同的根，所以在某个区间内求根对于初值的选取有很大的关系。
产生上述结果的原因是区间不满足大范围收敛的条件。
习题3
39．（上机题）列主元Gauss消去法
对于某电路的分析，归结为求解线性方程组RI=V。
其中,
31-13000-10000
-1335-90-110000
0-931-1000000
00-1079-30000-9
000-3057-70-50
0000-747-3000
00000-304100
0000-50027-2
000-9000-229
R=
VT=（-15,27,-23,0,-20,12,-7,7,10）T
（1）编制解n阶线性方程组Ax=b的列主元三角分解法的通用程序；
（2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留五位有效数；
（3）本编程之中，你提高了哪些编程能力？
程序为：
#include
#include
voidmain（void）
{
inti,j,n,k,q;
floata[10][11],s[10],s1[10];
cout<<"请输入n的值:
";
cin>>n;
cout<<"输入数组a:
"<for（i=1;i<=n;i++）
for（j=1;j<=（n+1）;j++）
cin>>a[i][j];//给矩阵a赋值
for（i=1;i<=n;i++）
{
for（j=1;j<=（n+1）;j++）
cout<cout<<'\n';
}//输出数组a
cout<<"'''''''''''''''''''''''''"<<'\n';
//进行第一行和第一列元素的求取'''''''''''''''''''''''''//
intt=1;
for（i=1;i<=n;i++）
{s[i]=a[i][1];}
floatmax=fabs（s[1]）;
for（i=2;i<=n;i++）
if（fabs（s[i]）>max）
{
max=fabs（s[i]）;
t=i;
}
for（j=1;j<=（n+1）;j++）
{
floatb=a[1][j];
a[1][j]=a[t][j];
a[t][j]=b;
}//进行第一列主元互换
for（i=2;i<=n;i++）
a[i][1]=a[i][1]/max;//第一列除以a[1][1]
for（i=1;i<=n;i++）
{
for（j=1;j<=（n+1）;j++）
cout<cout<<'\n';
}
//输出进行第一步变换的数组a
cout<<"'''''''''''''''''''''''''"<<'\n';
//进行第k步分解'''''''''''''''''''''''''''''''''''''''''//
for（k=2;k<=n;k++）
{
for（i=k;i<=n;i++）
{
floatsum=0;
for（q=1;qsum+=a[i][q]*a[q][k];
s1[i]=a[i][k]-sum;
}
intl=k;
floatm=fabs（s1[k]）;
for（i=k;i<=n;i++）
//比较第k步分解的第k列值的大小
{
if（fabs（s1[i]）>m）
{
m=fabs（s1[i]）;
l=i;//返回行值
}
}
for（j=1;j<=n+1;j++）//交换两行元素
{
floats2=a[k][j];
a[k][j]=a[l][j];
a[l][j]=s2;
}
for（j=k;j<=n+1;j++）
//算出第k行行元素的值
{
floatsum1=0;
for（q=1;qsum1+=a[k][q]*a[q][j];
a[k][j]=a[k][j]-sum1;
}
for（i=k+1;i<=n;i++）
//算出第k列列元素的值
{
floatsum2=0;
for（q=1;qsum2+=a[i][q]*a[q][k];
a[i][k]=（a[i][k]-sum2）/（a[k][k]）;
}
}//第k步分解结束
for（i=1;i<=n;i++）
{
for（j=1;j<=（n+1）;j++）
cout<cout<<'\n';
}//输出改变后的数组
//输出解'''''''''''''''''''''''''''''''''''''''''''''''''''''''//
floatx[10];
for（i=n-1;i>=1;i--）
{
x[n]=a[n][n+1]/a[n][n];
floatsum3=0;
for（j=i+1;j<=n;j++）
sum3+=a[i][j]*x[j];
x[i]=（a[i][n+1]-sum3）/a[i][i];
}//回代过程
for（i=1;i<=n;i++）
{
cout<<'x'<}//输出解向量
}
结果：
方程的解为：
x1=-0.28923，x2=0.34544，x3=-0.71281，
x4=-0.22061，x5=-0.43040，x6=0.15431，
x7=-0.057823，x8=0.20105，x9=0.29023。
分析：
我感觉是提高了查错误点的能力和编循环语句的能力，即利用很规整的迭代公式进行编程。
另外列主元三角分解法的阶梯步骤有了更深的了解!
习题4
37.（上机题）3次样条插值函数
（1）编制求第一型3次样条插值函数的通用程序;
（2）已知汽车曲线型值点的数据如下：
0
1
2
3
4
5
6
7
8
9
10
2.51
3.30
4.04
4.70
5.22
5.54
5.78
5.40
5.57
5.70
5.80
端点条件为
=0.8，
=0.2。
用所编制程序求车门的3次样条插值函数S（x）,并打印出S（i+0.5）（i=0,1,…9）。
解：
通用程序：
#include
voidmain（void）
{
floatx[11];//存放数组x[j]
floaty[11];//存放数组y[j]
floath[11];//存放数组h[j]
floatu[11];//存放数组u[j]
floatv[11];//存放数组v[j]
floatd[11];//存放数组d[j]
floatM[11];//存放数组M[j]
floatb[11];//存放数组b[j]
floatt[11],l[11],yy[11],s[4],aa1,aa2,aa3,aa4;
floats1[10];
inti,j,n;
floatxx;//x为区间值
//将初值初始化
cout<<"请输入n的值:
\n";
cin>>n;
cout<<"输入数组x:
\n";
for（i=0;i<=n;i++）cin>>x[i];
cout<<"输入数组y:
\n";
for（i=0;i<=n;i++）cin>>y[i];
//输入端点值
floatdf[2];
cout<<"输入两个端点值:
\n";
for（i=0;i<2;i++）cin>>df[i];
//利用书本上的算法求出所需要的值
//求出h[j]的值
for（j=0;j<=n-1;j++）
{
h[j]=x[j+1]-x[j];
cout<<'h'<<'['<}
cout<//求出u[j]和v[j]的初值
v[0]=1;
u[n]=1;
for（j=1;j<=n-1;j++）
{
u[j]=h[j-1]/（h[j-1]+h[j]）;
v[j]=h[j]/（h[j-1]+h[j]）;
}
//求出d[j]的值
for（j=1;j{d[j]=6*（（y[j+1]-y[j]）/h[j]-（y[j]-y[j-1]）/h[j-1]）/（h[j]+h[j-1]）;}
d[0]=6*（（y[1]-y[0]）/h[0]-df[0]）/h[0];
d[n]=6*（df[1]-（y[n]-y[n-1]）/h[n-1]）/h[n-1];
for（j=1;j<=n;j++）
{
cout<<'u'<<'['<}
cout<for（j=0;j{
cout<<'v'<<'['<}
cout<for（j=0;j<=n;j++）
{
cout<<'d'<<'['<}
cout<//利用书本上的追赶法求解方程组
for（i=0;i<=n;i++）
{b[i]=2;}
cout<t[0]=b[0];
yy[0]=d[0];
//消元过程
for（i=1;i<=n;i++）
{
l[i]=u[i]/t[i-1];
t[i]=b[i]-l[i]*v[i-1];
yy[i]=d[i]-l[i]*yy[i-1];
}
//回代过程
M[n]=yy[n]/t[n];
for（i=n-1;i>=0;i--）
{
M[i]=（yy[i]-v[i]*M[i+1]）/t[i];
}
//将M[j]的值输出
for（i=0;i<=n;i++）
{cout<<'M'<<'['<//输出插值多项式的系数
for（j=0;j{
s[0]=y[j];
s[1]=（y[j+1]-y[j]）/h[j]-（M[j]/3+M[j+1]/6）*h[j];
s[2]=M[j]/2;
s[3]=（M[j+1]-M[j]）/（6*h[j]）;
cout<<"当x的值在区间"<<'x'<<'['<\n";
for（intk=0;k<4;k++）
{
cout<<'s'<<'['<}
cout<}
}
（2）编制的程序求车门的3次样条插值函数S（x）：
x属于区间[0,1]时；S（x）=2.51+0.8（x）-0.0014861（x）（x）-0.00851395（x）（x）（x）
x属于区间[1,2]时；S（x）=3.3+0.771486（x-1）-0.027028（x-1）（x-1）-0.00445799（x-1）（x-1）（x-1）
x属于区间[2,3]时；S（x）=4.04+0.704056（x-2）-0.0404019（x-2）（x-2）-0.0036543（x-2）（x-2）（x-2）
x属于区间[3,4]时；S（x）=4.7+0.612289（x-3）-0.0513648（x-3）（x-3）-0.0409245（x-3）（x-3）（x-3）
x属于区间[4,5]时；S（x）=5.22+0.386786（x-4）-0.174138（x-4）（x-4）+0.107352（x-4）（x-4）（x-4）
x属于区间[5,6]时；S（x）=5.54+0.360567（x-5）+0.147919（x-5）（x-5）-0.268485（x-5）（x-5）（x-5）
x属于区间[6,7]时；S（x）=5.78-0.149051（x-6）-0.657537（x-6）（x-6）+0.426588（x-6）（x-6）（x-6）
x属于区间[7,8]时；S（x）=5.4-0.184361（x-7）+0.622227（x-7）（x-7）-0.267865（x-7）（x-7）（x-7）
x属于区间[8,9]时；S（x）=5.57+0.256496（x-8）-0.181369（x-8）（x-8）+0.0548728（x-8）（x-8）（x-8）
x属于区间[9,10]时；S（x）=5.7+0.058376（x-9）-0.0167508（x-9）（x-9）+0.0583752（x-9）（x-9）（x-9）
S（0.5）=2.90856S（1.5）=3.67843S（2.5）=4.38147
S（3.5）=4.98819S（4.5）=5.38328S（5.5）=5.7237
S（6.5）=5.59441S（7.5）=5.42989S（8.5）=5.65976
S（9.5）=5.7323
习题6
23．（上机题）常微分方程初值问题数值解
（1）编制RK4方法的通用程序；
（2）编制AB4方法的通用程序（由RK4提供初值）；
（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；
（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;
（5）对于初值问题
取步长
应用
（1）~（4）中的四种方法进行计算，并将计算结果和精确解
作比较；
（6）通过本上机题，你能得到哪些结论？
解：
程序为：
#include
#include
#include
#include
ofstreamoutfile（"data.txt"）;
//此处定义函数f（x,y）的表达式
//用户可以自己设定所需要求得函数表达式
doublef1（doublex,doubley）
{
doublef1;
f1=（-1）*x*x*y*y;
returnf1;
}
//此处定义求函数精确解的函数表达式
doublef2（doublex）
{
doublef2;
f2=3/（1+x*x*x）;
returnf2;
}
//此处为精确求函数解的通用程序
voidaccurate（doublea,doubleb,doubleh）
{
doublex[100],accurate[100];
x[0]=a;
inti=0;
outfile<<"输出函数准确值的程序结果:
\n";
do{
x[i]=x[0]+i*h;
accurate[i]=f2（x[i]）;
outfile<<"accurate["<i++;
}while（i<（b-a）/h+1）;
}
//此处为经典Runge-Kutta公式的通用程序
voidRK4（doublea,doubleb,doubleh,doublec）
{
inti=0;
doublek1,k2,k3,k4;
doublex[100],y[100];
y[0]=c;
x[0]=a;
outfile<<"输出经典Runge-Kutta公式的程序结果:
\n";
do
{
x[i]=x[0]+i*h;
k1=f1（x[i],y[i]）;
k2=f1（（x[i]+h/2）,（y[i]+h*k1/2））;
k3=f1（（x[i]+h/2）,（y[i]+h*k2/2））;
k4=f1（（x[i]+h）,（y[i]+h*k3））;
y[i+1]=y[i]+h*（k1+2*k2+2*k3+k4）/6;
outfile<<"y"<<"["<i++;
}while（i<（b-a）/h+1）;
}
//此处为4阶Adams显式方法的通用程序
voidAB4（doublea,doubleb,doubleh,doublec）
{
doublex[100],y[100],y1[100];
doublek1,k2,k3,k4;
y[0]=c;
x[0]=a;
outfile<<"输出4阶Adams显式方法的程序结果:
\n";
for（inti=0;i<=2;i++）
{
x[i]=x[0]+i*h;
k1=f1（x[i],y[i]）;
k2=f1（（x[i]+h/2）,（y[i]+h*k1/2））;
k3=f1（
升级会员 