数值分析上机作业最强版.docx

《数值分析上机作业最强版.docx》由会员分享，可在线阅读，更多相关《数值分析上机作业最强版.docx（25页珍藏版）》请在三一办公上搜索。

1、精选优质文档-倾情为你奉上 数值分析上机作业姓 名：唐 皓学 号：专 业：道路与铁道工程院 系： 交通学院授课教师：吴 宏 伟日 期：2015年1月习题一1 题目17（上机题）舍入误差与有效数设，其精确值为。（1）编制按从大到小的顺序，计算的通用程序；（2）编制按从小到大的顺序，计算的通用程序；（3）按两种顺序分别计算，并指出有效位数。（编制程序时用单精度）；（4）通过本上机题你明白了什么？2 通用程序代码2.1 按从小到大的顺序计算void AscendSum(unsigned long int N)/ 计算从大到小的总和for(unsigned long int j=2;j=N;j+)as

2、cendSum+=(float)1.0/(j*j-1);coutSum From 1 to N (Ascend) is: ascendSum=2;j-)descendSum+=(float)1.0/(j*j-1);coutSum From N to 1 (Descend) is: descendSumendl;Error(descendSum); Delimiter();3 计算结果展示图1 N=100时的计算结果图2 N=10000时的计算结果图3 N=时的计算结果专心-专注-专业表1-1 计算结果汇总精确值按从小到大按从大到小值有效位数值有效位数0.0.100.60.0.40.100.0.

3、20.24 计算结果分析（1）如果采用单精度数据结构进行计算，则相较于双精度的数据结果，由于数据存储字长的限制导致计算机存在较大的舍入误差，因此本程序采用的是双精度数据存储方式。（2）由计算结果可知，正序计算和逆序计算的精度是不稳定的。由计算结果可以发现，当N=100时，正序计算（1-N）的精度较高；当N=10000时，逆序计算（N-1）的精度较高；当N=时，正序计算和逆序计算的精度一样。当然，和其他同学做出来的结果对比，结论并不一致。我个人分析这是因为电脑的硬件、软件（位数）不同等原因导致的。但总体而言，在N较小时，正序计算精度高于逆序计算的精度。当N较大时，正序和逆序计算的精度接近。（3）

4、由于计算机的实际计算过程是一种舍入机制，故对于我们计算所采用的加法交换律是不成立的。计算机中若干数相加时，先要进行对阶操作，即将两数的阶数统一为绝对值较大的数的阶数。这样一来将导致绝对值较小的数的有效数字可能会大量损失，增大舍入误差，即所谓的“大数吃小数”现象。为了避免这种现象的出现，在进行加减法的时候应该先将绝对值较小的数相加，再与绝对值较大的数相加这样按阶逐步递增的相加。5 完整代码#include #include #include using namespace std;float accurateSum=0,ascendSum=0,descendSum=0;void Delimite

5、r()/输出一系列星号以间隔for(int i=1;i=50;i+) cout*;coutendl;void Error(float Sum)/计算绝对误差float error;error=fabs(Sum-accurateSum);int flag;for(flag=0;flag0.5)break;coutThere are flag Valid numbers.n;void AccurateSum(unsigned long int N)/计算精确值accurateSum=0.5*(1.5-(float)1/N-(float)1/(N+1);coutAccurate sum is: se

6、tprecision(10)accurateSumendl;Delimiter();void AscendSum(unsigned long int N)/计算从大到小的总和for(unsigned long int j=2;j=N;j+)ascendSum+=(float)1.0/(j*j-1);coutSum From 1 to N (Ascend) is: ascendSum=2;j-)descendSum+=(float)1.0/(j*j-1);coutSum From N to 1 (Descend) is: descendSumendl;Error(descendSum); Del

