《数值计算方法实验程序(0)jsp》由会员分享，可在线阅读，更多相关《数值计算方法实验程序(0)jsp（26页珍藏版）》请在金锄头文库上搜索。
1、1第一部分 C 语言程序部分一、误差与数值稳定性一、误差与数值稳定性例 1 数值方法计算积分收敛和发散实例dxexeIxn n101法 1In=1-nIn-1 #include #include main() int n,j; float i; i=1-exp(-1); printf(“%2d %-20.8f“,0,i); for(j=1;j #include #define N 20 main() int j; float I; I=0;printf(“%2d %-20.8f“,20,I); for(j=19;j=1;j-)I=1.0/j*(1-I);printf(“%2d %-20.8f“
2、,j,I);if(j%3=0) printf(“n“); printf(“n“); 二、大数吞掉小数的例子二、大数吞掉小数的例子#include“conio.h“ main() double x=7e+33,y=0.001; clrscr();2printf(“1=%dn“,x #include #include #define EPS 1e-2 #define DELTA 1e-6 float Bisection(float,float,float(*f)(float); float f(float); main() float a=1.5,b=2; float x; clrscr(); x
3、=Bisection(a,b,f); printf(“nThe root of the equation is%f“,x); getch(); float Bisection(float a,float b,float(*f)(float)float c,fc,fa=f(a),fb=f(b);int n=1;printf(“total bisection ctt f(c)n“);printf(“n a b c fa fb fcn“);while(1) if (fa*fb0) exit();c=(a+b)/2,fc=f(c);if(fabs(fc) #include float diedai(f
4、loat); main() int i;float x0,x1,e;scanf(“%f %f“,clrscr();for (i=1;i+) x1=diedai(x0); if(fabs(x1-x0) #include #include #define N 100 #define EPS 1e-6 #define ETA 1e-8 main() float f(float); float f1(float); float x0,y0; float newton(float(*)(float),float(*)(float),float); clrscr(); printf(“please inp
5、ut x0n“); scanf(“%f“, printf(“x(0)=%fn“,x0);4y0=newton(f,f1,x0); printf(“nthe root of the equation is x=%fn“,y0); getch(); float newton(float (*f)(float),float(*f1)(float),float x0) float x1,d; int k=0; dox1=x0-f(x0)/f1(x0);if (k+N|(fabs(f1(x1)EPS return x1; float f(float x) return x-cos(x); float f
6、1(float x) return 1+sin(x); 3弦割法弦割法例：用弦割法求方程 f(x)=x3-x-1=0 的根 #include #include #include #define EPS 1e-6 main()float f(float);float x0,x1,y0;float xuan(float(*)(float),float,float);clrscr();printf(“please input x0,x1:n“);scanf(“%f %f“,printf(“x(0)=%f,x(1)=%fn“,x0,x1);xuan(f,x0,x1);printf(“nthe root
7、 of the equation is x=%fn“,x1);getch();5 float xuan(float (*f)(float),float x0,float x1)float x2;int k=0;do x2=x1-f(x1)/(f(x1)-f(x0)*(x1-x0);if(fabs(x1-x0)(fabs(*(c+k*(n+1)+i) k=j;if(k!=i)for(j=i;j=0;i-)for(j=n-1;j=i+1;j-)(*(c+i*(n+1)+n)-=xj*(*(c+i*(n+1)+j);xi=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i);return x
8、; 第二部分第二部分 EXCEL 实验部分实验部分一、一、 积分积分收敛和发散实例收敛和发散实例dxexeIxn n101具体操作说明： 1） 第 A 列，第一、二行文字是直接输入的 2） 在 B3 中输入=1-exp(-1)即为 I0 3） 在 B4 中输入=1-A4*B3 即为 I1=1-1*I0 4） 选中 B4 当鼠标成黑十字星时向下拖拽即可计算出 In7二、二、 多项式计算多项式计算1 1直接法直接法 例求多项式 Pn(x)=2+3x+5x2+6x3-3x4-x5+4x6+8x7在 x=1.4、1.5 的值 f(1.4)、f(1.5) f(1.4)=130.00991，具体算法说明：
9、 1)数据的输入：第一行中 A1 中的 1.4 及 0，1，2，3，7 直接输入的分别表示 x=1.4 和 i，第 三行中的元素是多项式的系数也是直接输入的，在 b2 中输入“=1”,b4 中输入“=B2*B3” ，即就是 a0*x0，在 J4 中输入“=sum(B4:I4)即就是。 70ii ixa2)在 C2 中输入=$A$1*B2，选定 C2 将光标放在右下角，当鼠标变成黑十字星后，向右拖拽到 I7，释放鼠标，Excel 将进行自动填充，这也叫关系复制，其实质是计算 xi-1*x。 3)选定 B4 操作如 2）其实质是计算 ai*xi。 4)在 A1 中输入 1.5 回车即可实现 f(1
10、.5)=197.46875 的计算。82 2秦九韶算法秦九韶算法 同样以 1 的多项式 Pn(x)=2+3x+5x2+6x3-3x4-x5+4x6+8x7为例，求 f(x)在 x=1.4、1.5 的值 f(1.4)、 f(1.5)， 操作说明： 1)A、B、C 列中的数据是直接输入的 2)在 D2 输入“=B2”，即就是 u0=a7，在 D3 中输入“=D2*$B$1+B3” ，即就是 u1=u0*x+a6 3)选 定 D3 向下进行关系复制拖拽到 D9 为止 （这里实质是 ui=uk-1*x+an-k,k=1,2, 7，即得 u7=130.00991，即有 f(1.4)= u7三、二分法三、
11、二分法例：例：例用二分法求方程4sin)(2xxxf的非零实根的近似值，使误差不超过 10-2具体操作如下： 1 第一行的各单元格中的符号是直接输入的，第一列的序号也是直接输入的代表二分次数 2B2，C2 中的 1.5，2 直接输入的，在 D2 中输入=（b2+c2）/2,在 E2 中输入=sin(b2)-b22/4,选 中 E2 鼠标变成黑十字星后拖拽到 F2，G2 即得到相当于 f(b1),f(x1)的值， 3.在 H2 中输入=if(abs(c2-b2)/21 时，同时用到 f(xk)1 时的误差控制方法，当|xk+1-xk|emgfor i=1:nsum=0;for j=1:nif i
12、=jsum=sum+A(i,j)*x1(j);endendx2(i)=(b(i)-sum)/A(i,i);endr=max(abs(x2-x1);x1=x2;k=k+1;if kNdisp(迭代失败,返回);return;end end x=x1; 例 8：用迭代法解线性方程组 52532121 xxxxA=3 1;1 2 b=5 5 N=10 x0=0 0 emg=0.001 Jacobimethod(A,b,N,x0,emg)2Gauss_Seidel 迭代法迭代法function x,k=Gaussmethod(A,b,N,x0,emg)19%A 是线性方程组的系数矩阵 %B 是线性方程
13、组的右端向量 %x0 是迭代初始值 %N 是表示迭代次数上限若迭代次数,则迭代失败 %emg 表示控制精度 %k 表示迭代次数 %x 表示线性方程组的近似解 n=length(b); x1=zeros(n,1); x2=zeros(n,1); x1=x0; k=0; r=max(abs(b-A*x1); while remgfor i=1:nsum=0;for j=1:nif jisum=sum+A(i,j)*x1(j);elseif jNdisp(迭代失败,返回);return;end end x=x1;3超松驰迭代法超松驰迭代法例 11 用 SOR 法解线性方程组 731091041017
14、2424321xxxfunction x,k=SORmethod(A,b,N,x0,emg,w) %A 是线性方程组的系数矩阵 %B 是线性方程组的右端向量 %x0 是迭代初始值 %N 是表示迭代次数上限若迭代次数,则迭代失败20%emg 表示控制精度 %k 表示迭代次数 %x 表示线性方程组的近似解 n=length(b); x1=zeros(n,1); x2=zeros(n,1); x1=x0; k=0; r=max(abs(b-A*x1); while remgfor i=1:nsum=0;for j=1:nif jisum=sum+A(i,j)*x1(j);elseif jNdisp(
15、迭代失败,返回);return;end end x=x1; A=4 -2 -4;-2 17 10;-4 10 9 b=10 3 -7 N=25 x0=0 0 0 w=1.46 SORmethod(A,b,N,x0,emg,w)四、四、Lagrange 插值插值例 给定函数表如下：试用线性插值与抛物插值求 e0.285 的近似值，并估计截断误差。function y0,N =Laugrange_eval(X,Y,x0) %X,Y 是已知插值点坐标21%x0 是插值点 %y0 是 Laugrange 插值多项式在 xo 点处的值 %N 是 laugrange 插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:mN(i)=1;for j=1:mif j=iN(i)=N(i)*(x0-X(j)/(X(i)-X(j);endendy0=y0+Y(i)*N(i); end1.求插值值求插值值的程序的程序0.285eX=0.2 0.3 Y=1.2214 1.3499 x0=0.285 y0,N=Laugrange(X,Y,x0) X=0.2 0.3 0.4 Y=1.2214 1.3499 1.492 x0=0.285 y0,N=Laugrange(X,
