数值分析上机题概况

1、数值分析上机题习题 117（上机题）舍入误差与有效数设 SN21 ，其精确值为 1 3 1 1 。N j 2 j2 1 2 2 N N 1（1）编制按从大到小的顺序 SN212121 ，计算 SN 的通用程序。22 1 32 1 N 2 1 （2）编制按从小到大的顺序 SN211 221 ，计算 SN 的通用程序。N 2 1 （N 1）2 122 1（3）按两种顺序分别计算 S102 ， S104 ， S106 ，并指出有效位数。 （编制程序时用单精度）（4）通过本上机题你明白了什么？按从大到小的顺序计算SN 的通用程序为：按从小到大的顺序计算SN 的通用程序为：#include<ios

2、tream.h>#include<iostream.h>16float 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.749944

3、S1060.7498520.7499990.74999936通过本上机题， 看出按两种不同的顺序计算的结果是不相同的， 的值与精确值有较大的误差， 计算得到的结果的有效位数少。 致计算结果的精度有所降低，按从大到小的顺序计算 而按从小到大的顺序计算的值与精确值吻合。 从大到小的顺序 计算机在进行数值计算时会出现“大数吃小数”的现象，导 我们在计算机中进行同号数的加法时， 采用绝对值较小者先加的算法，其结果的相对误差较小。习题 220(上机题) Newton 迭代法(1)给定初值 x0及容许误差 ，编制 Newton 法解方程 f(x) 0根的通用程序。(2)给定方程 f(x) x3/3 x 0

4、，易知其有三个根 x13，x2 0， x33。1由 Newton 方法的局部收敛性可知存在0，当 x0 ( , ) 时， Newton 迭代序列收敛于根 x2 。试确定尽可能大的。2试取若干初始值，观察当x0 ( , 1)，( 1, )，( , )，( ,1) ， (1, ) 时Newton 序列是否收敛以及收敛于哪一个根。 (3)通过本上机题，你明白了什么？ 解：(1)编制的通用程序： #include<iostream.h> #include<math.h> #define eps 0.000001 float f(float x) float f; f=x*x*x

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

6、gt;eps);cout<<k<<'t'<<x0; /输出迭代的次数和根值 的函数为：#include<iostream.h>#include<math.h> void delay(int n) for(n=10000;n>0;n-); #define 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) 的导函

7、数float df;df=x*x-1; /f(x) 导函数的表达式 ; return (df);int judgement(float z)int count=5;float x0,x1,type,type1;x0=z;while(count->0)x1=x0-f(x0)/df(x0); type=fabs(x1);type1=fabs(x1-x0); / 调试值用 cout<<"count="<<count<<'t'<<"type= "<<type<<'

8、;t'<<"type1="<<type1<<'n'if(fabs(x1-x0)<eps)return 1;x0=x1;delay(30000);/调试值用return 0;void main(void)float delta=0;int flag=1;while(flag=1)cout<<" 方程的根为 :"<<'n' delta+=eps; flag=judgement(delta);cout<<" 输 出 方 程 根 收 敛

9、 的 区 间 值:n"cout<<delta-eps; / 输出收敛的区间值当步长为 0.001 时，程序计算出的 的为 =0.774 ，即在区间（ -0.774，0.774）内迭代序列 收敛于 0。对于不同得初始值收敛于不同的根 , x0在（ -， -1）内收敛于 x1* ，在（ -0.774，0.774）内 收敛于 x2，在（ 1，+）内收敛于 x3 ，但在内（ 0.774，1）和（ 1，0.774）均可能收敛于 x1*和x3。 x1* ， x2， x3分别为方程的精确解。分析： 对于不同的初值， 迭代序列会收敛于不同的根， 所以在某个区间内求根对于初值的选 取有很大

10、的关系。产生上述结果的原因是区间不满足大范围收敛的条件。习题 339 （上机题） 列主元 Gauss 消去法对于某电路的分析，归结为求解线性方程组RI=V。VT=(-15,27, -23,0,-20,12,-7,7,10)其中,31-13000-10000-1335-90-1100000-931-100000000-1079-30000-9000-3057-70-500000-747-300000000-3041000000-50027-2000-9000-229T1）编制解 n 阶线性方程组 Ax=b 的列主元三角分解法的通用程序；2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留

