数值方法课后上机.docx
《数值方法课后上机.docx》由会员分享，可在线阅读，更多相关《数值方法课后上机.docx（27页珍藏版）》请在冰豆网上搜索。
数值方法课后上机
第二章插值法
1、题目
编制通用程序，对n+1个节点
（1）n次拉格朗日插值计算公式
（2）n次牛顿向前插值计算公式
（3）n次牛顿向后插值计算公式
已知
，取
，
，
……，10。
用通用程序
（1）、
（2）、（3）计算ln1.54及ln1.98的近似值。
2.拉格朗日插值法的计算流程图、源程序以及计算结果
（1）流程图
图1拉格朗日插值计算公式流程图
（2）源程序
#include
#include
main（）
{floata=1.0;//左端点
floatb=2.0;//右端点
floatJ=0.1;//步长
floatdaiqiu=1.54;//待求的值
intNUM;
inti,k=0;
floatx[100],y[100];
floatm,l[100],L=0,fenzi,fenmu,result;
NUM=（int）（（b-a）/J+1）;
for（m=a;m<=b;m=m+J,k++）
{x[k]=m;
y[k]=log（m）;
}
for（i=0;i<=NUM-1;i++）
{
fenzi=1;
fenmu=1;
for（k=0;k<=NUM-1;k++）
{
if（k==i）
continue;
else
{
fenzi=fenzi*（daiqiu-x[k]）;
fenmu=fenmu*（x[i]-x[k]）;
}
}
l[i]=fenzi/fenmu;
}
for（i=0;i<=NUM-1;i++）
{
L=L+l[i]*y[i];
}
printf（"result=%f",L）;
}
（3）输出结果
ln1.54
ln1.98
3.牛顿向前差值的计算流程图、源程序以及计算结果
（1）流程图
图2牛顿向前插值计算公式流程图
（2）源程序
#include
#include
main（）
{floata=1;//左端点
floatb=2;//右端点
floatj=0.1;//步长
floatX=1.94;//待求值
intNUM;
floatI,t,T[100],x[100],y[100],N[100],delt_y[100][100],delt_y1[100],Niu,result;
inti,k=0,n;
NUM=（int）（（b-a）/j+1）;
t=（X-a）/j;
for（I=a;I<=b+j;I=I+j,k++）
{x[k]=I;
y[k]=log（I）;//求节点对应函数值
}
for（i=0;i<=NUM-2;i++）
{for（k=0;k<=NUM-2-i;k++）
{delt_y[i][k]=y[k+1]-y[k];
y[k]=delt_y[i][k];//构造向前差分表
}
}
for（i=0;i<=NUM-2;i++）
{k=1;//用于阶乘
T[i]=1;
for（n=1;n<=i+1;n++）
{k=k*n;
T[i]=T[i]*（t-n+1）;}
T[i]=T[i]/k;
}
Niu=y[0];
for（i=0;i<=NUM-2;i++）
{Niu=Niu+（float）（T[i]*delt_y[i][0]）;}
printf（"result=%f\n",Niu）;
}
（3）输出结果
ln1.54
ln1.98
4.牛顿向后差值法的计算流程图、源程序以及计算结果
（1）流程图
图3牛顿向后插值计算公式流程图
（2）源程序
#include
#include
#definea（float）1//左端点
#defineb（float）2//右端点
#definej0.1//步长
#defineX1.98//待求值
#defineNUM（int）（（b-a）/j+1）
main（）
{floatI,t,T[NUM],x[NUM],y[NUM],N[NUM],delt_y[NUM][NUM],delt_y1[NUM],Niu,result;
inti,k=0,n;
t=（X-b）/j;
for（I=b;I>=a-j;I=I-j,k++）//用“a-j”
{x[k]=I;
y[k]=log（I）;//求节点对应函数值
}
Niu=y[0];//向后差值公式中的yn
for（i=0;i<=NUM-2;i++）
{for（k=0;k<=NUM-2-i;k++）
{delt_y[i][k]=y[k]-y[k+1];
y[k]=delt_y[i][k];//构造向后差分表
}
}
for（i=0;i<=NUM-2;i++）
{k=1;//用于阶乘
T[i]=1;
for（n=1;n<=i+1;n++）
{k=k*n;
T[i]=T[i]*（t+n-1）;}
T[i]=T[i]/k;
}
for（i=0;i<=NUM-2;i++）
{Niu=Niu+（T[i]*delt_y[i][0]）;}
printf（"\nln（%f）=%f\n",X,Niu）;
}
（3）输出结果
ln1.54
ln1.98
结果分析：
有以上结果可以看出，差值点函数值能较好的与真实值相吻合。
但编程实现的难易程度以及对计算机内存的使用程度是不同的。
第五章解线性方程组的直接法
1、题目
编制解
通用子程序
（1）列主元消去法；
（2）改进平方根法；
用列主元消去法程序解方程组
用改进平方根法程序解方程组
2.列主元消去法计算流程图、源程序及输出结果
（1）流程图
图4列主元消去法流程图
（2）源程序
#include
#include
#defineA{{1.1348,3.8326,1.1651,3.4017},{0.5301,1.7875,2.5330,1.5435},{3.4129,4.9317,8.7643,1.3142},{1.2371,4.9998,10.6721,0.0147}}//系数矩阵
#defineB{9.5342,6.3941,18.4231,16.9237}//右端列向量
#defineR4//系数矩阵行数
#defineL4//系数矩阵列数
main（）
{floatxs[R][L]=A;
floatb[R]=B;
floatzg[R][L+1],temp,X[R],m;//增广矩阵、temp用于记录主元及交换
inti=0,j=0,k,r,l;//r、l分别用于记录主元坐标,i、j、k分别用于记录循环次数
for（i=0;i<=R-1;i++）
{for（j=0;j<=L-1;j++）
zg[i][j]=xs[i][j];}
for（i=0;i<=R-1;i++）
zg[i][L]=b[i];
for（i=0;i<=R-2;i++）
{temp=abs（zg[i][i]）;
r=i;
for（j=i+1;j<=R-1;j++）
{if（tempr=j;//选取列主元
}
for（k=i;k<=L;k++）
{
temp=zg[i][k];
zg[i][k]=zg[r][k];
zg[r][k]=temp;
}//交换行
for（j=i+1;j<=R-1;j++）
{m=zg[j][i]/zg[i][i];
for（k=i;k<=L;k++）
zg[j][k]=zg[j][k]-zg[i][k]*m;}//消元
for（k=0;k<=R-1;k++）
for（j=0;j<=L;j++）
{printf（"%f",zg[k][j]）;
if（j==L）
printf（"\n"）;//输出中间过程矩阵
}
printf（"\n"）;
}
X[R-1]=zg[R-1][L]/zg[R-1][L-1];
for（i=R-2;i>=0;i--）
{temp=0;
{for（j=i+1;j<=L-1;j++）
temp=temp+X[j]*zg[i][j];
X[i]=（zg[i][L]-temp）/zg[i][i];}//回代
}
for（i=0;i<=R-1;i++）
printf（"X[%d]=%f\n",i+1,X[i]）;}
（3）用列主元消去法解方程组
（1）输出结果
3.改进平方根法计算流程图、源程序及输出结果
（1）流程图
图5改进平方根法流程图
（2）源程序
#include
#include
#definexs{{4,2.4,2,3},{2.4,5.44,4,5.8},{2,4,5.21,7.45},{3,5.8,7.45,19.66}}//系数矩阵
#defineB{12.280,16.928,22.957,50.945}//右端向量
#definer4//系数矩阵的行数
#definel4//系数矩阵的列数
main（）
{floatA[r][l]=xs,b[r]=B;//A：
系数矩阵，b:
右端向量
floatL[r][l]={0},D[r][l]={0},y[r],x[r],tl;//L、D分别用于存储分解式中的矩阵
inti,ii,j,jj,k,m;//用于计数
//以下为分解计算
D[0][0]=A[0][0];
for（i=0;iL[i][0]=A[i][0]/D[0][0];
for（i=0;iL[i][i]=1;
for（i=1;ifor（j=i;j{tl=0;
for（k=0;k{tl=tl+L[j][k]*D[k][k]*L[i][k];}
if（j==i）
D[i][i]=A[i][i]-tl;
else
L[j][i]=（A[j][i]-tl）/D[i][i];
}
//以下输出分解得到的矩阵
printf（"D=\n"）;
for（ii=0;ii{for（jj=0;jjprintf（"%f",D[ii][jj]）;
if（jj==r）
printf（"\n"）;
}
printf（"\n"）;
printf（"L=\n"）;
for（ii=0;ii{for（jj=0;jjprintf（"%f",L[ii][jj]）;
if（jj==r）
printf（"\n"）;
}
//解Ly=b
y[0]=b[0]/L[0][0];
for（i=1;i{tl=0;
for（j=0;jtl=tl+L[i][j]*y[j];
y[i]=b[i]-tl;
}
//解L'x=inv（D）y
x[r-1]=y[r-1]/D[r-1][r-1];
for（i=r-2;i>=0;i--）
{tl=0;
for（j=r-1;j>i;j--）
tl=tl+x[j]*L[j][i];
x[i]=y[i]/D[i][i]-tl;
}
for（i=0;iprintf（"x[%d]=%f\n",i+1,x[i]）;
}
（3）输出结果
用改进平方根法程序
（2）解方程组
结果分析
直接法对方程组的求解过程在不考虑计算机精度误差下是能实现方程组根的准确算，都包含两个过程：
消去过程、回代过程；改进平方根法的计算量和对计算机内存的使用相对较小，但更容易产生误差。
第六章解线性方程组的迭代法
1、题目
编制通用子程序
（1）雅克比迭代格式；
（2）高斯-塞德尔迭代格式；
（3）SOR方法格式；
对方程组
，
取初值
，要求
（1）用雅克比方法计算；
（2）用高斯-塞德尔方法计算；
（3）用SOR迭代法计算当
时的解，最后输出近似解及迭代次数k。
2.雅克比迭代格式计算流程图、源程序及输出结果
（1）流程图
（2）源程序
#include
#include
#definexs{{4,-1,0,-1,0,0},{-1,4,-1,0,-1,0},{0,-1,4,0,0,-1},{-1,0,0,4,-1,0},{0,-1,0,-1,4,-1},{0,0,-1,0,-1,4}}//系数矩阵
#defineB{0,5,0,6,-2,6}//右端向量
#defineX{1,1,1,1,1,1}//初值
#defineyuzhi1e-5//允许误差
#defineR6//系数矩阵行数
#defineL6//系数矩阵列数
main（）
{floatx_0[R]=X,x_1[R],A[R][L]=xs,b[R]=B,duijiao[R],rep=0,fanshu;
inti,j,k;
do
{for（i=0;i{rep=0;
for（j=0;j{if（j==i）
continue;
rep=-A[i][j]/A[i][i]*x_0[j]+rep;}
x_1[i]=rep+b[i]/A[i][i];
}
rep=0;
for（i=0;i{rep=rep+pow（（x_1[i]-x_0[i]）,2）;
x_0[i]=x_1[i];}
for（i=0;iprintf（"x%d=%f",i+1,x_0[i]）;
fanshu=rep;
printf（"\n"）;
}
while（sqrt（fanshu）>yuzhi）;
for（i=0;iprintf（"\nx%d=%f\n",i+1,x_1[i]）;
}
（3）输出结果
3.高斯-塞德尔迭代格式计算流程图、源程序及输出结果
（1）流程图
（2）源程序
#include
#include
#definexs{{4,-1,0,-1,0,0},{-1,4,-1,0,-1,0},{0,-1,4,0,0,-1},{-1,0,0,4,-1,0},{0,-1,0,-1,4,-1},{0,0,-1,0,-1,4}}//系数矩阵
#defineB{0,5,0,6,-2,6}//右端向量
#defineX{1,1,1,1,1,1}//初值
#defineyuzhi1e-5//允许误差
#defineR6//系数矩阵行数
#defineL6//系数矩阵列数
main（）
{floatx_0[R]=X,x_1[R],x_2[R],A[R][L]=xs,b[R]=B,duijiao[R],rep=0,fanshu;
inti,j,k;
//以下为迭代过程
do
{for（i=0;i{rep=0;
for（j=0;j{if（j==i）
continue;
rep=-A[i][j]/A[i][i]*x_0[j]+rep;}
x_1[i]=rep+b[i]/A[i][i];//得到的x的新值
x_2[i]=x_0[i];//保留上次得到的值
x_0[i]=x_1[i];//实现高斯赛德尔得带格式中的要求
}
rep=0;
for（i=0;i{rep=rep+pow（（x_1[i]-x_2[i]）,2）;
x_0[i]=x_1[i];}
for（i=0;iprintf（"x%d=%f",i+1,x_0[i]）;
fanshu=rep;
printf（"\n"）;
}
while（sqrt（fanshu）>yuzhi）;
for（i=0;iprintf（"\nx%d=%f\n",i+1,x_1[i]）;
}
（3）输出结果
4．SOR方法的计算流程图、源程序及输出结果
（1）流程图
（2）源程序
#include
#include
#definexs{{4,-1,0,-1,0,0},{-1,4,-1,0,-1,0},{0,-1,4,0,0,-1},{-1,0,0,4,-1,0},{0,-1,0,-1,4,-1},{0,0,-1,0,-1,4}}//系数矩阵
#defineB{0,5,0,6,-2,6}//右端向量
#defineX{1,1,1,1,1,1}//初值
#definew1.3
#defineyuzhi1e-5//允许误差
#defineR6//系数矩阵行数
#defineL6//系数矩阵列数
main（）
{floatx_0[R]=X,x_1[R],x_2[R],A[R][L]=xs,b[R]=B,duijiao[R],rep=0,fanshu;
inti,j,k;
//以下为迭代过程
do
{for（i=0;i{rep=0;
for（j=0;j{
rep=-A[i][j]/A[i][i]*x_0[j]+rep;}
x_1[i]=x_0[i]+w*（rep+b[i]/A[i][i]）;//利用w,得到的x的新值
x_2[i]=x_0[i];//保留上次得到的值
x_0[i]=x_1[i];//实现迭带格式中的要求
}
rep=0;
for（i=0;i{rep=rep+pow（（x_1[i]-x_2[i]）,2）;
x_0[i]=x_1[i];}
for（i=0;iprintf（"x%d=%f",i+1,x_0[i]）;
fanshu=rep;
printf（"\n"）;
}
while（sqrt（fanshu）>yuzhi）;
for（i=0;iprintf（"\nx%d=%f\n",i+1,x_1[i]）;
}
（3）输出结果
结果分析：
不同的迭代格式对同一个方程组的求解结果不同，在满足相应结果的计算精度下，需要迭代的次数是不同的；即使是用同一种方法，选取的因子不同（SOR中的松弛因子ω）求解结果也是不同，在满足相应结果的计算精度下，需要迭代的次数是不同的；由结果分析可知当ω（0<ω<2）接近2时，迭代次数太多，已经不能得出方程组的近似解。
从收敛速度上看，选取合适松弛因子的SOR比GS快，GS比Jacobi的快。
松弛因子w的选取对于SOR的影响表现在w>=2时，SOR不收敛，并且在这里w=1.1时，收敛效果最好，比GS的收敛速度更快。
w趋向0和2时收敛速度变慢。
迭代法的优点是占有存储空间少，程序实现简单，尤其适用于大型稀疏矩阵，不足之处是要判断迭代是否收敛；相比之下线性方程组的高斯列主元消去法可以解决一般的线性方程组，但舍入误差对结果的影响不容忽视，对于高阶线性方程组容易受计算机容量的炼制，所以适用于中小型线性方程组。
二者各有优缺点。
升级会员 