计算方法程序集及例程.docx
- 文档编号:12798129
- 上传时间:2023-06-08
- 格式:DOCX
- 页数:45
- 大小:222.65KB
计算方法程序集及例程.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;j if(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;i f[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;j m[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;j printf("(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(count return0; } 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 程序 例程