数值分析上机题概况

1、数值分析上机题习题117.（上机题）舍入误差与有效数设Sn21 ，其精确值为-3 。j j -12 2 N N 1（1 ）编制按从大到小的顺序 Sn川-，计算Sn的通用程序。2131N -1（2）编制按从小到大的顺序 SN二亠，计算Sn的通用程序。N -1 （N 1） 12 -1（3） 按两种顺序分别计算 §02，Sd，S106，并指出有效位数。（编制程序时用单精度）（4 ）通过本上机题你明白了什么？按从大到小的顺序计算Sn的通用程序为:按从小到大的顺序计算Sn的通用程序为:#in clude<iostream.h>#in clude<iostream.h>f

2、loat sum(float N) float j,s,sum=0;for(j=N;j>=2;j-) s=1/(j*j-1); sum+=s;return sum;float sum(float N)float j,s,sum=0; for(j=2;j<=N;j+) s=1/(j*j-1); sum+=s;return sum;从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大S1020.7400490.740050.74004965S1040.7498520.74990.749944S1060.7498520.7499990.74999936通过本上机题，看出按两

3、种不同的顺序计算的结果是不相同的, 的值与精确值有较大的误差， 计算得到的结果的有效位数少。 致计算结果的精度有所降低，按从大到小的顺序计算 而按从小到大的顺序计算的值与精确值吻合。从大到小的顺序计算机在进行数值计算时会出现“大数吃小数”的现象，导 我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。习题220.(上机题)Newton迭代法(1 )给定初值X。及容许误差:，编制Newt on法解方程f(x) =0根的通用程序。(2)给定方程 f (x) = X3 / 3 - x = 0,易知其有三个根 X； - - 3， x2 = 0， x3 二 祸。1.由Ne

4、wt on方法的局部收敛性可知存在心 0，当x时，Newton迭代序列收敛于根x2。试确定尽可能大的。2 .试取若干初始值，观察当 X0(-：，-1) , (-',、) , C' ,1) , (1,：)时Newton序列是否收敛以及收敛于哪一个根。(3)通过本上机题，你明白了什么？ 解：(1)编制的通用程序：#in clude<iostream.h>#in clude<math.h>#defi ne eps 0.000001float f(float x)float f; f=x*x*x/3-x; return(f);float df(float x)计

5、算f(x)的导函数float df;df=x*x-1;/f(x)导函数的表达式；return (df);void mai n(void)(2)计算迭代序列收敛于根x2的尽可能大的/给定容许误差定义函数f(x)/f(x)的表达式；定义函数df(x),float x0,x1,a;int k=0;cout<<"请输入初值x0:"cin> >x0;doa=-f(x0)/df(x0); x1=x0+a;k+;x0=x1;while(fabs(a)>eps); cout<<k<<'t'<<x0;输出迭代的

6、次数和根值:的函数为：#in clude<iostream.h>#in clude<math.h> void delay(i nt n) for(n=10000 ;n >0; n-); #defi ne eps 0.000001 float f(float x)float f; f=x*x*x/3-x;/定义延时函数定义函数f(x)/f(x)的表达式；return(f);float df(float x)/定义函数 df(x)，计算 f(x)的导函数float df;df=x*x-1;/f(x)导函数的表达式；return (df);int judgeme nt(

7、float z)int coun t=5;float x0,x1,type,type1;xO=z;while(co un t->0)x仁 xO-f(xO)/df(xO);type=fabs(x1);type 仁 fabs(x1-x0);/ 调试值用cout<<"co un t="<<co un t<<'t'<<"type= "<<type<<'t'<<"type 仁"<<type1<<

8、9;n'if(fabs(x1-x0)<eps)return 1;x0=x1;delay(30000);调试值用return 0;void mai n(void)float delta=0;int flag=1;while(flag=1)cout<<"方程的根为:"<<'n' delta+=eps;flag=judgeme nt(delta);cout<<"输出方程根收敛的区间 值:n"cout<<delta-eps; /输出收敛的区间值当步长为0.001时，程序计算出的 3的为