11、五位有效数；/ 输出数组 a<<'n'/ 进 行 第 一 行 和 第 一 列 元 素 的 求 取(3)本编程之中，你提高了哪些编程能力？ 程序为：#include<iostream.h> #include<math.h> void main(void) int i,j,n,k,q;float a1011,s10,s110; cout<<" 请输入 n 的值 :"cin>>n;cout<<" 输入数组 a:"<<endl; for(i=1;i<=n;i

12、+) for(j=1;j<=(n+1);j+) cin>>aij;/给矩阵 a 赋值for(i=1;i<=n;i+) for(j=1;j<=(n+1);j+) cout<<aij<<'t'cout<<'n'cout<<"'/int 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&

13、lt;=(n+1);j+)<<'n'/进行第 k 步分解'/返回行值/交换两行元素x2= 0.34544 ， x3=x5= -0.43040 ， x6=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<=(n+1);j+)cout<<aij<<'t'cout<<'n'/输出进行第一步变换的数组 a cout<<

14、;" 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;int l=k;float m=fabs(s1k);for(i=k;i<=n;i+)/比较第 k步分解的第 k 列值的大小if(fabs(s1i)>m)m=fabs(s1i);l=i;for(j=1;j<=n+1;j+)float s2=akj; akj=alj; 结果：方程的解为：x1= -0.28923 ，-0.71281，x4= -0.22061 ， alj=s2; f

15、or(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;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<<

16、'n' / 输出改变后的数组/ 输出解 ''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

17、39;/float x10; for(i=n-1;i>=1;i-) xn=ann+1/ann; float sum3=0;for(j=i+1;j<=n;j+) sum3+=aij*xj;xi=(ain+1-sum3)/aii;/ 回代过程for(i=1;i<=n;i+) cout<<'x'<<i<<'='<<xi<<endl;/ 输出解向量0.15431，x7= -0.057823 ，x8= 0.20105 ，x9= 0.29023。分析：我感觉是提高了查错误点的能力和编循环语句的能力

18、， 即利用很规整的迭代公式 进行编程。 另外列主元三角分解法的阶梯步骤有了更深的了解习题 437. (上机题) 3 次样条插值函数1)编制求第一型 3 次样条插值函数的通用程序2) 已知汽车曲线型值点的数据如下：xi012345678910yi2.513.304.044.705.225.545.785.405.575.705.80端点条件为 y0' =0.8， y1' 0 =0.2 。用所编制程序求车门的 3 次样条插值函数 S(x), 并打印出S(i+0.5)(i=0,1,9) 。解：通用程序： #include<iostream.h> void main(voi

19、d) float x11;/ 存放数组 xjfloat y11;/ 存放数组 yjfloat h11;/ 存放数组 hjfloat u11;/ 存放数组 ujfloat v11;/ 存放数组 vjfloat d11;/ 存放数组 djfloat M11;/ 存放数组 Mjfloat b11;/ 存放数组 bjfloat t11,l11,yy11,s4,aa1,aa2,aa3,aa4;float s110;int i,j,n;float xx;/x 为区间值/将初值初始化cout<<" 请输入 n 的值 :n" cin>>n;cout<<&

20、quot; 输入数组 x:n" for(i=0;i<=n;i+)cin>>xi;cout<<" 输入数组 y:n"for(i=0;i<=n;i+)cin>>yi;/输入端点值float df2;cout<<" 输入两个端点值 :n" for(i=0;i<2;i+)cin>>dfi;/利用书本上的算法求出所需要的值 /求出 hj 的值 for(j=0;j<=n-1;j+)hj=xj+1-xj;cout<<'h'<<'

21、'<<j<<''<<'='<<hj<<'t'cout<<endl;/求出 uj和 vj 的初值v0=1;un=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+)dj=6*(yj+1-yj)/hj-(yj-yj-1)/h j-1)/(hj+hj-1);d0=6*(y1-y0)/h0-df0)/h0;dn=6*(df1-(yn-yn-1)/hn-1)/h

22、n-1 ;for(j=1;j<=n;j+)cout<<'u'<<''<<j<<''<<'='<<uj<<'t'cout<<endl;for(j=0;j<n;j+)cout<<'v'<<''<<j<<''<<'='<<vj<<'t'S(x) ：0,1

