计算方法程序集及例程.docx

1、计算方法程序集及例程 计算方法程序集及例程 （C语言部分）研究生计算方法课程学习辅助教材江 西 理 工 大 学 理 学 院实验一：拉格朗日插值多项式命名（源程序.cpp及工作区.dsw）：lagrange 问题：0.40.550.800.910.410750.578150.888111.026521.17520 用四次拉格朗日多项式求的函数近似值/Lagrange.cpp#include #include #define N 4int checkvalid(double x, int n);void printLag (double x, double y, double varx, int

2、n);double Lagrange(double x, double y, double varx, int n);void main ()double xN+1 = 0.4, 0.55, 0.8, 0.9, 1; double yN+1 = 0.41075, 0.57815, 0.88811, 1.02652, 1.17520; double varx = 0.5; if (checkvalid(x, N) = 1) printf(nn插值结果: P(%f)=%fn, varx, Lagrange(x, y, varx, N); else printf(结点必须互异); getch();i

3、nt checkvalid (double x, int n) int i,j; for (i = 0; i n; i+) for (j = i + 1; j n+1; j+) if (xi = xj)/若出现两个相同的结点，返回-1 return -1; return 1;double Lagrange (double x, double y, double varx, int n) double fenmu; double fenzi; double result = 0; int i,j; printf(Ln(x) =n); for (i = 0; i n+1; i+) fenmu =

4、1; for (j = 0; j n+1; j+) if (i != j) fenmu = fenmu * (xi - xj); printf(t%f, yi / fenmu); fenzi = 1; for (j = 0; j n+1; j+) if (i != j) printf(*(x-%f), xj); fenzi = fenzi * (varx - xj); if (i != n) printf(+n); result += yi / fenmu * fenzi; return result;运行及结果显示： 实验二：牛顿插值多项式命名（源程序.cpp及工作区.dsw）：newton

5、_cz问题：0.40.550.800.910.410750.578150.888111.026521.17520 用牛顿插值多项式求 /newton_cz.cpp#include#include#include#define N 4int checkvaild(double x,int n) int i,j; for(i=0;in+1;i+) for(j=i+1;jn+1;j+) if(xi=xj) return -1; return 1;void chashang(double x,double y,double fN+1) int i,j,t=0; for(i=0;iN+1;i+) fi0

6、=yi;/f00=y0,f10=y1;f20=y2;f30=y3;f40=y4 for(j=1;jN+1;j+)/ 阶数j for(i=0;iN-j+1;i+) /差商个数i fij=(fi+1j-1-fij-1)/(xt+i+1-xi);/一阶为f01f31;二阶为f02f22 /三阶为f03f13;四阶为f04 t+; double compvalue(double tN+1,double x,double varx) int i,j,r=0; double sum=t00,mN+1=sum,1,1,1,1; printf(itXit F(Xi)t 1阶t 2阶tt3阶t 4阶n); pr

7、intf(-); for(i=0;iN+1;i+) r=i; printf(x%dt%f ,i,xi); for(j=0;j=i;j+) printf(%f ,trj); r-; printf(n); printf(-);/*/ printf(N(x) =n); printf( %fn,t00); for(i=1;iN+1;i+) for(j=0;ji;j+) mi=mi*(varx-xj); mi=t0i*mi; sum=sum+mi; for(i=1;iN+1;i+) printf( +%f*,t0i); for(j=0;ji;j+) printf(x-%f),xj); printf(n)

8、; return sum;void main() double varx,fN+1N+1; double xN+1=0.4,0.55,0.8,0.9,1; double yN+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(nn牛顿插值结果: P(%f)=%fn,varx,compvalue(f,x,varx); else printf(输入的插值节点的x值必须互异!); get

9、ch();运行及结果显示：实验三：自适应梯形公式命名（源程序.cpp及工作区.dsw）：autotrap问题：计算的近似值，使得误差（EPS）不超过/autotrap.cpp#include#include#include#define EPS 1e-6double f(double x) return 4/(1+x*x);double AutoTrap(double(*f)(double),double a,double b,double eps) double t2,t1,sum=0.0; int i,k; t1=(b-a)*(f(a)+f(b)/2; printf(T(%4d)=%fn,

10、1,t1); for(k=1;k11;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)=%fn,(int)pow(2,k),t2); if(fabs(t2-t1)EPS) break; else t1=t2; sum=0.0; return sum;void main() double s; s=AutoTrap(f,0.0,1.0,EPS); getch();运行及结果显示：实验四：龙贝格算法命名（源程序.cpp及工作区.dsw）：ro

11、mberg问题：求的近似值，要求误差不超过/romberg.cpp#include #include #include #define N 20#define EPS 1e-6double f(double x)return 4/(1+x*x);double Romberg(double a,double b,double (*f)(double),double eps) double TNN,p,h; int k=1,i,j,m=0; T00=(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=

12、p+f(a+(2*i-1)*h/2); T0k=T0k-1/2+p*h/2; for(i=1;i=EPS); k-; while(m=k) for(i=0;i=m;i+) printf(T(%d,%d)=%f ,i,m-i,Tim-i); m+; printf(n); return Tk0;void main() printf(nn积分结果 =%f,Romberg(0,1,f,EPS); getch();运行及结果显示：实验五：牛顿切线法问题：求方程在，附近的根（精度）/newton_qxf.cpp#include #include#include #define MaxK 1000#defi

13、ne EPS 0.5e-3double f(double x) return x*x*x-x-1;double f1(double x) return 3*x*x-1;int newton(double (*f)(double), double (*f1)(double), double &x0, double eps) double xk, xk1; int count = 0; printf(ktxkn); printf(-n); xk = x0; printf(%dt%fn, count, xk); do if(*f1)(xk)=0) return 2; xk1 = xk - (*f)(

14、xk) / (*f1)(xk); if (fabs(xk - xk1) eps) count+; xk = xk1; printf(%dt%fn, count, xk); x0 = xk1; return 1; count+; xk = xk1; printf(%dt%fn, count, xk); while(count MaxK); return 0;void main() for(int i=0;i2;i+) double x=0.6; if(i=1) x=1.5; printf(x0初值为%f:n,x); if (newton(f, f1, x, EPS) = 1) printf(-n

15、); printf(the root is x=%fnnn, x); else printf(the method is fail!); getch();实验六：牛顿下山法命名（源程序.cpp及工作区.dsw）：newton_downhill问题：求方程在，附近的根（精度）/newton_downhill.cpp#include #include #include #include #define Et 1e-3/下山因子下界#define E1 1e-3/根的误差限#define E2 1e-3/残量精度double f(double x) return x * x * x - x - 1;

16、double f1(double x) return 3 * x * x - 1;void errormess(int b) char *mess; switch(b) case -1: mess = f(x)的导数为0!; break; case -2: mess = 下山因子已越界，下山处理失败; break; default: mess = 其他类型错误!; printf(the method has faild! because %s, mess);int Newton(double (*f)(double), double (*f1)(double), double &x0) int

17、 k = 0; double t; double xk, xk1; double fxk, fxk_temp; printf(k t xk f(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) / (*

18、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 Et) return -2; else printf(下山成功n); break; if (fabs(xk-xk1)E1) x0 = xk1; return 1; xk = xk1; void main() int b; double x0; x0 = 0.6; b = Newton(f, f1, x0

19、); if (b = 1) printf(nthe root x=%fn, x0); else errormess(b); getch();运行及结果显示：实验七：埃特金加速算法命名（源程序.cpp及工作区.dsw）：aitken问题：求方程在，附近的根（精度）/aitken.cpp#include #include #include #define MaxK 100#define EPS 0.5e-3double g(double x) return x * x * x - 1;int aitken (double (*g)(double), double &x, double eps) i

20、nt k; double xk = x, yk, zk, xk1; printf(k xk yk zk xk+1n); printf(-n); for (k = 0;kMaxK; k+) yk = (*g)(xk); zk = (*g)(yk); xk1 = xk - (yk - xk) * (yk - xk) / (zk - 2 * yk + xk); printf(%-10d%-15f%-15f%-15f%-15fn, k, xk, yk, zk, xk1); if (fabs(yk-xk)=eps) x = xk1; return k + 1; xk = xk1; return -1;v

21、oid main () double x = 1.5; int k; k = aitken(g, x, EPS); if (k = -1) printf(迭代次数越界!n); else printf(n 经k=%d次迭代，所得方程根为:x=%fn, k, x); getch();运行及结果显示：实验八：正割法问题：求方程在附近的根（精度0.5e-8）/ZhengGe.cpp#include #include #include #define MaxK 1000#define EPS 0.5e-8double f(double x) return x*x*x-x-1;int ZhengGe(do

22、uble (*f)(double), double x0, double &x1, double eps) printf(k xk f(xk)n); printf(-n); int k; double xk0, xk, xk1; xk0 = x0; printf(%-10d%-15f%-15fn, 0, x0, (*f)(x0); xk = x1; for (k=1;kMaxK;k+) if(*f)(xk)-(*f)(xk0)=0) return -1; xk1 = xk - (*f)(xk) / ( (*f)(xk) - (*f)(xk0) ) * (xk - xk0); printf(%-

23、10d%-15f%-15fn, k, xk, (*f)(xk); if (fabs(xk1 - xk)=eps) printf(%-10d%-15f%-15fnn, k + 1, xk1, (*f)(xk1); printf(-nn); x1 = xk1; return 1; xk0 = xk; xk = xk1; return -2;void main () double x0 = 1, x1 = 2; if (ZhengGe(f, x0, x1, EPS) = 1) printf(the root is x = %fn, x1); else printf(the method is fail!); getch();实验九：高斯列选主元算法命名（源程序.cpp及工作区.dsw）：colpivot问题：求解方程组并计算系数矩阵行列式值/colpivot.cpp#include #in