数值计算方法编程作业C语言版讲解

1、1:第二章(1)二分法求解非线性方程： #include#include #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)*

2、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,-5The root is x=1.324721, k=20Press any key to continue 总结：本题关键在于两个端点的取值和误差的判断，此程序较容易。二分法收敛速度较快， 但缺点是只能求解单根。(2)

3、牛顿法求解非线性方程：#include #include 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)e

4、ps); /* 当满足精度 ,输出近似根 */ printf(x=%fn,x1);程序运行结果：x=1.265328总结：关键是牛顿迭代的应用， 程序中最大缺点是函数及其导数已唯一给出确定不可求的随 意函数的根，牛顿法比二分法收敛快，可以求重根。2:第三章(1)列主元素消去法求解线性方程： #include#include#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;coutm;coutendl;cout 请按顺序输入增广矩

5、阵 a:endl;load();for(i=0;im;i+)for(j=i;jfabs(aii)?j:i;/* 找列最大元素 */for(n=0;nm+1;n+)s=ain; ain=acn; acn=s;/* 将列最大数防在对角线上 */for(p=0;pm+1;p+)coutaipt;coutendl;for(k=i+1;km;k+)lki=aki/aii;for(r=i;r=0;i-)d=0;for(j=i+1;jm;j+)d=d+aij*xj; xi=(aim-d)/aii;cout 该方程组的解为 :endl; for(i=0;im;i+) coutxi=xit; /system(p

6、ause); return 0;void load()int i,j;for(i=0;im;i+)for(j=0;jaij;运行结果：/* 化成三角阵 */*求解 */1 2 3 45 1 0 84 6 9 246920-6.5-11.255.50-1.86265e-008-0.115385该方程组的解为请按顺序输入增广矩阵 a:3.92308下面请输入未知数的个数 m=3 x0=-9.99999 x1=58 x2=-34 Press any key to continue 总结：列主元素消去法的目的是为了防止减去一个较小的数时大数淹没小数， 而使结果产生 较大误差， 本程序关键在每次消元时找

7、到相应列中的最大项， 然后交换两行位置， 在进行计 算。(2)LU 分解法求解线性方程：#includevoid 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=1;i-) /* 第二次回代过程开始 */ s2=0;for(j=n;ji;j-)t=-uij; s2=s2+t*xj; xi=(yi+s2)/uii;void main()float a100100,l100100,u100

8、100,x100,b100;int i,j,n,r,k;float s1,s2;1*/for(i=1;i=99;i+)/* 将所有的数组置零，同时将 L 矩阵的对角值设为 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(

9、i=1;i=n;i+)scanf(%f,&bi);for(r=1;r=n;r+)/* 求解矩阵 L 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);/* 输出矩阵for(i=1;i=n;i+)for(j=1;j=n;j+)printf(%7.3f ,lij);printf(n);printf(array U:n);/* 输出矩阵for(i=1;

10、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 5和 U*/L*/U*/input array B:3 1 -7array L:1.0000.0000.0002.0001.0000.000-1.0002.0001.000array U:2.0002.0003.0000.0003.0001.0000.0000.0006.

11、000解为 :x1=2.000000 x2=-2.000000x3=1.000000Press any key to continue总结：关键是把矩阵分解为 L、 U 两个三角矩阵，回代过程比较简单。3:第四章(1) 拉格朗日差值多项式 ; #include#include#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

12、); scanf(%d,&N);printf(n); for(i=0;iN;i+)scanf(%f,%f,&xi,&yi); for(i=0;i=m;i+)for(j=i;j=m;j+)tmp=0; for(k=0;kN;k+) tmp=tmp+pow(xk,(i+j);Xij=tmp;Xji=Xij;for(i=0;i=m;i+)tmp=0;for(k=0;kN;k+)tmp=tmp+yk*pow(xk,i);Yi=tmp;for(j=0;jm;j+)for(i=j+1,mi=j,mx=fabs(Xjj);imx)mi=i;mx=fabs(Xij);if(jmi)tmp=Yj;Yj=Ymi;

13、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=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.9804

14、17+(0.282759)*x1+(-0.299937)*x2+(0.022071)*x3+(0.016624)*x4+(-0.0 01934)*x5Press any key to continue 总结：拉格朗日计算公式中，只需要知道各个点即可4:第五章(1)曲线拟合：#include#include#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 输入给定点的个数

15、 n 及坐标 (x,y):n);scanf(%d,&N);printf(n);for(i=0;iN;i+)scanf(%f,%f,&xi,&yi);for(i=0;i=m;i+) for(j=i;j=m;j+)tmp=0; for(k=0;kN;k+) tmp=tmp+pow(xk,(i+j); Xij=tmp;Xji=Xij;for(i=0;i=m;i+)tmp=0;for(k=0;kN;k+)tmp=tmp+yk*pow(xk,i);Yi=tmp;for(j=0;jm;j+)for(i=j+1,mi=j,mx=fabs(Xjj);imx)mi=i;mx=fabs(Xij);if(jmi)t

16、mp=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=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所求的二次多项

17、式为 :P(x)=3.952280+(-0.506315)*x1+(0.050877)*x2Press any key to continue5:第六章(1)辛普生求积方法：#include #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=

18、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,3 s=1.854590Press any key to continue总结：辛普生算法是一种积分方法，采用三点法插

19、值，如果 h 较小的话，误差很小，因为它 b a h 4的插值余项 R( f ) b a h f (4)( ) ，辛普生算法比较精确，程序关键是对所取的点 180 2的取和，注意6:第七章(1)改进欧拉法求解常微分方程的初值问题#include 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+

20、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:20 input x0,xn:1,6 input y0:2 y(2.000000)=7.0000Press any

21、key to continue (2)四阶龙格 库塔法 #include 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