9、S =0.774，即在区间（-0.774, 0.774）内迭代序列 收敛于0。对于不同得初始值收敛于不同的根 ，x0在（-8, -1 ）内收敛于X；，在（-0.774 , 0.774）内 收敛于x2,在（1 , + 内收敛于 x3 ,但在内（0.774, 1）和（一1 , 0.774）均可能收敛于 X；和x3。 X； , x2 , x3分别为方程的精确解。分析：对于不同的初值，迭代序列会收敛于不同的根，所以在某个区间内求根对于初值的选取有很大的关系。产生上述结果的原因是区间不满足大范围收敛的条件。习题339.(上机题)列主元Gauss消去法VT=( -15,27, -23,0,-20,12,

10、-7,7,10)对于某电路的分析，归结为求解线性方程组RI=V。31-13000-10000-1335-90-1100000-931-100000000-1079-30000-9000-3057-70-500000-747-300000000-3041000000-50027-2000-9000-229 _T(1) 编制解n阶线性方程组 Ax=b的列主元三角分解法的通用程序；(2 )用所编制的程序解线性方程组RI=V，并打印出解向量，保留五位有效数;II输出数组aII进行第一行和第一列元素的求取(3)本编程之中，你提高了哪些编程能力? 程序为：|#in clude<iostream.h&

11、gt;#in clude<math.h>void mai n(void)in t i,j,n ,k,q;float a1011,s10,s110;cout<<"请输入n的值:"；cin»n;cout<<"输入数组 a:"<<endl; for(i=1;i<=n ;i+)for(j=1;j<=(n+1);j+)cin>>aij;II给矩阵 a 赋值for(i=1;i<=n ;i+)for(j=1;j<=(n+1);j+) cout<<aij<&l

12、t;'t'cout<<'n'cout<<'IIint t=1;for(i=1;i<=n ;i+) si=ai1; float max=fabs(s1);for(i=2;i<=n ;i+) if(fabs(si)>max) max=fabs(si); t=i;for(j=1;j<=(n+1);j+)float b=a1j;a1j=atj;atj=b;进行第一列主元互换for(i=2;i<=n ;i+)ai1=ai1/max;/第一列除以 a11for(i=1;i<=n ;i+) for(j=1;j

13、<=(n+1);j+) cout<<aij<<'t' cout<<'n：进行第k步分解'/输出进行第一步变换的数组a cout<<" for(k=2;k<=n; k+)for(i=k;i<=n ;i+)float sum=0; for(q=1;q<k;q+) sum+=aiq*aqk; s1i=aik-sum;in t l=k;float m=fabs(s1k);for(i=k;i<=n ;i+)/比较第k步分解的第k列值的大小 if(fabs(s1i)>m)m=fab

14、s(s1i);l=i;for(j=1;j<=n+1;j+)float s2=akj; akj=alj;结果：方程的解为：x1= -0.28923，-0.71281，/返回行值交换两行元素x2= 0.34544，x3=x5= -0.43040， x6=x4= -0.22061， alj=s2;for(j=k;j<=n+1;j+)/算出第k行行元素的值float sum1=0;for(q=1;q<k;q+)sum1+=akq*aqj; akj=akj-sum1;for(i=k+1;i<=n ;i+)算出第k列列元素的值float sum2=0;for(q=1;q<k;

15、q+)sum2+=aiq*aqk; aik=(aik-sum2)/(akk);/第k步分解结束for(i=1;i<=n ;i+)for(j=1;j<=(n+1);j+)cout<<aij<<'t'cout<<'n'/输出改变后的数组float x10;for(i=n-1;i>=1;i-)xn =a n n+1/a nn;float sum3=0;for(j=i+1;j<=n ;j+)sum3+=aij*xj;xi=(ai n+1-sum3)/aii;/回代过程for(i=1;i<=n ;i+)co

16、ut<<'x'<<i<<'='<<xi<<e ndl;输出解向量0.15431，x7= -0.057823，x8= 0.20105，x9=0.29023。分析：我感觉是提高了查错误点的能力和编循环语句的能力，即利用很规整的迭代公式 进行编程。另外列主元三角分解法的阶梯步骤有了更深的了解习题437.(上机题)3次样条插值函数(1) 编制求第一型3次样条插值函数的通用程序(2) 已知汽车曲线型值点的数据如下：Xi012345678910Yi2.513.304.044.705.225.545.785.405.

17、575.705.80端点条件为y0 =0.8，y；0=O.2。用所编制程序求车门的 3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,9)。解：通用程序：#in clude<iostream.h>void mai n(void)float x11; 存放数组 xjfloat y11; 存放数组 yjfloat h11; 存放数组 hjfloat u11; 存放数组 ujfloat v11;/ 存放数组 vjfloat d11; 存放数组 djfloat M11; 存放数组 Mjfloat b11; 存放数组 bj floatt11,l11,yy11,s4,aa1,aa