7、imiter();void main()/主程序long int N;while(1)cout=2):;cinN;if (N2)cout=eps)x1=x2;x2=x1-f(x1)/df(x1);return x1;2.2 求解尽可能大double MaximalDeviateRange() /求解尽可能大的范围double step=1e-5; /step lengthint cnt=1; /step count double delta;cout*Newton Iteration (eps=1e*5)*endl;while(fabs(NewtonIteration(step*cnt,eps

8、)=eps) cnt+;delta=step*cnt;coutThe maximal deviate range for x converging to x2*=0 is (-delta,delta)endl;return delta;3 计算结果展示图2-4 计算结果在取步长为10-5的情况下，允许误差eps=10-5时有轴上的一个小区间 (-0.7746,0.7746)为Newton迭代序列在x2*=0处的尽可能大的局部收敛区间。当x0=(-,1)，(-1,-)，(-,)，(,1)，(1,)时牛顿迭代序列分别收敛于-1.,1.,0,-1.,1.。4 计算结果分析（1）通过本次上机编程并通过

9、多次的调试，可以发现运行结果很好的验证了教材上牛顿迭代法具有局部收敛性这一重要性质。（2）选择不同的初值区间，迭代序列会收敛于不同的根。所以为了收敛到需要的计算结果，就需要选择合适的牛顿迭代法大范围收敛区间。（3）产生上述结果的原因是所选取的部分区间不满足大范围收敛的条件，导致并没有收敛到理想的结果。5 完整代码#include#include#include using namespace std;double eps=1e-5;double f(double x)/函数f(x)的录入return (x*x*x)/3-x;double df(double x)/函数f(x)的导数return

10、 (x*x)-1;double NewtonIteration(double x0,double eps)/Newton迭代法求解子程序double x1,x2;x1=x0;x2=x1-f(x1)/df(x1);while (fabs(x1-x2)=eps)x1=x2;x2=x1-f(x1)/df(x1);return x1;double MaximalDeviateRange() /求解尽可能大的范围double step=1e-5; /step lengthint cnt=1; /step count double delta;cout*Newton Iteration (eps=1e*5

11、)*endl;while(fabs(NewtonIteration(step*cnt,eps)=eps) cnt+;delta=step*cnt;coutThe maximal deviate range for x converging to x2*=0 is (-delta,delta)endl;return delta;void Calculate(double x0)/计算Newton迭代法和相应的子程序printf(If x0=% 11.4f, x converges to the value of: % 8.6fn,x0,NewtonIteration(x0,eps);void m

12、ain()/主程序double delta=MaximalDeviateRange(); Calculate(-10000);Calculate(-1-delta)/2);Calculate(-delta/2);Calculate(delta/2);Calculate(1+delta)/2);Calculate(10000);习题三1 题目39（上机题）列主元Guass消去法对于某电路的分析，归结为求解线性方程组RI=V。其中（1）编制解阶线性方程组Ax=b的列主元Guass消去法的通用程序；（2）用所编程序解线性方程组RI=V，并打印出解向量，保留5位有效数；（3）本题编程之中，你提高了哪些

13、编程能力？2 通用程序代码void GaussElimination(float matrixROWCOL)/高斯列主元消去法主程序int row,col;int tmpRow,maxRow,tmpCol;float ratio;float rootROW;/Calculate the echelon matrixfor(row=0;rowROW-1;row+)col=row;maxRow=row;for(tmpRow=row;tmpRowfabs(matrixmaxRowcol)maxRow=tmpRow;SwapRow(matrixrow,matrixmaxRow);for(tmpRow=

14、row+1;tmpRowROW;tmpRow+)ratio=matrixtmpRowcol/matrixrowcol;if(ratio=0)continue;for(tmpCol=0;tmpCol=0;tmpRow-)sum=0;for(cnt=tmpRow+1;cntROW;cnt+)sum+=AugmentedMatrixtmpRowcnt*rootcnt;roottmpRow=(AugmentedMatrixtmpRowCOL-1-sum)/AugmentedMatrixtmpRowtmpRow;printf(The Solution Vector is:n);for(cnt=0;cnt

15、ROW;cnt+)printf(% 8.5f ,rootcnt);printf(n);3 计算结果展示图 3-5计算结果解向量为：4 结果分析和心得体会（1）在本题编程过程中，由于考虑到列主元Gauss消去法的算法特点：程序的过程重复性高，故本程序采用模块化程序设计方法。这样就可以将每步都要重复进行的消去、交换列主元等过程的功能以模块函数来实现。这样可以使程序结构清晰，从而调试变得相对容易。（2）在输入/输出的控制方面，本程序采用格式化输出，使输出结果自动对齐，提高了程序的结果的可阅读性。（3）为了能更方便快捷地进行矩阵数据的输入。将需要求解的线性方程增广矩阵存储在文件中，程序运行时自动调用并

16、输出。因此，对于任意线性矩阵，只需将其增广矩阵存储于其中，并修改ROW和COL为实际行数与列数即可。（4）通过该程序的编写，提高了我对循环语句的运用、子程序调用以及找bug（错位）的能力。此外，我也对列主元Gauss消去法的步骤有了更加深刻的理解。5 完整代码#include#include#include#includeconst int ROW=9,COL=10;float AugmentedMatrixROWCOL;void PrintMatrix(float matrixROWCOL)/输出解子程序int row,col;for(row=0;rowROW;row+)for(col=0;

17、colCOL;col+)printf(% 9.5f ,matrixrowcol);printf(n);void PrintDelimiter()for(int i=0;i100;i+)printf(*);printf(n);void ReadData(char filename, float matrixROWCOL)/读入数据程序ifstream infile(filename);if (!infile)coutCannot open the Augmented Matrix data file!n;coutThe Augmented Matrix is:endl;for(int i=0;i

18、ROW;i+)for(int j=0;jmatrixij;coutsetw(3)matrixij ;coutendl;PrintDelimiter();infile.close();void SwapRow(float rowA, float rowB)/交换两行的子程序int i;float tmp;for(i=0;iCOL;i+)tmp=rowAi;rowAi=rowBi;rowBi=tmp;void GaussElimination(float matrixROWCOL)/高斯列主元消去法主程序int row,col;int tmpRow,maxRow,tmpCol;float rati

19、o;float rootROW;/Calculate the echelon matrixfor(row=0;rowROW-1;row+)col=row;maxRow=row;for(tmpRow=row;tmpRowfabs(matrixmaxRowcol)maxRow=tmpRow;SwapRow(matrixrow,matrixmaxRow);for(tmpRow=row+1;tmpRowROW;tmpRow+)ratio=matrixtmpRowcol/matrixrowcol;if(ratio=0)continue;for(tmpCol=0;tmpCol=0;tmpRow-)sum=

20、0;for(cnt=tmpRow+1;cntROW;cnt+)sum+=AugmentedMatrixtmpRowcnt*rootcnt;roottmpRow=(AugmentedMatrixtmpRowCOL-1-sum)/AugmentedMatrixtmpRowtmpRow;printf(The Solution Vector is:n);for(cnt=0;cntROW;cnt+)printf(% 8.5f ,rootcnt);printf(n);void main()/主程序ReadData(AugmentedMatrix.txt,AugmentedMatrix);GaussElim

21、ination(AugmentedMatrix);题目四1 题目37（上机题）3次样条插值函数（1）编制求第一型3次样条插值函数的通用程序；（2）已知汽车门曲线型值点的数据如下：0123456789100123456789102.513.304.044.705.225.545.785.405.575.705.80端点条件为y0=0.8，y10=0.2，用所编制程序求车门的3次样条插值函数S(x),并打印出Si+0.5，i=0,1,9。2 计算结果展示通过运行程序，求得S(X)的表达式后，分别求各个点的值，得到如下结果：图 4-6 运行结果即有：S(0.5)= 2.S(1.5)= 3.S(2.5

22、)= 4.S(3.5)= 4.S(4.5)= 5.S(5.5)= 5.S(6.5)= 5.S(7.5)= 5.S(8.5)= 5.S(9.5)= 5.3 完整代码#include#include#includeusing namespace std;const n=10;class SPLINE/样条类定义private:double dy0,dyn;double mn+1,hn; void INSERT(); void CHASE(double an+1,double bn+1,double cn,double dn+1);public:double xn+1,fn+1;void INIT(

23、double x_1n+1,double f_1n+1,double dy0_1,double dyn_1); double SX(double xk); spline;void SPLINE:INIT(double x_1n+1,double f_1n+1,double dy0_1,double dyn_1)for(int i=0;i=n;i+)xi=x_1i;fi=f_1i;dy0=dy0_1; dyn=dyn_1;INSERT();void SPLINE:CHASE(double an+1,double bn+1,double cn,double dn+1)for(int i=1;i=0

24、;i-) mi=(di-ci*mi+1)/bi;void SPLINE:INSERT() double yn+1;for(int i=0;i=1;i-) yi=(yi-yi-1)/(xi+1-xi-1);y0=(y0-dy0)/(x1-x0);for(i=0;i=n;i+) yi*=6;double un+1,ln;for(i=0;i=n-1;i+) hi=xi+1-xi; for(i=1;i=n-1;i+) ui=hi-1/(hi-1+hi); li=1-ui; l0=1; un=1; double bn+1; for(i=0;i=xn|xk=x0) coutError!=xi+1) i+;

25、double s=0;s=mi*pow(xi+1-xk),3)/(6*hi) + mi+1*pow(xk-xi),3)/(6*hi)+(mi*hi*hi/6-fi)*(xk-xi+1)/hi + (fi+1-mi+1*hi*hi/6)*(xk-xi)/hi;return s;void Init()double xn+1=0,1,2,3,4,5,6,7,8,9,10, fn+1=2.51,3.3,4.04,4.7,5.22,5.54, 5.78,5.4,5.57,5.7,5.8;double dy0=0.8,dyn=0.2; spline.INIT(x,f,dy0,dyn);cout= Trig

