数值分析上机实验指导书

《数值分析》上机指导书 谭立云 丁智斌 华北科技学院基础部信息与计算科学专业教研室 实验一：舍入误差与数值稳定性 1.目的与要求： 1） 通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； 2） 通过上机计算，了解舍入误差所引起的数值不稳定性。 2.内容概要： 舍入误差在计算方法中是一个很重要的概念。在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。因此，选取稳定的算法，在实际计算中是十分重要的。 3.程序与实例： 例1 对n=0,1,2,…20计算定积分 算法1 利用递推公式 n=1,2,…,20 取 算法2 利用递推公式 n=20,19,.…,1 注意到 取 算法1的程序和输出结果如下： /*数值不稳定算法*/ #include "stdafx.h" #include #include main() { float y_0=log(6.0)-log(5.0),y_1; int n=1; printf("y[0]=%-20f",y_0); while(1) { y_1=1.0/n-5*y_0; printf("y[%d]=%-20f",n,y_1); if(n>=20)break; y_0=y_1; n++; if(n%3==0)printf("\n"); } } 运行结果： y[0]=0.182322 y[1]=0.088392 y[2]=0.058039 y[3]=0.043138 y[4]=0.034310 y[5]=0.028448 y[6]=0.024428 y[7]=0.020719 y[8]=0.021407 y[9]=0.004076 y[10]=0.079618 y[11]=-0.307181 y[12]=1.619237 y[13]=-8.019263 y[14]=40.167745 y[15]=-200.772052 y[16]=1003.922745 y[17]=-5019.554824 y[18]=25097.828993 y[19]=-125489.087993 y[20]=627445.479688 算法2的程序和输出结果如下： /*稳定算法*/ #include "stdafx.h" #include #include main() { float y_0=(1/105.0+1/126.0)/2,y_1; int n=20; printf("y[20]=%-20f",y_0); while(1) { y_1=1/(5.0*n)-y_0/5.0; printf("y[%d]=%-20f",n-1,y_1); if(n<=1)break; y_0=y_1; n--; if(n%3==0)printf("\n"); } } 运行结果： y[20]=0.008730 y[19]=0.008254 y[18]=0.008876 y[17]=0.009336 y[16]=0.009898 y[15]=0.010520 y[14]=0.011229 y[13]=0.012040 y[12]=0.012977 y[11]=0.014071 y[10]=0.015368 y[9]=0.016926 y[8]=0.018837 y[7]=0.021233 y[6]=0.024325 y[5]=0.028468 y[4]=0.034306 y[3]=0.043139 y[2]=0.058039 y[1]=0.088392 y[0]=0.182322

说明

关于失联党员情况说明岗位说明总经理岗位说明书会计岗位说明书行政主管岗位说明书

