合工大计算方法实验报告

1、计算方法试验报告班级：学号：姓名：实验一、牛顿下山法1 实验目的(1) 熟悉非线性方程求根简单迭代法，牛顿迭代及牛顿下山法(2) 能编程实现简单迭代法，牛顿迭代及牛顿下山法(3) 认识选择迭代格式的重要性(4) 对迭代速度建立感性的认识；分析实验结果体会初值对迭代的影响2 实验内容（1）用牛顿下山法解方程（初值为0.6）输入：初值，误差限，迭代最大次数，下山最大次数输出：近似根各步下山因子 （2）设方程f(x)=x- 3x 1=0 有三个实根 x=1.8793 , x=-0.34727 ,x=-1.53209现采用下面六种不同计算格式，求 f(x)=0的根 x 或x x = ； x = ；x

2、= ；x = ；x = ；x = x - 输入：初值，误差限，迭代最大次数输出：近似根、实际迭代次数3 算法基本原理求非线性方程组的解是科学计算常遇到的问题，有很多实际背景各种算法层出不穷，其中迭代是主流算法。只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。因此设计算法之前，对于一般迭代进行收敛性的判断是至关重要的。牛顿法也叫切线法，是迭代算法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代至关重要。当初值选取不当可以采用牛顿下山算法进行纠正。一般迭代： 牛顿公式：牛顿下山公式： 图3.1一般迭代算法流程图下山因子下山条件4 算法描述一般迭代算法见流程图图3

3、.2牛顿下山算法流程图牛顿下山算法见流程图：5、代码： #include <iostream> #include <fstream> #include <cmath> using namespace std; class srrt private: int n; double *a, *xr, *xi; public: srrt (int nn) n = nn; a = new doublen+1; /动态分配内存 xr = new doublen; xi = new doublen; void input (); /由文件读入代数方程系数 void srr

4、t_root (); /执行牛顿下山法 void output (); /输出根到文件并显示 srrt () delete a, xr, xi; ; void srrt:input () /由文件读入代数方程系数 int i; char str120; cout <<"n输入文件名: " cin >>str1; ifstream fin (str1); if (!fin) cout <<"n不能打开这个文件 " <<str1 <<endl; exit(1); for (i=n; i>=0;

5、 i-) fin >>ai; /读入代数方程系数 fin.close (); void srrt:srrt_root () /执行牛顿下山法 int m,i,jt,k,is,it; double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c; double g,u,v,pq,g1,u1,v1; m=n; while (m>0)&&(fabs(am)+1.0=1.0) m=m-1; if (m<=0) cout <<"n程序工作失败！" <<endl; return; for (i=0; i&l

6、t;=m; i+) ai=ai/am; for (i=0; i<=m/2; i+) w=ai; ai=am-i; am-i=w; k=m; is=0; w=1.0; jt=1; while (jt=1) pq=fabs(ak); while (pq<1.0e-12) xrk-1=0.0; xik-1=0.0; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; return; pq=fabs(ak); q=log(pq); q=q/(1.0*k); q=exp(q); p=q; w=w*p; for (i=1; i<=k; i+) ai=ai/q;

7、q=q*p; x=0.0001; x1=x; y=0.2; y1=y; dx=1.0; g=1.0e+37; l40: u=a0; v=0.0; for (i=1; i<=k; i+) p=u*x1; q=v*y1; pq=(u+v)*(x1+y1); u=p-q+ai; v=pq-p-q; g1=u*u+v*v; if (g1>=g) if (is!=0) it=1; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (dd>1.0) dd=1.0; dc=6.28/(4.5*k); c=0.0; while(1=1) c=c+dc; dx=

8、dd*cos(c); dy=dd*sin(c); x1=x+dx; y1=y+dy; if (c<=6.29) it=0; break; dd=dd/1.67; if (dd<=1.0e-07) it=1; break; c=0.0; if (it=0) goto l40; else it=1; while (it=1) t=t/1.67; it=0; x1=x-t*dx; y1=y-t*dy; if (k>=50) p=sqrt(x1*x1+y1*y1); q=exp(85.0/k); if (p>=q) it=1; if (t>=1.0e-03) goto l

