东南大学数值分析上机题C++参考答案.DOC

1、数值分析上机题姓名：陈作添 学号：040816习题120（上机题）舍入误差与有效数设，其精确值为。（1）编制按从大到小的顺序，计算的通用程序。（2）编制按从小到大的顺序，计算的通用程序。（3）按两种顺序分别计算，并指出有效位数。（编制程序时用单精度）（4）通过本上机题，你明白了什么？第 6 页 共 16 页按从大到小的顺序计算的通用程序为：#includefloat sum(float N)float j,s,sum=0;for(j=2;j=N;j+) s=1/(j*j-1);sum+=s;return sum; 按从小到大的顺序计算的通用程序为：#includefloat sum(float

2、 N)float j,s,sum=0;for(j=N;j=2;j-)s=1/(j*j-1);sum+=s;return sum;从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.7400490.740050.740049650.7498520.74990.7499440.7498520.7499990.74999936通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。从大到小的顺序计算得到的结果的有效位数少。计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低

3、，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。习题220（上机题）Newton迭代法（1）给定初值及容许误差，编制Newton法解方程根的通用程序。（2）给定方程，易知其有三个根，。1由Newton方法的局部收敛性可知存在，当时，Newton迭代序列收敛于根。试确定尽可能大的。2试取若干初始值，观察当，时Newton序列是否收敛以及收敛于哪一个根。（3）通过本上机题，你明白了什么？解：（1）编制的通用程序：#include#include#define eps 0.000001 /给定容许误差float f(float x) /定义函数f(x)float

4、 f;f=x*x*x/3-x; /f(x)的表达式;return(f);float df(float x) /定义函数df(x)，计算f(x)的导函数float df;df=x*x-1; /f(x)导函数的表达式;return (df);void main(void)float x0,x1,a;int k=0;coutx0;doa=-f(x0)/df(x0);x1=x0+a;k+;x0=x1;while(fabs(a)eps);coutktx0; /输出迭代的次数和根值（2）计算迭代序列收敛于根的尽可能大的的函数为：#include#includevoid delay(int n) /定义延时

5、函数for(n=10000;n0;n-);#define eps 0.000001float f(float x) /定义函数f(x)float f;f=x*x*x/3-x; /f(x)的表达式;return(f);float df(float x) /定义函数df(x)，计算f(x)的导函数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)

6、; type1=fabs(x1-x0); /调试值用coutcount=countttype=typettype1=type1n;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输出方程根收敛的区间值:n;coutdelta-eps; /输出收敛的区间值当步长为0.001时，程序计算出的的为=0.774，即在区间（-

7、0.774，0.774）内迭代序列收敛于0。对于不同得初始值收敛于不同的根, 在（-，-1）内收敛于，在（-0.774，0.774）内收敛于，在（1，+）内收敛于，但在内（0.774，1）和（1，0.774）均可能收敛于和。，分别为方程的精确解。分析：对于不同的初值，迭代序列会收敛于不同的根，所以在某个区间内求根对于初值的选取有很大的关系。产生上述结果的原因是区间不满足大范围收敛的条件。习题335（上机题）列主元三角分解法对于某电路的分析，归结为求解线性方程组RI=V。（1）编制解n阶线性方程组Ax=b的列主元三角分解法的通用程序；（2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留

8、五位有效数；（3）本编程之中，你提高了哪些编程能力？程序为：#include#includevoid main(void)int i,j,n,k,q;float a1011,s10,s110;coutn;cout输入数组a:endl;for(i=1;i=n;i+)for(j=1;jaij; /给矩阵a赋值for(i=1;i=n;i+)for(j=1;j=(n+1);j+)coutaijt;coutn; /输出数组acoutn;/进行第一行和第一列元素的求取/int t=1;for(i=1;i=n;i+)si=ai1;float max=fabs(s1);for(i=2;imax)max=fab

9、s(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=(n+1);j+)coutaijt;coutn;/输出进行第一步变换的数组acoutn;/进行第k步分解/for(k=2;k=n;k+)for(i=k;i=n;i+)float sum=0;for(q=1;qk;q+)sum+=aiq*aqk;s1i=aik-sum;int l=k;float m=fabs(s1k);for(i=k;im)m

10、=fabs(s1i);l=i; /返回行值for(j=1;j=n+1;j+) /交换两行元素float s2=akj;akj=alj;alj=s2;for(j=k;j=n+1;j+)/算出第k行行元素的值float sum1=0;for(q=1;qk;q+)sum1+=akq*aqj;akj=akj-sum1;for(i=k+1;i=n;i+)/算出第k列列元素的值float sum2=0;for(q=1;qk;q+)sum2+=aiq*aqk;aik=(aik-sum2)/(akk); /第k步分解结束for(i=1;i=n;i+)for(j=1;j=(n+1);j+)coutaijt;co

11、ut=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+)coutxi=xiendl; /输出解向量结果：方程的解为：x1= -0.28923，x2= 0.34544，x3= -0.71281，x4= -0.22061，x5= -0.43040，x6= 0.15431，x7= -0.057823，x8= 0.20105，x9= 0.29023。分析：我感觉是提高了查错误点的能力和编循环语句的能力，即利用很规整的迭代公式进行编程。另外列主元三角分

