计算方法程序集及例程.docx
《计算方法程序集及例程.docx》由会员分享，可在线阅读，更多相关《计算方法程序集及例程.docx（45页珍藏版）》请在冰豆网上搜索。
计算方法程序集及例程
《计算方法》程序集及例程
（C语言部分）
研究生《计算方法》课程学习辅助教材
江西理工大学理学院
实验一：
拉格朗日插值多项式
命名（源程序.cpp及工作区.dsw）：
lagrange
问题：
0.4
0.55
0.80
0.9
1
0.41075
0.57815
0.88811
1.02652
1.17520
用四次拉格朗日多项式求
的函数近似值
//Lagrange.cpp
#include
#include
#defineN4
intcheckvalid（doublex[],intn）;
voidprintLag（doublex[],doubley[],doublevarx,intn）;
doubleLagrange（doublex[],doubley[],doublevarx,intn）;
voidmain（）
{doublex[N+1]={0.4,0.55,0.8,0.9,1};
doubley[N+1]={0.41075,0.57815,0.88811,1.02652,1.17520};
doublevarx=0.5;
if（checkvalid（x,N）==1）
{printf（"\n\n插值结果:
P（%f）=%f\n",varx,Lagrange（x,y,varx,N））;}
else
{printf（"结点必须互异"）;}
getch（）;}
intcheckvalid（doublex[],intn）
{
inti,j;
for（i=0;i{
for（j=i+1;j{
if（x[i]==x[j]）//若出现两个相同的结点，返回-1
{
return-1;
}
}
}
return1;
}
doubleLagrange（doublex[],doubley[],doublevarx,intn）
{
doublefenmu;
doublefenzi;
doubleresult=0;
inti,j;
printf（"Ln（x）=\n"）;
for（i=0;i{
fenmu=1;
for（j=0;j{
if（i!
=j）
{
fenmu=fenmu*（x[i]-x[j]）;
}
}
printf（"\t%f",y[i]/fenmu）;
fenzi=1;
for（j=0;j{
if（i!
=j）
{
printf（"*（x-%f）",x[j]）;
fenzi=fenzi*（varx-x[j]）;
}
}
if（i!
=n）
{
printf（"+\n"）;
}
result+=y[i]/fenmu*fenzi;
}
returnresult;
}运行及结果显示：
实验二：
牛顿插值多项式
命名（源程序.cpp及工作区.dsw）：
newton_cz
问题：
0.4
0.55
0.80
0.9
1
0.41075
0.57815
0.88811
1.02652
1.17520
用牛顿插值多项式求
//newton_cz.cpp
#include
#include
#include
#defineN4
intcheckvaild（doublex[],intn）
{
inti,j;
for（i=0;i{
for（j=i+1;jif（x[i]==x[j]）
return-1;
}
return1;
}
voidchashang（doublex[],doubley[],doublef[][N+1]）
{
inti,j,t=0;
for（i=0;i{
f[i][0]=y[i];//f[0][0]=y[0],f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];f[4][0]=y[4]
}
for（j=1;j{
for（i=0;if[i][j]=（f[i+1][j-1]-f[i][j-1]）/（x[t+i+1]-x[i]）;//一阶为f[0][1]~f[3][1];二阶为f[0][2]~f[2][2]
//三阶为f[0][3]~f[1][3];四阶为f[0][4]
t++;
}
}
doublecompvalue（doublet[][N+1],doublex[],doublevarx）
{
inti,j,r=0;
doublesum=t[0][0],m[N+1]={sum,1,1,1,1};
printf（"i\tXi\tF（Xi）\t1阶\t2阶\t\t3阶\t4阶\n"）;
printf（"--------------------------------------------------------------------------------"）;
for（i=0;i{
r=i;
printf（"x%d\t%f",i,x[i]）;
for（j=0;j<=i;j++）
{
printf（"%f",t[r][j]）;
r--;
}
printf（"\n"）;
}
printf（"--------------------------------------------------------------------------------"）;/**/
printf（"N（x）=\n"）;
printf（"%f\n",t[0][0]）;
for（i=1;i{
for（j=0;jm[i]=m[i]*（varx-x[j]）;
m[i]=t[0][i]*m[i];
sum=sum+m[i];
}
for（i=1;i{printf（"+%f*",t[0][i]）;
for（j=0;jprintf（"（x-%f）",x[j]）;
printf（"\n"）;
}
returnsum;
}
voidmain（）
{
doublevarx,f[N+1][N+1];
doublex[N+1]={0.4,0.55,0.8,0.9,1};
doubley[N+1]={0.41075,0.57815,0.88811,1.02652,1.17520};
checkvaild（x,N）;
chashang（x,y,f）;
varx=0.5;
if（checkvaild（x,N）==1）
{
chashang（x,y,f）;
printf（"\n\n牛顿插值结果:
P（%f）=%f\n",varx,compvalue（f,x,varx））;
}
else
printf（"输入的插值节点的x值必须互异!
"）;
getch（）;
}
运行及结果显示：
实验三：
自适应梯形公式
命名（源程序.cpp及工作区.dsw）：
autotrap
问题：
计算
的近似值，使得误差（EPS）不超过
//autotrap.cpp
#include
#include
#include
#defineEPS1e-6
doublef（doublex）
{
return4/（1+x*x）;
}
doubleAutoTrap（double（*f）（double）,doublea,doubleb,doubleeps）
{
doublet2,t1,sum=0.0;
inti,k;
t1=（b-a）*（f（a）+f（b））/2;
printf（"T（%4d）=%f\n",1,t1）;
for（k=1;k<11;k++）
{
for（i=1;i<=pow（2,k-1）;i++）
sum+=f（a+（2*i-1）*（b-a）/pow（2,k））;
t2=t1/2+（b-a）*sum/pow（2,k）;
printf（"T（%4d）=%f\n",（int）pow（2,k）,t2）;
if（fabs（t2-t1）break;
else
t1=t2;
sum=0.0;
}
returnsum;
}
voidmain（）
{
doubles;
s=AutoTrap（f,0.0,1.0,EPS）;
getch（）;
}
运行及结果显示：
实验四：
龙贝格算法
命名（源程序.cpp及工作区.dsw）：
romberg
问题：
求
的近似值，要求误差不超过
//romberg.cpp
#include
#include
#include
#defineN20
#defineEPS1e-6
doublef（doublex）{return4/（1+x*x）;}
doubleRomberg（doublea,doubleb,double（*f）（double）,doubleeps）
{
doubleT[N][N],p,h;
intk=1,i,j,m=0;
T[0][0]=（b-a）/2*（f（a）+f（b））;
do{
p=0;
h=（b-a）/pow（2,k-1）;
for（i=1;i<=pow（2,k-1）;i++）
p=p+f（a+（2*i-1）*h/2）;
T[0][k]=T[0][k-1]/2+p*h/2;
for（i=1;i<=k;i++）
{j=k-i;
T[i][j]=（pow（4,i）*T[i-1][j+1]-T[i-1][j]）/（pow（4,i）-1）;
}
k++;
p=fabs（T[k-1][0]-T[k-2][0]）;
}while（p>=EPS）;
k--;
while（m<=k）
{for（i=0;i<=m;i++）
printf（"T（%d,%d）=%f",i,m-i,T[i][m-i]）;
m++;
printf（"\n"）;
}
returnT[k][0];
}
voidmain（）
{
printf（"\n\n积分结果= %f",Romberg（0,1,f,EPS））;
getch（）;
}
运行及结果显示：
实验五：
牛顿切线法
问题：
求方程
在
，
附近的根（精度
）
//newton_qxf.cpp
#include
#include
#include
#defineMaxK1000
#defineEPS0.5e-3
doublef（doublex）
{
returnx*x*x-x-1;
}
doublef1（doublex）
{
return3*x*x-1;
}
intnewton（double（*f）（double）,double（*f1）（double）,double&x0,doubleeps）
{
doublexk,xk1;
intcount=0;
printf（"k\txk\n"）;
printf（"-----------------------\n"）;
xk=x0;
printf（"%d\t%f\n",count,xk）;
do
{
if（（*f1）（xk）==0）
return2;
xk1=xk-（*f）（xk）/（*f1）（xk）;
if（fabs（xk-xk1）{
count++;
xk=xk1;
printf（"%d\t%f\n",count,xk）;
x0=xk1;
return1;
}
count++;
xk=xk1;
printf（"%d\t%f\n",count,xk）;
}while（countreturn0;
}
voidmain（）
{
for（inti=0;i<2;i++）
{
doublex=0.6;
if（i==1）
x=1.5;
printf（"x0初值为%f:
\n",x）;
if（newton（f,f1,x,EPS）==1）
{
printf（"-----------------------\n"）;
printf（"therootisx=%f\n\n\n",x）;
}
else
{
printf（"themethodisfail!
"）;
}
}
getch（）;
}
实验六：
牛顿下山法
命名（源程序.cpp及工作区.dsw）：
newton_downhill
问题：
求方程
在
，
附近的根（精度
）
//newton_downhill.cpp
#include
#include
#include
#include
#defineEt1e-3//下山因子下界
#defineE11e-3//根的误差限
#defineE21e-3//残量精度
doublef（doublex）
{
returnx*x*x-x-1;
}
doublef1（doublex）
{
return3*x*x-1;
}
voiderrormess（intb）
{
char*mess;
switch（b）
{
case-1:
mess="f（x）的导数为0!
";
break;
case-2:
mess="下山因子已越界，下山处理失败";
break;
default:
mess="其他类型错误!
";
}
printf（"themethodhasfaild!
because%s",mess）;
}
intNewton（double（*f）（double）,double（*f1）（double）,double&x0）
{
intk=0;
doublet;
doublexk,xk1;
doublefxk,fxk_temp;
printf（"ktxkf（xk）\n"）;
printf（"----------------------------------------------------------\n"）;
printf（"%-20d",k）;
xk=x0;
printf（"%-15f",x0）;
fxk=（*f）（xk）;
printf（"%-20f",fxk）;
printf（"\n"）;
for（k=1;;k++）
{
t=1;
while
（1）
{
printf（"%-10d",k）;
printf（"%-10f",t）;
if（（*f1）（xk）!
=0）
{
xk1=xk-t*（*f）（xk）/（*f1）（xk）;
}
else
{
return-1;
}
printf（"%-15f",xk1）;
fxk_temp=（*f）（xk1）;
printf（"%-20f",fxk_temp）;
if（fabs（fxk_temp）>fabs（fxk））
{
t=t/2;
printf（"\n"）;
if（t{
return-2;
}
}
else
{
printf（"下山成功\n"）;
break;
}
}
if（fabs（xk-xk1）{
x0=xk1;
return1;
}
xk=xk1;
}
}
voidmain（）
{
intb;
doublex0;
x0=0.6;
b=Newton（f,f1,x0）;
if（b==1）
printf（"\ntherootx=%f\n",x0）;
else
errormess（b）;
getch（）;
}
运行及结果显示：
实验七：
埃特金加速算法
命名（源程序.cpp及工作区.dsw）：
aitken
问题：
求方程
在
，
附近的根（精度
）
//aitken.cpp
#include
#include
#include
#defineMaxK100
#defineEPS0.5e-3
doubleg（doublex）
{
returnx*x*x-1;
}
intaitken（double（*g）（double）,double&x,doubleeps）
{
intk;
doublexk=x,yk,zk,xk1;
printf（"kxkykzkxk+1\n"）;
printf（"-------------------------------------------------------------------\n"）;
for（k=0;k{
yk=（*g）（xk）;
zk=（*g）（yk）;
xk1=xk-（yk-xk）*（yk-xk）/（zk-2*yk+xk）;
printf（"%-10d%-15f%-15f%-15f%-15f\n",k,xk,yk,zk,xk1）;
if（fabs（yk-xk）<=eps）
{
x=xk1;
returnk+1;
}
xk=xk1;
}
return-1;
}
voidmain（）
{
doublex=1.5;
intk;
k=aitken（g,x,EPS）;
if（k==-1）
printf（"迭代次数越界!
\n"）;
else
printf（"\n经k=%d次迭代，所得方程根为:
x=%f\n",k,x）;
getch（）;
}
运行及结果显示：
实验八：
正割法
问题：
求方程
在
附近的根（精度0.5e-8）
//ZhengGe.cpp
#include
#include
#include
#defineMaxK1000
#defineEPS0.5e-8
doublef（doublex）
{
returnx*x*x-x-1;
}
intZhengGe（double（*f）（double）,doublex0,double&x1,doubleeps）
{
printf（"kxkf（xk）\n"）;
printf（"---------------------------------------------\n"）;
intk;
doublexk0,xk,xk1;
xk0=x0;
printf（"%-10d%-15f%-15f\n",0,x0,（*f）（x0））;
xk=x1;
for（k=1;k{
if（（*f）（xk）-（*f）（xk0）==0）
return-1;
xk1=xk-（*f）（xk）/（（*f）（xk）-（*f）（xk0））*（xk-xk0）;
printf（"%-10d%-15f%-15f\n",k,xk,（*f）（xk））;
if（fabs（xk1-xk）<=eps）
{
printf（"%-10d%-15f%-15f\n\n",k+1,xk1,（*f）（xk1））;
printf（"---------------------------------------------\n\n"）;
x1=xk1;
return1;
}
xk0=xk;
xk=xk1;
}
return-2;
}
voidmain（）
{
doublex0=1,x1=2;
if（ZhengGe（f,x0,x1,EPS）==1）
{
printf（"therootisx=%f\n",x1）;
}
else
{
printf（"themethodisfail!
"）;
}
getch（）;
}
实验九：
高斯列选主元算法
命名（源程序.cpp及工作区.dsw）：
colpivot
问题：
求解方程组并计算系数矩阵行列式值
//colpivot.cpp
#include
#in
升级会员 