9、40; if (g>1.0e-18) it=0; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (dd>1.0) dd=1.0; dc=6.28/(4.5*k); c=0.0; while(1=1) c=c+dc; dx=dd*cos(c); dy=dd*sin(c); x1=x+dx; y1=y+dy; if (c<=6.29) it=0; break; dd=dd/1.67; if (dd<=1.0e-07) it=1; break; c=0.0; if (it=0) goto l40; if (fabs(y)<=1.0e-

10、06) p=-x; y=0.0; q=0.0; else p=-2.0*x; q=x*x+y*y; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i<=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; else g=g1; x=x1; y=y1; is=0; if (g<=1.0e-22) if (fabs(y)<=1.0e-06) p=-x; y=0.0; q=0.0; else p=-2.0

11、*x; q=x*x+y*y; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i<=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; else u1=k*a0; v1=0.0; for (i=2; i<=k; i+) p=u1*x; q=v1*y; pq=(u1+v1)*(x+y); u1=p-q+(k-i+1)*ai-1; v1=pq-p-q; p=u1*u1+v1*v1; if (p<=1.0

12、e-20) it=0; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (dd>1.0) dd=1.0; dc=6.28/(4.5*k); c=0.0; while(1=1) c=c+dc; dx=dd*cos(c); dy=dd*sin(c); x1=x+dx; y1=y+dy; if (c<=6.29) it=0; break; dd=dd/1.67; if (dd<=1.0e-07) it=1; break; c=0.0; if (it=0) goto l40; if (fabs(y)<=1.0e-06) p=-x; y=0.0;

13、 q=0.0; else p=-2.0*x; q=x*x+y*y; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i<=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; else dx=(u*u1+v*v1)/p; dy=(u1*v-v1*u)/p; t=1.0+4.0/k; it=1; while (it=1) t=t/1.67; it=0; x1=x-t*dx; y1=y-t*dy; if (k>

14、=50) p=sqrt(x1*x1+y1*y1); q=exp(85.0/k); if (p>=q) it=1; if (t>=1.0e-03) goto l40; if (g>1.0e-18) it=0; if (it=0) is=1; dd=sqrt(dx*dx+dy*dy); if (dd>1.0) dd=1.0; dc=6.28/(4.5*k); c=0.0; while(1=1) c=c+dc; dx=dd*cos(c); dy=dd*sin(c); x1=x+dx; y1=y+dy; if (c<=6.29) it=0; break; dd=dd/1.

15、67; if (dd<=1.0e-07) it=1; break; c=0.0; if (it=0) goto l40; if (fabs(y)<=1.0e-06) p=-x; y=0.0; q=0.0; else p=-2.0*x; q=x*x+y*y; xrk-1=x*w; xik-1=-y*w; k=k-1; for (i=1; i<=k; i+) ai=ai-ai-1*p; ai+1=ai+1-ai-1*q; xrk-1=x*w; xik-1=y*w; k=k-1; if (k=1) xr0=-a1*w/a0; xi0=0.0; if (k=1) jt=0; else

16、 jt=1; void srrt:output () /输出根到文件并显示 int k; char str220; cout <<"n输出文件名: " cin >>str2; ofstream fout (str2); if (!fout) cout <<"n不能打开这个文件 " <<str2 <<endl; exit(1); for (k=0; k<n; k+) fout <<xrk <<" " <<xik <<endl

17、; cout <<xrk <<" +j " <<xik <<endl; fout.close (); void main () /主函数 srrt root(6); root.input (); /由文件读入代数方程系数 root.srrt_root (); /执行牛顿下山法 root.output (); /输出根到文件并显示 6、输入输出输出结果如下：7、分析体会牛顿下山法作为计算方法课程中重要的知识点，在看书学习时较易大致理解其思路，但上级编写代码时却是有难度的。在学习借鉴后用如上代码能够实现基本功能，在实验过程中也加深

18、了我对牛顿下山法的理解。实验二 高斯消去法1 实验目的(1) 熟悉求解线性方程组的有关理论和方法；(2) 能编程实现雅可比及高斯-塞德尔迭代法、列主元高斯消去法、约当消去，追赶法(3) 通过测试，进一步了解各种方法的优缺点(4) 根据不同类型的方程组，选择合适的数值方法2 实验内容（1）用Gauss - Seidel 迭代法求解方程组输入：系数矩阵A，最大迭代次数N，初始向量，误差限e输出：解向量（2）一个城镇有三个主要生产企业：煤矿、电厂和地方铁路作为它的经济系统.煤矿：生产价值1元的煤，需消耗0.25元的电费和0.35元的运输费；电厂：生产价值1元的电，需消耗0.40元的煤费、0.05元的

