数值计算方法编程作业(C语言知识版).doc
1:第二章(1)二分法求解非线性方程：#include<stdio.h>#include <math.h>#define f(x) (x*x-1)*x-1) void main() float a,b,x,eps; int k=0; printf("intput epsn");/* 容许误差 */ scanf("%f",&eps); printf("a,b=n"); for(;) scanf("%f, %f",&a ,&b); if(f(a)*f(b)>=0) /* 判断是否符合二分法使用的条件 */ printf("二分法不可使用,请重新输入:n"); else break; do x=(a+b)/2; k+; if(f(a)*f(x)<0) /* 如果f(a)*f(x)<0，则根在区间的左半部分 */ b=x; else if(f(a)*f(x)>0) /* 否则根在区间的右半部分 */ a=x; else break; while(fabs(b-a)>eps);/*判断是否达到精度要求,若没有达到,继续循环*/ x=(a+b)/2; /* 取最后的小区间中点作为根的近似值 */ printf("n The root is x=%f, k=%dn",x,k);运行结果：intput eps0.00001a,b=2,-5 The root is x=1.324721, k=20Press any key to continue总结：本题关键在于两个端点的取值和误差的判断，此程序较容易。二分法收敛速度较快，但缺点是只能求解单根。 (2)牛顿法求解非线性方程：#include <stdio.h>#include <math.h>float f(float x) /* 定义函数f(x) */ return(-3*x+4)*x-5)*x+6; float f1(float x) /* 定义函数f(x)的导数 */ return (-9*x+8)*x-5; void main() float eps,x0,x1=1.0; printf("input eps:n"); scanf("%f",&eps); /* 输入容许误差 */ do x0=x1; /* 准备下一次迭代的初值 */ x1=x0-f(x0)/f1(x0); /* 牛顿迭代 */ while(fabs(x1-x0)>eps); /*当满足精度,输出近似根*/ printf("x=%fn",x1);程序运行结果：x=1.265328总结：关键是牛顿迭代的应用，程序中最大缺点是函数及其导数已唯一给出确定不可求的随意函数的根，牛顿法比二分法收敛快，可以求重根。2:第三章(1)列主元素消去法求解线性方程：#include<iostream>#include<cmath>#define N 20using namespace std;void load();float aNN;int m;int main()int i,j;int c,k,n,p,r;float xN,lNN,s,d;cout<<"下面请输入未知数的个数m="cin>>m;cout<<endl;cout<<"请按顺序输入增广矩阵a:"<<endl;load();for(i=0;i<m;i+) for(j=i;j<m;j+) c=(fabs(aji)>fabs(aii)?j:i; /*找列最大元素*/for(n=0;n<m+1;n+) s=ain; ain=acn; acn=s; /*将列最大数防在对角线上*/for(p=0;p<m+1;p+)cout<<aip<<"t"cout<<endl;for(k=i+1;k<m;k+) lki=aki/aii; for(r=i;r<m+1;r+) /*化成三角阵*/ akr=akr-lki*air; xm-1=am-1m/am-1m-1;for(i=m-2;i>=0;i-) d=0;for(j=i+1;j<m;j+)d=d+aij*xj;xi=(aim-d)/aii; /*求解*/cout<<"该方程组的解为:"<<endl;for(i=0;i<m;i+)cout<<"x"<<i<<"="<<xi<<"t" /system("pause");return 0;void load()int i,j;for(i=0;i<m;i+)for(j=0;j<m+1;j+)cin>>aij;运行结果：下面请输入未知数的个数m=3请按顺序输入增广矩阵a:1 2 3 45 1 0 84 6 9 24 6 9 20 -6.5 -11.25 5.50 -1.86265e-008 -0.115385 3.92308该方程组的解为:x0=-9.99999 x1=58 x2=-34 Press any key to continue总结：列主元素消去法的目的是为了防止减去一个较小的数时大数淹没小数，而使结果产生较大误差，本程序关键在每次消元时找到相应列中的最大项，然后交换两行位置，在进行计算。(2)LU分解法求解线性方程：#include<stdio.h>void solve(float l100,float u100,float b,float x,int n)int i,j;float t,s1,s2;float y100;for(i=1;i<=n;i+) /* 第一次回代过程开始 */ s1=0; for(j=1;j<i;j+) t=-lij; s1=s1+t*yj; yi=(bi+s1)/lii; for(i=n;i>=1;i-) /* 第二次回代过程开始 */ s2=0; for(j=n;j>i;j-) t=-uij; s2=s2+t*xj; xi=(yi+s2)/uii; void main()float a100100,l100100,u100100,x100,b100;int i,j,n,r,k;float s1,s2;for(i=1;i<=99;i+)/*将所有的数组置零，同时将L矩阵的对角值设为1*/ for(j=1;j<=99;j+) lij=0,uij=0; if(j=i) lij=1; printf ("input n:n");/*输入方程组的个数*/scanf("%d",&n);printf ("input array A:n");/*读取原矩阵A*/for(i=1;i<=n;i+) for(j=1;j<=n;j+) scanf("%f",&aij);printf ("input array B:n");/*读取列矩阵B*/for(i=1;i<=n;i+) scanf("%f",&bi);for(r=1;r<=n;r+)/*求解矩阵L和U*/for(i=r;i<=n;i+) s1=0; for(k=1;k<=r-1;k+)s1=s1+lrk*uki;uri=ari-s1; for(i=r+1;i<=n;i+)s2=0; for(k=1;k<=r-1;k+)s2=s2+lik*ukr;lir=(air-s2)/urr;printf("array L:n");/*输出矩阵L*/for(i=1;i<=n;i+) for(j=1;j<=n;j+) printf("%7.3f ",lij); printf("n");printf("array U:n");/*输出矩阵U*/for(i=1;i<=n;i+) for(j=1;j<=n;j+) printf("%7.3f ",uij); printf("n");solve(l,u,b,x,n);printf("解为:n");for(i=1;i<=n;i+)printf("x%d=%fn",i,xi);运行结果：input n:3input array A:2 2 34 7 7-2 4 5input array B:3 1 -7array L: 1.000 0.000 0.000 2.000 1.000 0.000 -1.000 2.000 1.000array U: 2.000 2.000 3.000 0.000 3.000 1.000 0.000 0.000 6.000解为:x1=2.000000x2=-2.000000x3=1.000000Press any key to continue总结：关键是把矩阵分解为L、U两个三角矩阵，回代过程比较简单。3:第四章(1)拉格朗日差值多项式;#include<stdio.h>#include<math.h>#define MAX 100void main() int i,j,k,m,n,N,mi; float tmp,mx; float XMAXMAX,YMAX,xMAX,yMAX,aMAX; printf("n 输入拟合多项式的次数:n"); scanf("%d",&m); printf("n 输入给定点的个数n及坐标(x,y):n"); scanf("%d",&N); printf("n"); for(i=0;i<N;i+) scanf("%f,%f",&xi,&yi); for(i=0;i<=m;i+) for(j=i;j<=m;j+) tmp=0; for(k=0;k<N;k+) tmp=tmp+pow(xk,(i+j); Xij=tmp; Xji=Xij; for(i=0;i<=m;i+) tmp=0; for(k=0;k<N;k+) tmp=tmp+yk*pow(xk,i); Yi=tmp; for(j=0;j<m;j+) for(i=j+1,mi=j,mx=fabs(Xjj);i<=m;i+) if(fabs(Xij)>mx) mi=i; mx=fabs(Xij); if(j<mi) tmp=Yj;Yj=Ymi;Ymi=tmp; for(k=j;k<=m;k+) tmp=Xjk; Xjk=Xmik; Xmik=tmp; for(i=j+1;i<=m;i+) tmp=-Xij/Xjj; Yi+=Yj*tmp; for(k=j;k<=m;k+) Xik+=Xjk*tmp; am=Ym/Xmm; for(i=m-1;i>=0;i-) ai=Yi;for(j=i+1;j<=m;j+) ai-=Xij*aj; ai/=Xii;printf("n 所求的二次多项式为:n");printf("P(x)=%f",a0);for(i=1;i<=m;i+)printf("+(%f)*x%d",ai,i);运行结果： 输入拟合多项式的次数:5 输入给定点的个数n及坐标(x,y):31,25,34,2 所求的二次多项式为:P(x)=1.980417+(0.282759)*x1+(-0.299937)*x2+(0.022071)*x3+(0.016624)*x4+(-0.001934)*x5Press any key to continue总结：拉格朗日计算公式中，只需要知道各个点即可4:第五章(1)曲线拟合：#include<stdio.h>#include<math.h>#define MAX 100void main() int i,j,k,m,n,N,mi; float tmp,mx; float XMAXMAX,YMAX,xMAX,yMAX,aMAX; printf("n 输入拟合多项式的次数:n"); scanf("%d",&m); printf("n 输入给定点的个数n及坐标(x,y):n"); scanf("%d",&N); printf("n"); for(i=0;i<N;i+) scanf("%f,%f",&xi,&yi); for(i=0;i<=m;i+) for(j=i;j<=m;j+) tmp=0; for(k=0;k<N;k+) tmp=tmp+pow(xk,(i+j); Xij=tmp; Xji=Xij; for(i=0;i<=m;i+) tmp=0; for(k=0;k<N;k+) tmp=tmp+yk*pow(xk,i); Yi=tmp; for(j=0;j<m;j+) for(i=j+1,mi=j,mx=fabs(Xjj);i<=m;i+) if(fabs(Xij)>mx) mi=i; mx=fabs(Xij); if(j<mi) tmp=Yj;Yj=Ymi;Ymi=tmp; for(k=j;k<=m;k+) tmp=Xjk; Xjk=Xmik; Xmik=tmp; for(i=j+1;i<=m;i+) tmp=-Xij/Xjj; Yi+=Yj*tmp; for(k=j;k<=m;k+) Xik+=Xjk*tmp; am=Ym/Xmm; for(i=m-1;i>=0;i-) ai=Yi;for(j=i+1;j<=m;j+) ai-=Xij*aj; ai/=Xii;printf("n 所求的二次多项式为:n");printf("P(x)=%f",a0);for(i=1;i<=m;i+)printf("+(%f)*x%d",ai,i); 输入拟合多项式的次数:2 输入给定点的个数n及坐标(x,y):51,25,32,48,3-1,5 所求的二次多项式为:P(x)=3.952280+(-0.506315)*x1+(0.050877)*x2Press any key to continue5:第六章(1)辛普生求积方法：#include <stdio.h>#define N 16 /* 等分数 */float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float h) int i; for(i=0;i<=N;i+) yi=func(a+i*h);float simpson(float y,float h) float s,s1,s2; int i; s1=y1; s2=0.0; for(i=2;i<=N-2;i=i+2) s1+=yi+1; /* 计算奇数项的函数值之和 */ s2+=yi; /* 计算偶数项的函数值之和 */ s=y0+yN+4.0*s1+2.0*s2; return(s*h/3.0);main() float a,b,h,s,fN+1; scanf("%f,%f",&a,&b); h=(b-a)/( float)N; gedianzhi(f,a,h); s=simpson(f,h); printf("s=%fn",s);运行结果：1,3s=1.854590Press any key to continue总结：辛普生算法是一种积分方法，采用三点法插值，如果h较小的话，误差很小，因为它的插值余项，辛普生算法比较精确，程序关键是对所取的点的取和，注意6:第七章(1)改进欧拉法求解常微分方程的初值问题#include <stdio.h>float func(float x,float y) return(y-x);float euler(float x0,float xn,float y0,int N) float x,y,yp,yc,h; int i; x=x0; y=y0; h=(xn-x0)/(float)N; for(i=1;i<=N;i+) yp=y+h*func(x,y); x=x0+i*h; yc=y+h*func(x,yp); y=(yp+yc)/2.0; return(y);main() float x0,xn,y0,e; int n; printf("ninput n:n "); scanf("%d",&n); printf("input x0,xn:n "); scanf("%f,%f",&x0,&xn); printf("input y0:n "); scanf("%f",&y0); e=euler(x0,xn,y0,n); printf("y(%f)=%6.4f",y0,e);input n: 20input x0,xn: 1,6input y0: 2y(2.000000)=7.0000Press any key to continue(2)四阶龙格库塔法#include <stdio.h>float func(float x,float y) return(x-y);float runge_kutta(float x0,float xn,float y0,int N) float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/(float)N; for(i=1;i<=N;i+) xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; return(y); main() float x0,xn,y0,e; int N; printf("ninput n:n "); scanf("%d",&N); printf("input x0,xn:n "); scanf("%f,%f",&x0,&xn); printf("input y0:n "); scanf("%f",&y0); e=runge_kutta(x0,xn,y0,N); printf("y(%f)=%8.6f",y0,e);input n: 10input x0,xn: 1,2input y0: 5y(5.000000)=2.833863Press any key to continue2-2 Gauss-Seidel方法#include <math.h>#include<stdio.h>int gsdl(a,b,n,x,eps)int n;double a,b,x,eps;int i,j,u,v;double p,t,s,q;for(i=0;i<=n-1;i+)u=i*n+i;p=0.0;xi=0.0;for(j=0;j<=n-1;j+)v=i*n+j;p=p+fabs(av);if(p>=fabs(au)printf(“failn”);return(-1);p=eps+1.0;while(p>=eps)for(i=0;i<=n-1;i+)t=xi;s=0.0;for(j=0;j<=n-1;j+)if(j!=i)s=s+ai*n+j*xj;xi=(bi-s)/ai*n+j;q=fabs(xi-t)/(1.0+fabs(xi);if(q>p)p=q;return(1);main()int i;double eps;static double a44=7,2,1,-29,15,3,-2-2,-2,11,51,3,2,13;static double x5,b4=4,7,-1,0;eps=0.000001;if(dsdl(a,b,4,x,eps)>0)for(i=0;i<=3;i+)printf(“x(%d)=%13.7en”,i,xi);