：从计算结果可以看出，算法1是不稳定的，而算法2是稳定的。 上机实习题 1 1.用2种不同的顺序计算 ，分析其误差的变化。 2.给出一个有效的算法和一个无效的算法计算积分 ，n=0,1,…,10 3.设 ，已知其精确值为 。 1）编制按从大到小的顺序计算 的程序； 2）编制按从小到大的顺序计算 的程序； 3）按2种顺序分别计算 ， ， ，并指出有效位数。 实验二：插值法 1.目的与要求： 1）熟悉拉格朗日插值多项式和牛顿插值多项式，注意其不同特点； 2）会用插值解决一些实际问题。 2.拉格朗日插值多项式 算法 a) 输入 b) 对 计算 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 程序与实例 例1 已知函数表 Xi 0.56160 0.56280 0.56401 0.56521 Yi 0.82741 0.82659 0.82577 0.82495 用三次拉格朗日插值多项式求x=0.5635时的函数近似值。 #include "stdafx.h" #include"stdio.h" #include"conio.h" #include"malloc.h" float Lagrange(float *x,float *y,float xx,int n) { int i,j; float *a,yy=0.0; a=(float *) malloc(n*sizeof(float)); for(i=0;i<=n-1;i++) { a[i]=y[i]; for(j=0;j<=n-1;j++) if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]); yy+=a[i]; } free(a); return yy; } void main() { float x[4]={0.56160,0.56280,0.56401,0.56521}; float y[4]={0.82741,0.82659,0.82577,0.8245}; float xx=0.5635,yy; float Lagrange(float *,float *,float ,int); yy=Lagrange(x,y,xx,4); printf("x=%f,y=%f\n",xx,yy); } 运行结果： x=0.563500, y=0.826116 3.牛顿插值多项式 算法 eq \o\ac(○,1)输入 ； eq \o\ac(○,2)对 计算各阶差商 ； 计算函数值 程序与实例 例2 已知函数表 Xi 0.4 0.55 0.65 0.8 0.9 Yi 0.41075 0.57815 0.69675 0.88811 1.02652 用牛顿插值多项式求 （0.596）和 （0.895）。 程序： #include "stdafx.h" #include"stdio.h" #include"conio.h" #include"malloc.h" #define N 4 void Difference(float *x,float *y,int n) { float *f; int k,i; f=(float *)malloc(n*sizeof(float)); for(k=1;k<=n;k++) { f[0]=y[k]; for(i=0;i=0;i--) b=b*(varx-x[i])+y[i]; printf("N(%f)=%f\n",varx,b); } 运行结果： N(0.596000)=0.631918 （注意：要求其他值，就只要varx=0.895，就可得 N (0.895)=1.019368） 上机实习题2 1.按下列数据 EMBED Equation.3 —3.0 —1.0 1.0 2.0 3.0 1.0 1.5 2.0 2.0 1.0 作二次插值，并求 时的函数近似值。 2.按下列数据 EMBED Equation.3 0.30 0.42 0.50 0.58 0.66 0.72 1.04403 1.08462 1.11803 1.15603 1.19817 1.23223 作五次插值，并求 时的函数近似值。 3.编写一个用牛顿前插值公式计算函数值的程序，并求先输出差分表，再计算x点的函数值，并应用于下面的问题： EMBED Equation.3 20 21 22 23 24 1.30103 1.32222 1.34242 1.36173 1.38021 求x=21.4时的三次插值多项式的值。 实验三 曲线拟合 1.目的与要求： 1) 了解最小二乘法的基本原理，通过计算机解决实际问题； 2) 了解超定方程组的最小二乘法。 2. 最小二乘法 算法 已知数据对 。求多项式 ，使得 为最小。注意到此时 多项式系数 满足下面的线性方程组： 其中 ， 然后只要调用解线性方程组的函数程序即可。 程序与实例 例1 由化学实验得到某物质浓度与时间的关系如下： 时间t 1 2 3 4 5 6 7 8 浓度y 4.00 6.40 8.00 8.80 9.22 9.50 9.70 9.86 时间t 9 10 11 12 13 14 15 16 浓度y 10.00 10.20 10.32 10.42 10.50 10.55 10.58 10.60 求浓度与时间的二次拟合曲线。 程序： // lyf618-11.cpp : Defines the entry point for the console application. #include "stdafx.h" #include"stdio.h" #include"malloc.h" #include"math.h" void main() { int i; float *a; float x[16]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16}; float y[16]={4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86,10.00,10.20,10.32,10.42,10.50,10.55,10.58,10.60}; float *Approx(float *,float *,int,int); a=Approx(x,y,16,2); for(i=0;i<=2;i++) printf("a[%d]=%f\n",i,a[i]); } float *Approx(float *x,float *y,int m,int n) { float *c,*a; int i,j,t; float power(int ,float); float *ColPivot(float *,int); c=(float *)malloc((n+1)*(n+2)*sizeof(float)); for(i=0;i<=n;i++) { for(j=0;j<=n;j++) { *(c+i*(n+2)+j)=0.0; for(t=0;t<=m-1;t++) *(c+i*(n+2)+j)+=power(i+j,x[t]); } *(c+i*(n+2)+n+1)=0.0; for(j=0;j<=m-1;j++) *(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]); } a=ColPivot((float *) c,n+1); return a; } float *ColPiovt(float *a,int n) { int i,j,t,k; float *x,*c,p; x=(float *)malloc(n*sizeof(float)); c=(float *)malloc(n*(n+1)*sizeof(float)); for(i=0;i<=n-1;i++) for(j=0;j<=n;j++) *(c+i*(n+1)+j)=*(a+i*(n+1)+j); for(i=0;i<=n-2;i++) { k=i; for(j=i+1;j<=n-1;j++) if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i)))) k=j; if(k!=i) for(j=i;j<=n;j++) { p=*(c+i*(n+1)+j); *(c+i*(n+1)+j)=*(c+k*(n+1)+j); *(c+k*(n+1)+j)=p; } for(j=i+1;j<=n-1;j++) { p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i)); for(t=i;t<=n-1;t++) *(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t)); *(c+j*(n+1)+n)-=*(c+i*(n+1)+n)*p; } } for(i=n-1;i>=0;i--) { for(j=n-1;j>=i+1;j--) (*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j)); x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i)); } free(c); return x; } float power(int i,float v) { float a=1.0; while(i--) a*=v; return a; } 运行结果： a[0]=4.387500 a[1]=1.065962 a[2]=-0.044466 因此二次拟合多项式为y(t)=4.3875+1.065962t-0.044466 例2 工程中经常遇到y= 类型的曲线，其中a与b是正数。试求最小二乘法拟合下面的数据： 5 10 15 20 23.75 22.50 35.75 42.50 44.00 42.00 令z= ,a0=lna,a1=-b,则原曲线化为z= +a1x。 计算结果为y=2.35074x 上机实习题3 1.试分别用抛物线 和指数曲线 拟合下列数据 EMBED Equation.3 1 1.5 2 2.5 3 3.5 4 4.5 33.4 79.50 122.65 159.05 189.05 214.15 238.65 252.50 EMBED Equation.3 5 5.5 6 6.5 7 7.5 8 267.55 280.50 296.65 301.40 310.40 318.75 325.15 比较两个拟合函数的优劣。 2. 已知实验数据如下： EMBED Equation.3 1.0 2.5 3.5 4.0 3.8 1.50 26.0 33.0 试用形如 的抛物线进行最小二乘拟合。 实验四 数值积分 1.目的与要求： 1) 通过实际计算体会各种方法的精确度； 2) 会编写用龙贝格算法求定积分的程序。 2. 复化梯形公式与复化辛卜生公式的自适应算法 算法 复化辛卜生公式为 ，计算过程为： eq \o\ac(○,1)令 ； eq \o\ac(○,2)对 计算 。 程序与实例 例1 用复化辛卜生公式计算积分I= 。 程序如下： #include "stdafx.h" #include #include void main() { int i,n=2; float s; float f(float); float Simpson(float(*)(float),float,float,int); for(i=0;i<=2;i++) { s=Simpson(f,0,1,n); printf("s(%d)=%f\n",n,s); n*=2; } } float Simpson(float (*f)(float),float a,float b, int n) { int k; float s,s1,s2=0.0; float h=(b-a)/n; s1=f(a+h/2); for(k=1;k<=n-1;k++) { s1+=f(a+k*h+h/2); s2+=f(a+k*h); } s=h/6*(f(a)+4*s1+2*s2+f(b)); return s; } float f(float x) { return 1/(1+x*x); } 运行结果为 s(2)=0.785392 s(4)=0.785398 s(8)=0.785398 说明：本例运行了三次，当n= =8时，就与n= =4时有6位数字相同,若用复化梯形法计算，当n=512时有此结果。 例2 用复化辛卜生法计算积分 I= 。 只要将函数定义为 #include "stdafx.h" #include #include #include void main() { int i,n=2; float s; float f(float); float Simpson(float(*)(float),float,float,int); for(i=0;i<=2;i++) { s=Simpson(f,0,1,n); printf("s(%d)=%f\n",n,s); n*=2; } getch(); } float Simpson(float (*f)(float),float a,float b, int n) { int k; float s,s1,s2=0.0; float h=(b-a)/n; s1=f(a+h/2); for(k=1;k<=n-1;k++) { s1+=f(a+k*h+h/2); s2+=f(a+k*h); } s=h/6*(f(a)+4*s1+2*s2+f(b)); return s; } float f(float x) { if(x==0) return 1; else return sin(x)/x; } 运行结果为 s(2)=0.946087 s(4)=0.946083 s(8)=0.946083 3自适应梯形公式 算法 变步长梯形算法依据公式 。计算时可按如下步骤： eq \o\ac(○,1)输入精度 ； eq \o\ac(○,2)s=0; 1 2 3 ⑥ 4 ②。 例3 用变不长梯形法积分I= 。 程序如下： #include "stdafx.h" #include #include #include int n; void main() { int i; float s; float f(float); float AutoTrap(float(*)(float),float,float); s=AutoTrap(f,0.0,1.0); printf("T(%d)=%f\n",n,s); getch(); } float AutoTrap(float (*f)(float),float a,float b) { int i; float x,s,h=b-a; float t1,t2=h/2.0*(f(a)+f(b)); n=1; do { s=0.0; t1=t2; for(i=0;i<=n-1;i++) { x=a+i*h+h/2; s+=f(x); } t2=(t1+s*h)/2.0; n*=2; h/=2.0; } while(fabs(t2-t1)>1e-6); return t2; } float f(float x) { return 1/(1+x*x); } 运行结果为 T(512)=0.785398 例4 用变步长梯形法积分I= 。 只要将函数定义为 float f(float x) { if(x==0) return 1; else return sin(x)/x; } 运行结果为 T(512)= 0.946083 4.龙贝格算法 算法 程序与实例 例5 用龙贝格方法计算积分I= 。 程序如下： #include "stdafx.h" #include #include #include float f(float x) { return 1/(1+x*x); } float Romberg(float a,float b,float (*f)(float),float epsilon) { int n=1,k; float h=b-a,x,temp; float T1,T2,S1,S2,C1,C2,R1,R2; T1=(b-a)/2*((*f)(a)+(*f)(b)); while(1) { temp=0; for(k=0;k<=n-1;k++) { x=a+k*h+h/2; temp+=(*f)(x); } T2=(T1+temp*h)/2; if(fabs(T2-T1)(fabs(*(c+k*(n+1)+i)))) k=j; if(k!=i) for(j=i;j<=n;j++) {p=*(c+i*(n+1)+j); *(c+i*(n+1)+j)=*(c+k*(n+1)+j); *(c+k*(n+1)+j)=p; } for(j=i+1;j<=n-1;j++) { p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i)); for(t=i;t<=n;t++) *(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t)); } } for(i=n-1;i>=0;i--) { for(j=n-1;j>=i+1;j--) (*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j)); x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i)); } return x; } 程序结果： x[0]=-0.398234 x[1]=0.013795 x[2]=0.335144 Press any key to continue 3.矩阵直接三角分解法 算法 将方程组AX=b中的A 分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，则方程组AX=b化为解2个方程组Ly=b，Ux=y,具体算法如下： 1 对j=1，2，3，…,n计算 对i=2,3,…,n计算 2 对k=1，2，3，…,n: a. 对j=k,k+1,…,n计算 b. 对i=k,k+1,…,n计算 3 对k=2,3,…,n计算 ④ 对k=n-1,n-2,…,2,1计算 注:由于计算u的公式与计算y的公式形式上一样,故可直接对增广矩阵 [A|b]= 施行算法②,③,此时U的第n+1列元素即为y. 程序与实例 例 2 求解方程组Ax=b ,b= 。 程序如下： #include "stdafx.h" #include"stdio.h" #include"math.h" #include"conio.h" #include"malloc.h" void main() { float *x; int i; float a[4][5]={1.0, 2.0,-12.0,8.0,27.0, 5.0,4.0,7.0,-2.0,4.0,-3.0,7.0,9.0,5.0,11.0, 6.0,-12.0,-8.0,3.0,49.0}; float *DirectLU(float *,int); x=DirectLU((float *)a,4); for(i=0;i<=3;i++) printf("x[%d]=%f\n",i,x[i]); } float *DirectLU(float *u,int n) { int i,r,k; float *x; x=(float*)malloc(n*sizeof(float)); for(r=0;r<=n-1;r++) { for(i=r;i<=n;i++) for(k=0;k<=r-1;k++) *(u+i*(n+1)+i)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+i)); for(i=r+1;i<=n-1;i++) { for(k=0;k<=r-1;k++) *(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+r)); *(u+i*(n+1)+r)/=*(u+r*(n+1)+r); } } for(i=n-1;i>=0;i--) { for(r=n-1;r>=i+1;r--) *(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*x[r]; x[i]=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i)); } return x; } 结果为： x[0]=28.665836 x[1]=-0.849526 x[2]=-0.272213 x[3]=-0.404167 Press any key to continue 上机实验题 5 1. 用列主元消去法解方程组: 2. 用LU分解法解方程组AX=b,其中 3.编写用改进的平方根法解方程组AX=b的程序,并解下列方程组: , 4.编写用追赶法解三对角线性方程组的程序,并解下列方程组: 实验六 解线性方程组的迭代法 1.目的与要求： 1） 熟悉求解线性方程组的有关理论哈方法。 2） 会编制雅可比迭代和高斯—塞得尔迭代法。 3） 通过实际计算，进一步了解各算法的优缺点，选择合适的数值方法。 2．雅可比迭代法 算法 设方程组AX=b的系数矩阵的对角元素 M为迭代次数容许的最大值, 为容许误差. 1 取初始向量 令k=0; 2 对 计算 3 如果 则输出 ,结束;否则执行④, 4 如果 则不收敛,终止程序;否则 转②. 程序与实例 例4.用雅克比迭代法解方程组 程序如下: #include "stdafx.h" #include "iostream" #include"stdio.h" #include"math.h" #include"conio.h" #include"malloc.h" #include #define EPS 1e-6 #define MAX 100 float *Jacobi(float a[3][4],int n) { float *x,*y,s; double epsilon; int i,j,k=1; x=(float *)malloc(n*sizeof(float)); y=(float *)malloc(n*sizeof(float)); for(i=0;iEPS) ; if(k>=MAX) { return y; } for(i=0;i=N) { printf("迭代发散\n"); exit(1); } nu++; } while(fabs(dx)<1e-6); return x; } 结果为： x[0]=2.500000 x[1]=2.090909 x[2]=1.227273 Press any key to continue 上机实验 6 1.分别用雅可比迭代法与高斯-塞德尔迭代法解下列方程组: EMBED Equation.3 实验七 方程求根 1.目的与要求： 1） 通过对二分法与牛顿迭代法作编程练习与上机运算，进一步体会二分法与牛顿迭代法的不同特点； 2） 编写割线迭代法的程序，求非线性方程的解，并与牛顿迭代法作比较。 2.二分法 算法 给定区间[a,b]，并设f(a)与f(b)符号相反，取为ε根的容许误差，δ为|f(x)|的容许误差。 1 令c=(a+b)/2 2 如果(c-a)<ε或|f(c)|<δ，则输出c，结束；否则执行③， 3 如果f(a)f(c)>0，则令a:=c；否则b:=c，重复①，②，③。 程序与实例 例1 求方程 在1.5附近的根。 程序如下： #include "stdafx.h" #include #include #include #define EPS 5e-6 #define DELTA 1e-6 float Bisection(float,float,float(*f)(float)); float f(float); void main() { float a=1,b=2; float x; x=Bisection(a,b,f); printf("\nThe root of the equation is %f\n",x); } float Bisection(float a,float b,float(*f)(float)) { float c,fc,fa=f(a),fb=f(b); int n=1; printf("二分次数 c\t\t f(c)\n"); while(1) { if(fa*fb>0) exit(0); c=(a+b)/2,fc=f(c); if(fabs(fc) #include #include #define N 100 #define EPS 1e-6 #define ETA 1e-8 void main() { float f(float); float f1(float); float x0,y0; float Newton(float (*)(float),float (*)(float),float); printf("Please input x0\n"); scanf("%f",&x0); printf("x(0)=%f\n",x0); y0=Newton(f,f1,x0); printf("\nThe root of the equation is x=%f\n",y0); } float Newton(float (*f)(float),float (*f1)(float),float x0) { float x1,d; int k=0; do { x1=x0-f(x0)/f1(x0); if((k++>N)||(fabs(f1(x1))EPS&&fabs(f(x1))>ETA); return x1; } float f(float x) { return x*x*x+x*x-3*x-3; } float f1(float x) { return 3.0*x*x+2*x-3; } 输出结果如下： 若取出之x(0)=1.000000,则 x(1)=3.000000, x(2)=2.200000, x(3)=1.830151, x(4)=1.737795 x(5)=1.732072, x(6)=1.732051, x(7)=1.732051 原方程的根为x=1.732051。 若取初值x(0)=1.500000，则 x(1)=1.7777778, x(2)1.733361, x(3)=1.732052, x(4)=1.732051 原方程的根为x=1.732051. 若取初值x(0)=2.500000，则 x(1)=1.951807, x(2)=1.758036, x(3)=1.732482, x(4)=1.732051 原方程的根为x=1.732051 说明：上面程序取3个不同初值，得到同样的结果，但迭代次数不同。初值越接近所求的根，迭代次数越少。 例3 求方程 在0.3附近的根。 只要将上述程序中的函数及导函数改为下面的函数： float f(float x) { return 3*x+sin(x)-exp(x); } float f1(float x) { return 3+cos(x)-exp(x); } 输出结果如下： 若取初值x(0)=0.50000，则 X(1)=