26、onometric Interpolation=endl;void main()/主程序 Init(); for(int i=0;i=9;i+) coutS(i+0.5)= setprecision(9)spline.SX(double)(i+0.5)endl; 习题六1 题目23（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（5）对于初值问题取步长h=0.1,应用（1）（4）中的

27、四种方法进行计算，并将计算结果和精确解yx=3/(1+x3)作比较；（6）通过本上机题，你能得到哪些结论？2 通用程序2.1 编制RK4方法的通用程序；void AB4(int j)/AB4方法的通用程序yj+1=yj+float(h/24)*(55*f(xj,yj)-59*f(xj-1,yj-1)+37*f(xj-2,yj-2)-9*f(xj-3,yj-3);2.2 编制AB4方法的通用程序（由RK4提供初值）；void AB4(int j)/AB4方法的通用程序yj+1=yj+float(h/24)*(55*f(xj,yj)-59*f(xj-1,yj-1)+37*f(xj-2,yj-2)-

28、9*f(xj-3,yj-3);2.3 编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；void AM4(int j)/AB4-AM4预测校正方法的通用程序yj+1=yj+float(h/24)*(9*f(xj+1,yj+1)+19*f(xj,yj)-5*f(xj-1,yj-1)+f(xj-2,yj-2);2.4 编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）；i=3;/带改进的AB4-AM4法带入计算（包括带改进的通用程序）while(xi1.5) AB4(i);AM4(i);i+;coutAB4-AM4ty1.5=y14ty(1.5)=AY(15)tAY(

29、15)-y14endl;3 计算结果展示图6-7 运行结果RK4y1.5=0.63238y(1.5)=0. 0.AB4y1.5=0.y(1.5)=0. 0.AM4y1.5=0.y(1.5)=0. 0.AB4-AM4y1.5=0.y(1.5)=0. 0.4 结论分析（1）采用单精度计算，可以发现计算结果误差都比较大，也比较接近。（2）若将计算数据结构改为双精度，可以发现求微分方程的数值方法的精度是不同的。如在本题中，经典的Runge-Kutta公式，AB4方法以及AB4-AM4预测校正方法求解公式的精度是不同的。若将所有每步计算结果做为向量求范数，可以发现在本案例中单步方法RK4的精度很好，多步

30、法中改进AB4-AM4得精度比AB4-AM4高，AB4的精度是最底的。5 完整代码#include #include using namespace std;float h=0.1;const int n=15;float xn,yn;float f(float x,float yy)float temp=0;temp=-1*x*x*yy*yy;return temp;void RK4(int j)/RK4方法的通用程序float k1,k2,k3,k4;k1=f(xj,yj);k2=f(xj+0.5*h,yj+0.5*h*k1);k3=f(xj+0.5*h,yj+0.5*h*k2);k4=f

31、(xj+h,yj+h*k3);yj+1=yj+(float)(h/6)*(k1+2*k2+2*k3+k4);void AB4(int j)/AB4方法的通用程序yj+1=yj+float(h/24)*(55*f(xj,yj)-59*f(xj-1,yj-1)+37*f(xj-2,yj-2)-9*f(xj-3,yj-3);void AM4(int j)/AB4-AM4预测校正方法的通用程序yj+1=yj+float(h/24)*(9*f(xj+1,yj+1)+19*f(xj,yj)-5*f(xj-1,yj-1)+f(xj-2,yj-2);float AY(int j) /精确值计算return(3

32、/(1+xj*xj*xj);void Outfit()int i;int k;for(k=0;kn;k+)xk=k*h; yk=0;y0=3;x15=1.5;RK4(0);RK4(1);RK4(2);RK4(3);i=3;/以RK4法带入计算while(xi1.5) RK4(i); i+;coutRK4ty1.5=y14ty(1.5)=AY(15)tAY(15)-y14endl;i=3;/以AB法带入计算while(xi1.5) AB4(i);i+;coutAB4ty1.5=y14ty(1.5)=AY(15)tAY(15)-y14endl;i=3;/以AM4法带入计算while(xi1.5) AM4(i); i+;coutAM4ty1.5=y14ty(1.5)=AY(15)tAY(15)-y14endl;i=3;/带改进的AB4-AM4法带入计算（包括带改进的通用程序）while(xi1.5) AB4(i);AM4(i);i+;coutAB4-AM4ty1.5=y14ty(1.5)=AY(15)tAY(15)-y14endl;void main()/主程序Outfit();