12、解法的阶梯步骤有了更深的了解!36逐次超松弛迭代法（1）编制解n阶线性方程组Ax=b的SOR方法的通用程序（要求）；（2）对于35题中所给的线性方程组，取松弛因子，容许误差，打印松弛因子、迭代次数、最佳松弛因子及解向量。程序为：#include#include#define eps 0.5e-5 /迭代误差void main(void)int i,j,l;float w,t;float m9;float sum;float a99=31,-13,0,0,0,-10,0,0,0,-13,35,-9,0,-11,0,0,0,0,0,-9,31,-10,0,0,0,0,0,0,0,-10,79,-3

13、0,0,0,0,-9,0,0,0,-30,57,-7,0,-5,0,0,0,0,0,-7,47,-30,0,0,0,0,0,0,0,-30,41,0,0,0,0,0,0,-5,0,0,27,-2,0,0,0,-9,0,0,0,-2,29;float b9=-15,27,-23,0,-20,12,-7,7,10;float max(float m9);for(t=1;t=99;t+)l=0;float x09=1,1,1,1,1,1,1,1,1;float x19=1,1,1,1,1,1,1,1,1;w=t/50;dofor(i=0;i9;i+)x0i=x1i;for(i=0;i9;i+)sum

14、=0;for(j=0;ji;j+)sum=sum+aij*x1j;for(j=i+1;j9;j+)sum=sum+aij*x0j;x1i=(1-w)*x0i+w*(bi-sum)/aii; /解出九个解for(i=0;i=eps);if(max(m)=eps)cout迭代次数=ltw=wn;for(i=0;i9;i+)coutx1i=x1it;cout-n;float max(float m9) /求出最大的迭代误差float k;k=(fabs(m0);for(int i=1;ik)k=fabs(mi);return k;结果为： 迭代次数 迭代次数 迭代次数 迭代次数0.02 12910.

15、04 7000.06 4860.08 3730.10 3030.12 2550.14 2210.16 1940.18 1730.20 1560.22 1420.24 1300.26 1200.28 1110.30 1030.32 960.34 900.36 850.38 800.40 750.42 710.44 680.46 640.48 610.50 580.52 550.54 530.56 510.58 480.60 460.62 440.64 420.66 410.68 390.70 370.72 360.74 340.76 330.78 320.80 300.82 290.84 280

16、.86 270.88 260.90 250.92 240.94 230.96 220.98 211.00 201.02 191.04 181.06 171.08 161.10 151.12 151.14 141.16 121.18 101.20 111.22 121.24 121.26 131.28 131.30 141.32 151.34 151.36 161.38 171.40 181.42 191.44 191.46 201.48 211.50 221.52 241.54 251.56 271.58 301.60 311.62 341.64 361.66 391.68 441.70 49

17、1.72 551.74 601.76 711.78 811.80 971.82 1211.84 1571.86 2281.88 3981.90 15821.92 54271.94 21781.96 13741.98 1009 从1.92到2.00均出现不合理的迭代次数，迭代次数偏大。不是合理的迭代值。当初值为x=（1,1,1,1,1,1,1,1,1,）T 时,从上表可以看出,最佳松弛因子为=1.18,迭代次数仅为10次,方程的解为: x1= -0.289231，x2= 0.345437，x3= -0.712811， x4= -0.220608，x5= -0.430400，x6= 0.15430

18、9，x7= -0.057823，x8= 0.201054，x9= 0.290229。数 值 分 析 上 机 题姓名：戴载星 学号：040139习题 438.（上机题）3次样条插值函数（1）编制求第一型3次样条插值函数的通用程序; （2） 已知汽车曲线型值点的数据如下：0123456789102.513.304.044.705.225.545.785.405.575.705.80端点条件为=0.8，=0.2。用所编制程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,9)。解：通用程序：第 8 页 共 16 页#includevoid main(void)float x1

19、1;/存放数组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为区间值/将初值初始化coutn;cout输入数组x:n;for(i=0;ixi;cout输入数组y:n;for(i=0;iyi;/输入端点值float df2;cout输入两个端点值:n;for(i=

20、0;idfi;/利用书本上的算法求出所需要的值/求出hj的值for(j=0;j=n-1;j+)hj=xj+1-xj;couthj=hjt;coutendl;/求出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;jn;j+)dj=6*(yj+1-yj)/hj-(yj-yj-1)/hj-1)/(hj+hj-1);d0=6*(y1-y0)/h0-df0)/h0;dn=6*(df1-(yn-yn-1)/hn-1)/hn-1;for(j=1;j=n;j+)coutuj=ujt;cout