19、电费和0.10元的运输费；铁路：而提供价值1元的铁路运输服务，则需消耗0.45元的煤、0.10元的电费和0.10元的运输费.在某个星期内，除了这三个企业间的彼此需求，煤矿得到50000元的订单，电厂得到25000元的电量供应要求，而地方铁路得到价值30000元的运输需求.提问：这三个企业在这星期各应生产多少产值才能满足内外需求？提示：“投入”和“产出”之间的平衡关系：物资消耗和新创造的价值等于它生产的总产值. 输入：物资消耗和新创造的价值输出：三个企业生产总值（3）一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立一个以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（1天文单位

20、为地球到太阳的平均距离：9300万里）。他五个不同时间对小行星作五次观测，得到轨道上五个点的坐标分别为（5.764,0.648）（6.286,1.202）（6.759,1.823）（7.168,2.562）与（7.408，3.360）。由开普勒第一定律知小行星轨道为一椭圆，试建立它的方程。（4）用选主元高斯消去求行列式值提示：A.B. 消元结果直接存储在系数矩阵中C. 当消元过程发生两行对调的情况为偶数次时，行列式值为对角线乘积，否则为对角线乘积的相反数（5）用选主元约当消去分别对矩阵求其逆矩阵，若不可逆输出奇异标志和提示：3 算法基本原理无论是三次样条还是拟合问题最终都归结为线性方程组，求解

21、线性方程组在数值分析中非常重要，在工程计算中也不容忽视。线性方程组大致分迭代法和直接法。只有收敛条件满足时，才可以进行迭代。雅可比及高斯-塞德尔是最基本的两类迭代方法，最大区别是迭代过程中是否引用新值进行剩下的计算。消元是最简单的直接法，并且也十分有效的，列主元高斯消去法对求解一般的线性方程组都适用，同时可以用来求矩阵对应的行列式。约当消去实质是经过初等行变换将系数矩阵化为单位阵，主要用来求矩阵的逆。在使用直接法，要注意从空间、时间两方面对算法进行优化。高斯-塞德尔迭代： 列主元高斯消去法：列主元消元 回代图4.2列主元的约当消去图4.1G-S迭代算法流程图约当消去 4 算法描述Gauss -

22、 Seidel算法见流程图选列主元高斯消去法见流程图选列主元约当消去法见流程图5计算用例的参考输出实验1输出参考书本p146页 表5-2实验2参考输出：三个企业产值X=114458, 65395, 85111实验3参考输出：椭圆曲线系数a=0.0508 -0.0702 0.0381 -0.4531 0.2643实验4参考输出：18实验5参考输出：图4.3列主元的高斯消去流程图不可逆6、代码#include<iostream>#include<fstream>#define N 3using namespace std;int main()ifstream file(&q

23、uot;d:/a.txt");double xNN+1;double resultN;for(int i=0;i<N;i+)for(int j=0;j<N+1;j+)file>>xij;for(i=0;i<N;i+)for(int j=i+1;j<N+1;j+)xij=xij/xii;for(j=i+1;j<N;j+)for(int m=i+1;m<N+1;m+)xjm=xjm-xji*xim;for(i=N-1;i>=0;i-)resulti=xiN;for(int j=N-1;j>i;j-)resulti=resulti

24、-resultj*xij;for(i=0;i<N;i+)cout<<resulti<<' 'return 0;7、输出实验结果如下：8、分析体会消元过程中每次都选择绝对值最大者作为主元，这是高斯消去法很重要的一点。用高斯消去法计算高次较为快捷方便，这种算法可以用来解决所有线性方程组。即使一个方程组不能被化为一个三角形的格式，高斯消元法仍可找出它的解。学好高斯消去法对于解线性方程组是十分重要的。实验三 高斯-赛德尔迭代1 实验目的(5) 熟悉求解线性方程组的有关理论和方法；(6) 能编程实现雅可比及高斯-塞德尔迭代法、列主元高斯消去法、约当消去，追赶

25、法(7) 通过测试，进一步了解各种方法的优缺点(8) 根据不同类型的方程组，选择合适的数值方法2 实验内容（1）用Gauss - Seidel 迭代法求解方程组输入：系数矩阵A，最大迭代次数N，初始向量，误差限e输出：解向量（2）一个城镇有三个主要生产企业：煤矿、电厂和地方铁路作为它的经济系统.煤矿：生产价值1元的煤，需消耗0.25元的电费和0.35元的运输费；电厂：生产价值1元的电，需消耗0.40元的煤费、0.05元的电费和0.10元的运输费；铁路：而提供价值1元的铁路运输服务，则需消耗0.45元的煤、0.10元的电费和0.10元的运输费.在某个星期内，除了这三个企业间的彼此需求，煤矿得到5