23、 时； 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)(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.38

24、6786(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(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)

25、-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<<endl; for(j=0;j<=n;j+) cout<<'d'<<''<<j<<&

26、#39;'<<'='<<dj<<'t' cout<<endl; /利用书本上的追赶法求解方程组 for(i=0;i<=n;i+)bi=2; cout<<endl;t0=b0; yy0=d0; /消元过程 for(i=1;i<=n;i+) li=ui/ti-1; ti=bi-li*vi-1; yyi=di-li*yyi-1; /回代过程 Mn=yyn/tn; for(i=n-1;i>=0;i-)Mi=(yyi-vi*Mi+1)/ti;(2)编制的程序求车门的 3 次样条插值函数

27、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'<<''<<i<<''<<'='<<Mi<<en dl;/输出插值多项式的系数for(j=0;j<n;j+)s0=

28、yj;s1=(yj+1-yj)/hj-(Mj/3+Mj+1/6) *hj;s2=Mj/2; s3=(Mj+1-Mj)/(6*hj);cout<<" 当 x 的 值 在 区 间 "<<'x'<<''<<j<<''<<" 到 "<<'x'<<''<<(j+1)<<''<<" 时 ,输出插值多项 式的系数 :n"fo

29、r(int k=0;k<4;k+)cout<<'s'<<''<<k<<''<<'='<<sk<<'t'cout<<endl;习题 623（上机题） 常微分方程初值问题数值解（1）编制 RK4 方法的通用程序；（2）编制 AB4 方法的通用程序（由 RK 4提供初值）；（3）编制 AB 4-AM 4预测校正方法的通用程序（由RK4提供初值）；（4）编制带改进的 AB 4-AM 4预测校正方法的通用程序（由 RK 4提

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

31、表达式 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=f

32、2(xi);outfile<<"accurate"<<i<<"="<<accuratei<<'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 公式

33、的程序结果 :n"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"<<""<<i<<"="<<yi<<'n' i+;while(i<(b-a)/h+1);/此处为 4 阶 Adams 显式方法的

34、通用程序 void AB4(double a,double b,double h,double c)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*

35、k3+k4)/6;int j=3;y10=y0;y11=y1;y12=y2;y13=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 main(void)double a,b,h,c;cout<<" 输 入上 下区 间、步 长和 初始 值:n"cin>>a>>b>>h>>c; accu

36、rate(a,b,h); RK4(a,b,h,c); AB4(a,b,h,c);结果为：y(xi)和精确值和由经典 Runge-Kutta 公式得出的结果列在下面的表格中，以及精确值 数值解的误差：Y10=3y11=2.997y12=2.97619y15=2.66467 y16=2.4652y17=2.23308y110=1.50219 y111=1.28876 y112=1.10072 y115=0.685335 通过本上机题我明白了各种求微分方程的数值方法，y13=2.92113y18=1.98495y113=0.93871y14=2.81839y19=1.73704y114=0.8011

37、35经典 Runge-Kutta 公式， AB4 方法以及ixiyiy(x i )|y(x i )-y i |0033010.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.7

38、35111.735111.15554e-007101.01.500011.55.80668e-006111.11.287011.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 阶精度。习题 7

39、10.抛物线方程 Crank-Nicolson 格式( 1)编制用 Crank-Nicolson 格式求抛物线方程t T)22u/ t -a u/ x = f(x,t) (0<x<1, 0u(x,0) = (x) (0 x 1)u(0,t) =x ， u(1,t )= t (0<t 1)数值解的通用程序。2 )就 a=1, f(x,t)=0,x =exp(x), t =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) 已知所给方

40、程的精确解为 u(x,t)=exp(x+t), 将步长反复二分， 从( 0.2 ，1.0 )，(0.4 ， 1.0)，(0.6，1.0)，(0.8 ，1.0)4点处精确解与数值解的误差观察当步长缩小一半时，误差以什么规律缩小。算法：本题选择空间步长h=0.025, 时间步长 =0.0251. 置初值；求解点 (I,k+1)2. 通过点 (i-1,k),(i-1,k+1),(i,k),(i+1,k),(i+1,k+1)的值；3. 输出结果原程序： #include <iostream>#include <math.h> float h=0.025,k=0.025;int m=40;int n=40;float y4040,r=a*k/(h*h);void Input()int i,j;cout<<"Loading Input Data."<<endl; for(i=0;i<