21、endl;for(j=0;jn;j+)coutvj=vjt;coutendl;for(j=0;j=n;j+)coutdj=djt;coutendl;/利用书本上的追赶法求解方程组for(i=0;i=n;i+)bi=2;coutendl;t0=b0;yy0=d0;/消元过程for(i=1;i=0;i-)Mi=(yyi-vi*Mi+1)/ti;/将Mj的值输出for(i=0;i=n;i+)coutMi=Miendl;/输出插值多项式的系数for(j=0;jn;j+)s0=yj;s1=(yj+1-yj)/hj-(Mj/3+Mj+1/6)*hj;s2=Mj/2;s3=(Mj+1-Mj)/(6*hj);

22、cout当x的值在区间xj到x(j+1)时,输出插值多项式的系数:n;for(int k=0;k4;k+)coutsk=skt;coutendl;（2）编制的程序求车门的3次样条插值函数S(x)：x属于区间0,1时；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(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)x属于区间2,3时；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543

23、(x-2)(x-2)(x-2)x属于区间3,4时；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3)x属于区间4,5时；S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4)x属于区间5,6时；S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5)x属于区间6,7时；S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.4

24、26588(x-6)(x-6)(x-6)x属于区间7,8时；S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7)x属于区间8,9时；S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8)x属于区间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(0.5)=2.90856 S(1.5)=3.67843 S (2.5)=4.38147S(3.5

25、)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976S(9.5)=5.7323习题五 重积分的计算23（上机题）重积分的计算题目：给定积分。取初始步长h和k，及精度。应用复化Simpson公式，采用逐次二分步长的方法，编制计算I(f)的通用程序。计算至相邻两次近似值之差的绝对值不超过为止。1） 用所编程序计算积分，取。算法概述初始时候只在x，y方向上各二分一次，根据复化Simpson公式计算积分值，然后再二分一次，仍然根据上式重新计算积分值，比较两次计算结果的差值，如果小于误差限，则

26、已求得满足要求的结果，否则继续二分区间直到满足误差要求为止。程序如下：#include#include#define PI 3.1415926#define error 0.5E-5double f(double x,double y) double answ; answ=tan(x*x+y*y); return answ;void main() double old,temp,a=0,c=0; int i,j,m=1,n=1; double b=PI/3;/重积分内层上限 double d=PI/6;/重积分外层上限 double h=(b-a)/(2*n);/ double k=(d-c)

27、/(2*m); double answ=0; do old=answ; answ=0; answ+=f(a,c); temp=0; for(i=1;in;i+) temp+=f(a+2*i*h,c); temp*=2; answ+=temp; temp=0; for(i=1;i=n;i+) temp+=f(a+(2*i-1)*h,c); temp*=4; answ+=temp; answ+=f(b,c); temp=0; for(j=1;jm;j+) temp+=f(a,c+2*j*k); temp*=2; answ+=temp; temp=0; for(j=1;jm;j+) for(i=1;

28、in;i+)temp+=f(a+2*i*h,c+2*j*k); temp*=4; answ+=temp; temp=0; for(j=1;jm;j+) for(i=1;i=n;i+)temp+=f(a+(2*i-1)*h,c+2*j*k); temp*=8; answ+=temp; temp=0; for(j=1;jm;j+) temp+=f(b,c+2*j*k); temp*=2; answ+=temp; temp=0; for(j=1;j=m;j+) temp+=f(a,c+(2*j-1)*k); temp*=4; answ+=temp; temp=0; for(j=1;j=m;j+) f

29、or(i=1;in;i+)temp+=f(a+2*i*h,c+(2*j-1)*k); temp*=8; answ+=temp; temp=0; for(j=1;j=m;j+) for(i=1;i=n;i+)temp+=f(a+(2*i-1)*h,c+(2*j-1)*k); temp*=16; answ+=temp; temp=0; for(j=1;j=m;j+) temp+=f(b,c+(2*j-1)*k); temp*=4; answ+=temp; answ+=f(a,d); temp=0; for(i=1;in;i+) temp+=f(a+2*i*h,d); temp*=2; answ+=

30、temp; temp=0; for(i=1;ierror); coutansw is:answ; coutnDivided intom/2partsn; return;程序的输出结果为：Result is:0.33652Divided into64parts即在x，y方向上各二分6次。本题公式较长，形式较为复杂，但是并不需要太多编程技巧，只需细心即可。随着二分的继续，计算的结果越来越趋向准确值。本题在二分6次后，经检验误差符合精度要求。习题621（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正

31、方法的通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;（5）对于初值问题取步长,应用（1）（4）中的四种方法进行计算，并将计算结果和精确解作比较；（6）通过本上机题，你能得到哪些结论？解：程序为：第 16 页 共 16 页#include#include#include#includeofstream outfile(data.txt);/此处定义函数f(x,y)的表达式/用户可以自己设定所需要求得函数表达式double f1(double x,double y)double f1;f1=(-1)*x*x*y*y;return f1;/

32、此处定义求函数精确解的函数表达式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);outfileaccuratei=accuratein;i+;while(i(b-a)/h+1);/此处为经典Runge-Kutta公式的通用程序void RK4(doub

33、le 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;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;outfileyi=yin;i+;while(i(b-a)/h+1);/此处为4阶Adams显式方法的通用程序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