18、2,aa3,aa4; float s110;int i,j,n;float xx;/x为区间值/将初值初始化cout<<"请输入n的值:n" cin»n;cout<<"输入数组 x:n" for(i=0;i<=n ;i+)ci n> >xi; cout<<"输入数组 y:n" for(i=0;i<=n ;i+)ci n> >yi;/输入端点值float df2;cout<<"输入两个端点值:n" for(i=0;i<

19、2;i+)ci n> >dfi;/利用书本上的算法求出所需要的值 /求出hj的值for(j=0;j<=n _1;j+)hj=xj+1-xj;cout<<'h'<<T<<j<<T<<'='<<hj<<'t'cout<<e ndl;/求出uj和vj的初值v0=1;u n=1;for(j=1;j<=n _1;j+)uj=hj-1/(hj-1+hj); vj=hj/(hj-1+hj);/求出dj的值for(j=1;j< n;j+)

20、dj=6*(yj+1-yj)/hj-(yj-yj-1)/h j-1)/(hj+hj-1);d0=6*(y1-y0)/h0-df0)/h0;d n=6*(df1-(y n-y n-1)/h n-1)/h n-1 ;for(j=1;j<=n ;j+)cout<<'u'<<''<<j<<''<<'='<<uj<<'t'cout<<e ndl;for(j=0;j< n;j+)cout<<'v

21、9;<<''<<j<<''<<'='<<vj<<'t'S(x):0,1时；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)1.2 时；S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1)2.3 时；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(

22、x-2)3.4 时；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3)4.5 时；S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4)5.6 时；S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5)6.7 时；S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6)7.8 时；S

23、(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7)8.9 时；S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8)9.10 时；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9)S(1.5)=3.67843S (2.5)=4.38147S(4.5)=5.38328S(5.5)=5.7237S(7.5)=5.42989S(8.5)=5.65976cout&

24、lt;<e ndl; for(j=0;j<=n ;j+) cout<<'d'<<T<<j<<T<<'='<<dj<<'t' cout<<e ndl;/利用书本上的追赶法求解方程组 for(i=0;i<=n ;i+)bi=2; cout<<e ndl; tO=bO; yyO=dO；/消元过程 for(i=1;i<=n ;i+)li=ui/ti-1; ti=bi-li*vi-1; yyi=di-li*yyi-l;/回代过

25、程M n=yy n/t n; for( i=n-1;i>=0;i_) Mi=(yyi-vi*Mi+1)/ti;(2) 编制的程序求车门的3次样条插值函数x属于区间 x属于区间 x属于区间 x属于区间 x属于区间 x属于区间 x属于区间 x属于区间 x属于区间 x属于区间S(0.5)=2.90856 S(3.5)=4.98819 S(6.5)=5.59441 S(9.5)=5.7323将Mj的值输出for(i=0;i<=n ;i+) cout<<'M'<<T<<i<<''<<'=

26、9;<<Mi<<e n dl;/输出插值多项式的系数for(j=0;j< n;j+)s0=yj;s1=(yj+1-yj)/hj-(Mj/3+Mj+1/6) *hj;s2=Mj/2;s3=(Mj+1-Mj)/(6*hj);cout<<"当 x 的值在区间 "<<'x'<<T<<j<<T<<"到"<<'x'<<T<<(j+1)<<''<<"时

27、,输出插值多项 式的系数:n"for(int k=0;k<4;k+)cout<<'s'<<T<<k<<T<<'='<<sk<<'t'cout<<e ndl;习题623.（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2） 编制AB4方法的通用程序（由 RK4提供初值）；（3） 编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（4） 编制带改进的 AB 4-AM 4预测校正方法的通用程序（由RK4提供初值

28、）；（5）对于初值问题V = -x2y2（0兰 x兰1.5）»（0）=3取步长h =0.1 ,应用（1） （ 4）中的四种方法进行计算，并将计算结果和精确解3y（x） =3/（1 x ）作比较；（6）通过本上机题，你能得到哪些结论？解：程序为：#in clude<iostream.h>#in clude<fstream.h>#in clude<stdlib.h>#in clude<math.h>ofstream outfile("data.txt");此处定义函数f(x,y)的表达式/用户可以自己设定所需要求得函数表

29、达式double f1(double x,double y)double f1;f1=(_1)*x*x*y*y;return f1;/此处定义求函数精确解的函数表达式double f2(double x)double f2;f2=3/(1+x*x*x);return f2;/此处为精确求函数解的通用程序void accurate(double a,double b,double h)double x100,accurate100;x0=a;int i=0;outfile<<"输出函数准确值的程序结果:n"doxi=x0+i*h; accuratei=f2(xi)

30、;outfile<<"accurate"<<i<<"="<<accurat ei<<'n'i+;while(i<(b-a)/h+1);/此处为经典Runge-Kutta公式的通用程序 void RK4(double a,double b,double h,double c)int i=0;double k1,k2,k3,k4;double x100,y100;y0=c;x0=a;outfile<<"输出经典 Runge-Kutta公式的程序结果：n&q