26、0000元的订单，电厂得到25000元的电量供应要求，而地方铁路得到价值30000元的运输需求.提问：这三个企业在这星期各应生产多少产值才能满足内外需求？提示：“投入”和“产出”之间的平衡关系：物资消耗和新创造的价值等于它生产的总产值. 输入：物资消耗和新创造的价值输出：三个企业生产总值（3）一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立一个以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（1天文单位为地球到太阳的平均距离：9300万里）。他五个不同时间对小行星作五次观测，得到轨道上五个点的坐标分别为（5.764,0.648）（6.286,1.202）（6.759,1.82

27、3）（7.168,2.562）与（7.408，3.360）。由开普勒第一定律知小行星轨道为一椭圆，试建立它的方程。（4）用选主元高斯消去求行列式值提示：D.E. 消元结果直接存储在系数矩阵中F. 当消元过程发生两行对调的情况为偶数次时，行列式值为对角线乘积，否则为对角线乘积的相反数（5）用选主元约当消去分别对矩阵求其逆矩阵，若不可逆输出奇异标志和提示：3 算法基本原理无论是三次样条还是拟合问题最终都归结为线性方程组，求解线性方程组在数值分析中非常重要，在工程计算中也不容忽视。线性方程组大致分迭代法和直接法。只有收敛条件满足时，才可以进行迭代。雅可比及高斯-塞德尔是最基本的两类迭代方法，最大区别

28、是迭代过程中是否引用新值进行剩下的计算。消元是最简单的直接法，并且也十分有效的，列主元高斯消去法对求解一般的线性方程组都适用，同时可以用来求矩阵对应的行列式。约当消去实质是经过初等行变换将系数矩阵化为单位阵，主要用来求矩阵的逆。在使用直接法，要注意从空间、时间两方面对算法进行优化。高斯-塞德尔迭代： 列主元高斯消去法：列主元消元 回代图4.2列主元的约当消去图4.1G-S迭代算法流程图约当消去 4 算法描述Gauss - Seidel算法见流程图选列主元高斯消去法见流程图选列主元约当消去法见流程图5计算用例的参考输出实验1输出参考书本p146页 表5-2实验2参考输出：三个企业产值X=1144

29、58, 65395, 85111实验3参考输出：椭圆曲线系数a=0.0508 -0.0702 0.0381 -0.4531 0.2643实验4参考输出：18实验5参考输出：图4.3列主元的高斯消去流程图不可逆6、代码：#include<iostream>#include<fstream>using namespace std;int main()double xishu34;double x3=0;ifstream in("D:/1.txt");for(int i=0;i<3;i+)for(int j=0;j<4;j+)in>>