31、uot;doxi=x0+i*h;k1= f1(xi,yi);k2=f1(xi+h/2),(yi+h*k1/2);k3=f1(xi+h/2),(yi+h*k2/2);k4=f1(xi+h),(yi+h*k3); yi+1=yi+h*(k1+2*k2+2*k3+k4)/6; outfile<<"y"<<T<<i<<"="<<yi<<'n'i+;while(i<(b-a)/h+1);此处为4阶Adams显式方法的通用程序void AB4(double a,double

32、 b,double h,doublec)double x100,y100,y1100;double k1,k2,k3,k4;y0=c;x0=a;outfile<<"输出4阶Adams显式方法的程序结果:n"for(int i=0;i<=2;i+)xi=x0+i*h;k1= f1(xi,yi);k2=f1(xi+h/2),(yi+h*k1/2);k3=f1(xi+h/2),(yi+h*k2/2);k4=f1(xi+h),(yi+h*k3);yi+1=yi+h*(k1+2*k2+2*k3+k4)/6;int j=3;y10=y0;y11=y1;y12=y2;y

33、13=y3;doxj=x0+j*h;y1j+1=y1j+(55*f1(xj,y1j)-59*f1(x j-1,y1j-1)+37*f1(xj-2,y1j-2)-9*f1(xj-3 ,y1j-3)*h/24;n'j+;while(j<(b-a)/h+1);/主函数void mai n(void)double a,b,h,c;cout<<"输入上下区间、步长和初始 值:n"cin> >a>>b»h»c;accurate(a,b,h);RK4(a,b,h,c);AB4(a,b,h,c);结果为：y（xj和精确值

34、和由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值 数值解的误差：Y10=3y11=2.997y12=2.97619y1 5=2.66467y1 6=2.4652y17=2.23308y110=1.50219y111=1.28876y112=1.10072y115=0.685335通过本上机题我明白了各种求微分方程的数值方法，y13=2.92113 y18=1.98495 y113=0.93871y14=2.81839y19=1.73704y114=0.801135经典Runge-Kutta公式，AB4方法以及iXiyiy(x i)|y(x i )-y i |003301

35、0.12.9972.9971.87138e-00720.22.976192.976193.91665e-00730.32.921132.921137.58342e-00740.42.819552.819551.61101e-00650.52.666662.666673.17735e-00660.62.46712.467115.00551e-00670.72.23382.23385.77233e-00680.81.984121.984134.12954e-00690.91.735111.735111.15554e-007101.01.500011.55.80668e-006111.11.2870

36、11.2871.13075e-005121.21.099721.099711.54242e-005131.30.9383970.938381.77272e-005141.40.80130.8012821.83754e-005151.50.6857320.6857141.78e-005由AB4方法得出的结果为:Runge-Kutta公式的精度，四AB4-AM4预测校正方法求解公式的精度是不同的。其中经典 阶Adams显式方法（AB4）具有4阶精度。习题710.抛物线方程Crank-Nicolson格式(1)编制用Crank-Nicolson格式求抛物线方程<t<T)2 2.:u/ -

37、 a u/ ： x = f(x,t) (0<x<1, 0Y u(x,0) =®(x)(0ExW1).u(0,t) =a(x) , u (1,t ) =P(t) (0<t <1)数值解的通用程序。(2) 就 a=1, f(x,t)=0,x =exp(x), xt =exp(t), : t =exp(1+t),M=40,N=40,输入点(0.2，1.0)，( 0.4，1.0)，( 0.6，1.0)，( .8，1.0)4 点处 u ( x,t )的近似值。(3) 已知所给方程的精确解为u(x,t)=exp(x+t),将步长反复二分，从(0.2，1.0)，(0.4，1

38、.0)，( 0.6，1.0)，( 0.8，1.0)4点处精确解与数值解的误差观察当步长缩小一半时，误差以什么规律缩小。算法:本题选择空间步长h=0.025,时间步长 Z =0.0251. 置初值;求解点(I,k+1)2. 通过点(i-1,k),(i-1,k+1),(i,k),(i+1,k),(i+1,k+1)的值；3. 输出结果原程序：#in clude <iostream>#in clude <math.h> float h=0.025,k=0.025;int m=40;i nt n=40;float y4040,r=a*k/(h*h);void In put()int i,j;cout<<"Loadi ng In put Data."<<e ndl;for(i=0;i<m;i+)for(j=0;j< n;j+)if (i=j) aij=1+r;for(j=0;j< n;j+)if (j=i+1)|(i=j+1)