30、;xishuij;in.close();/ofstream out("d:/2.txt");int t=0;while(t+<7)for(int i=0;i<3;i+)xi=xishui3;for(int j=0;j<3;j+)if(!(j=i)xi-=xishuij*xj;xi/=xishuii;out<<xi<<' 'out<<endl;return 0;7、输入与输出8、分析体会本次试验通过理论与实例对线性方程组的解法、收敛性及误差分析进行了探讨.在对线性方程组数值解法的讨论下用到了高斯-赛德尔迭代

31、法，让我进一步研究和总结了高斯-赛德尔迭代法的理论与应用，使我在分析问题与编辑程序时能更好的把握对高斯-赛德尔迭代法的应用。实验四 最小二乘法和三次样条插值1 实验目的(1) 明确插值多项式和分段插值多项式各自的优缺点；(2) 编程实现拉格朗日插值算法，分析实验结果体会高次插值产生的龙格现象；(3) 理解最小二乘拟合，并编程实现线性拟合，掌握非线性拟合转化为线性拟合的方法(4) 运用常用的插值和拟合方法解决实际问题。2 实验内容(1) 对于要求选取11个等距插值节点，分别采用拉格朗日插值和分段线性插值，计算x为0.5, 4.5处的函数值并将结果与精确值进行比较。输入：区间长度，n(即n+1个节

32、点)，预测点输出：预测点的近似函数值，精确值，及误差（2）对于给定的一元函数 的n+1个节点值 。数据如下：0.40.550.650.800.951.050.410750.578150.696750.901.001.25382求五次拉格朗日多项式和分段三次插值多项式，计算 的值。输入：数据点集及个数，预测点。输出：预测点的近似函数值（3）在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间的拟合曲线。分）0 5 10 15 20 25 30 35 40 45 50 550 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64输

33、入： 数据点集及个数。输出：拟合函数，并打印出与的误差提示：近似解析表达式为（4）用线性拟合编程实现给药方案背景知识：一种新药用于临床之前，必须设计给药方案.临床上，每种药物有一个最小有效浓度c1和一个最大有效浓度c2。设计给药方案时，要使血药浓度 保持在c1c2之间。 本题设c1=10，c2=25(ug/ml).在实验方面，对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血药,测得血药浓度c(ug/ml)如下表:t (h) 0.25 0.5 1 1.5 2 3 4 6 8c (mg/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.

34、45 5.24 3.01问题：B. 在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律C. 设计给药方案：每次注射剂量多大；间隔时间多长提示：A. 血药浓度的变化规律符合指数关系B. 将指数关系转化为线性关系3 算法基本原理当精确函数 y = f(x) 非常复杂或未知，在一系列节点 x0 xn 处测得函数值 y0 = f(x0), yn = f(xn)， 希望由此构造一个简单易算的近似函数 g(x) » f(x)，满足条件g(xi) = f(xi) (i = 0, n)。这里的 g(x) 称为f(x) 的插值函数，由插值函数可以去近似估计f(x)在一些未知

35、点处的函数值。最常用的插值函数是多项式插值。 所谓多项式插值即求 n 次多项式使得基于基函数的拉格朗日插值是构造多项式插值最基本方法。也是推导数值微积分和微分方程数值解的公式的理论基础。拉格朗日插值公式为： 其中当结点比较多，次数较高的插值多项式往往发生Runge现象，分段低次插值是避免Runge现象的重要手段。分段一次插值将整个区间分段，在每个小区间上，用一次多项式逼近 f (x)，直观上即用折线代替曲线，只要区间足够小就可以保证很好的逼近效果，但曲线缺乏光滑性。除了插值，逼近复杂函数f(x)的另一种方法是拟合。结点数据比较多，或者不够精确，往往采用拟合的方法。用来拟合的简单函数p(x)，不

36、需严格通过给定结点，但在这些结点处总体误差达到极小，即满足最小二乘拟合的条件。线性拟合是最简单的拟合方法，许多非线性问题都可以转化为线性拟合。令线性函数为a0+a1x，要确定a0和a1，最终求解如下线性方程组：图1.1拉格朗日插值算法4 算法描述 (2) 拉格朗日插值算法见流程图(3) 分段线性插值算法要注意以下关键点：A. 定位：预测点x在第几个区间？B. 假定x在第k个区间，即 一次插值公式为(3) 线性拟合算法在没有讨论线性方程组数值解法的前提下，用最原始的消元法求解线性方程组5 计算用例的参考输出（1）实验1参考输出如下X y(精确) y(拉格朗日) y(分段线性) 误差(拉) 误差(

37、分)0.5000000.8000000.8434080.750000 -0.043408 0.0500004.5000000.0470591.578721 0.048643 -1.537662 -0.001584（2）实验2参考输出如下模型中参数 k=0.2347 v=15.02故可制定给药方案：首次注射 375 mg， 其余每次注射 225 mg， 注射的间隔时间为 4 小时6 代码最小二乘法：#include <iostream>#include <fstream>using namespace std;void main()ifstream in("d:

38、/2.txt"); float *a,x,result; int n; float k,b; void smallest(float *a,int n,float*k,float *b);/ cout<<"输入点数：" in>>n; a=new float *n; for(int i=0;i<n;i+)ai=new float2;/ cout<<"输入各个点："<<endl; for(i=0;i<n;i+)for(int j=0;j<2;j+)in>>aij; sma

39、llest(a,n,&b,&k);/ cout<<"输入自变量的值：" in>>x; result=b+k*x; out<<"拟合曲线方程:nY="<<b<<"+"<<k<<"X"<<endl; out<<"x="<<x<<"时运算结果是："<<result<<endl;void smallest(floa

40、t *a,int n,float*k,float *b) float can23=n,0,0,0,0,0; for(int i=0;i<n;i+) can01+=ai0; can02+=ai1; can11+=ai0*ai0; can12+=ai0*ai1; can10=can01; for(i=2;i>=0;i-)can1i=can1i-can0i*can10/can00; *b=can12/can11; *k=(can02-can01*b)/can00;三次样条插值：#include<iostream>#include<iomanip>#include<fstream>using namespace std;const int MAX = 50;float xMAX, yMAX, hMAX;float cMAX, aMAX, fxymMAX;float f(int x1, int x2, int x3) float a = (yx3 - yx2) / (xx3 - xx2); float b = (yx2 - yx1