《计算方法》 试验报告

　 班级：

　学号：

　姓名：

　实验一、牛顿下山法 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 = ； x = ；x = ；x = x -

　输入：初值，误差限，迭代最大次数 输出：近似根、实际迭代次数 3

　算法基本原理 求非线性方程组的解是科学计算常遇到的问题，有很多实际背景．各种算法层出不穷，其中迭代是主流算法。只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。因此设计算法之前，对于一般迭代进行收敛性的判断是至关重要的。牛顿法也叫切线法，是迭代算法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代 至关重要。当初值选取不当可以采用牛顿下山算法进行纠正。

　一般迭代：

　 牛顿公式：

　牛顿下山公式：

　图3.1一般迭代算法流程图

　下山因子

　下山条件 4

　算法描述 一般迭代算法见流程图 图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 double[n+1];

　 //动态分配内存

　 xr = new double[n];

　 xi = new double[n];

　}

　 void input ();

　//由文件读入代数方程系数

　void srrt_root ();

　//执行牛顿下山法

　void output ();

　 //输出根到文件并显示

　~srrt ()

　{

　delete [] a, xr, xi;

　 }

　 };

　void srrt::input ()

　//由文件读入代数方程系数

　 {

　int

　i;

　char str1[20];

　cout <<“＼n输入文件名:

　“;

　cin >>str1;

　ifstream

　fin (str1);

　if (!fin)

　{ cout <<“＼n不能打开这个文件 “ <<str1 <<endl; exit(1); }

　for (i=n; i>=0; i--)

　fin >>a[i];

　//读入代数方程系数

　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(a[m])+1.0==1.0)) m=m-1;

　 if (m<=0)

　 {

　cout <<“＼n程序工作失败！“ <<endl;

　return;

　}

　 for (i=0; i<=m; i++) a[i]=a[i]/a[m];

　 for (i=0; i<=m/2; i++)

　 {

　w=a[i]; a[i]=a[m-i]; a[m-i]=w;

　}

　 k=m; is=0; w=1.0;

　 jt=1;

　 while (jt==1)

　 {

　pq=fabs(a[k]);

　while (pq<1.0e-12)

　 {

　 xr[k-1]=0.0; xi[k-1]=0.0; k=k-1;

　 if (k==1)

　 {

　xr[0]=-a[1]*w/a[0]; xi[0]=0.0;

　 return;

　 }

　 pq=fabs(a[k]);

　 }

　q=log(pq); q=q/(1.0*k); q=exp(q);

　 p=q; w=w*p;

　 for (i=1; i<=k; i++)

　 {

　a[i]=a[i]/q; q=q*p; }

　 x=0.0001; x1=x; y=0.2; y1=y; dx=1.0;

　 g=1.0e+37;

　l40:

　 u=a[0]; v=0.0;

　 for (i=1; i<=k; i++)

　 {

　 p=u*x1; q=v*y1;

　 pq=(u+v)*(x1+y1);

　 u=p-q+a[i]; 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=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 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.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;

　 xr[k-1]=x*w;

　 xi[k-1]=-y*w;

　 k=k-1;

　}

　 for (i=1; i<=k; i++)

　{

　a[i]=a[i]-a[i-1]*p;

　 a[i+1]=a[i+1]-a[i-1]*q;

　}

　 xr[k-1]=x*w; xi[k-1]=y*w;

　 k=k-1;

　 if (k==1)

　{ xr[0]=-a[1]*w/a[0]; xi[0]=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*x; q=x*x+y*y;

　 xr[k-1]=x*w;

　 xi[k-1]=-y*w;

　 k=k-1;

　 }

　 for (i=1; i<=k; i++)

　 {

　a[i]=a[i]-a[i-1]*p;

　 a[i+1]=a[i+1]-a[i-1]*q;

　 }

　 xr[k-1]=x*w; xi[k-1]=y*w;

　 k=k-1;

　 if (k==1)

　 { xr[0]=-a[1]*w/a[0]; xi[0]=0.0; }

　}

　 else

　 {

　u1=k*a[0]; 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)*a[i-1];

　 v1=pq-p-q;

　 }

　 p=u1*u1+v1*v1;

　 if (p<=1.0e-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; q=0.0; }

　 else

　{

　 p=-2.0*x; q=x*x+y*y;

　xr[k-1]=x*w;

　xi[k-1]=-y*w;

　k=k-1;

　}

　 for (i=1; i<=k; i++)

　{

　a[i]=a[i]-a[i-1]*p;

　 a[i+1]=a[i+1]-a[i-1]*q;

　}

　 xr[k-1]=x*w; xi[k-1]=y*w;

　 k=k-1;

　 if (k==1)

　{ xr[0]=-a[1]*w/a[0]; xi[0]=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>=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.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;

　xr[k-1]=x*w;

　xi[k-1]=-y*w;

　k=k-1;

　}

　 for (i=1; i<=k; i++)

　{

　a[i]=a[i]-a[i-1]*p;

　 a[i+1]=a[i+1]-a[i-1]*q;

　}

　 xr[k-1]=x*w; xi[k-1]=y*w;

　 k=k-1;

　 if (k==1)

　{ xr[0]=-a[1]*w/a[0]; xi[0]=0.0; }

　 }

　 }

　 }

　 if (k==1) jt=0;

　 else jt=1;

　 }

　 }

　void srrt::output ()

　 //输出根到文件并显示

　 {

　int

　k;

　char str2[20];

　cout <<“＼n输出文件名:

　“;

　cin >>str2;

　ofstream fout (str2);

　if (!fout)

　{ cout <<“＼n不能打开这个文件 “ <<str2 <<endl; exit(1); }

　for (k=0; k<n; k++)

　{

　 fout <<xr[k] <<“

　“ <<xi[k] <<endl;

　 cout <<xr[k] <<“

　+j

　“ <<xi[k] <<endl;

　}

　fout.close ();

　 }

　void main ()

　//主函数

　 {

　srrt

　root(6);

　root.input ();

　//由文件读入代数方程系数

　root.srrt_root ();

　//执行牛顿下山法

　root.output ();

　 //输出根到文件并显示

　 } 6、输入输出 输出结果如下：

　7、分析体会

　牛顿下山法作为计算方法课程中重要的知识点，在看书学习时较易大致理解其思路，但上级编写代码时却是有难度的。在学习借鉴后用如上代码能够实现基本功能，在实验过程中也加深了我对牛顿下山法的理解。

　实验二 高斯消去法 1 实验目的 (1) 熟悉求解线性方程组的有关理论和方法； (2) 能编程实现雅可比及高斯-塞德尔迭代法、列主元高斯消去法、约当消去，追赶法 (3) 通过测试，进一步了解各种方法的优缺点 (4) 根据不同类型的方程组，选择合适的数值方法 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元的运输费. 在某个星期内，除了这三个企业间的彼此需求，煤矿得到50000元的订单，电厂得到25000元的电量供应要求，而地方铁路得到价值30000元的运输需求. 提问：这三个企业在这星期各应生产多少产值才能满足内外需求？ 提示：

　“投入”和“产出”之间的平衡关系：物资消耗和新创造的价值等于它生产的总产值.

　 输入：物资消耗和新创造的价值 输出：三个企业生产总值

　（3）一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立一个以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（1天文单位为地球到太阳的平均距离：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

　算法基本原理 无论是三次样条还是拟合问题最终都归结为线性方程组，求解线性方程组在数值分析中非常重要，在工程计算中也不容忽视。

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

　 高斯-塞德尔迭代：

　 列主元高斯消去法：列主元 消元

　 回代 图4.2列主元的约当消去

　图4.1G-S迭代算法流程图

　约当消去

　 4

　算法描述 Gauss - 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 3 using namespace std; int main() {

　ifstream file(“d://a.txt“);

　double x[N][N+1];

　double result[N];

　for(int i=0;i<N;i++)

　{

　 for(int j=0;j<N+1;j++)

　 {

　file>>x[i][j];

　 }

　}

　for(i=0;i<N;i++)

　{

　 for(int j=i+1;j<N+1;j++)

　 {

　x[i][j]=x[i][j]/x[i][i];

　 }

　 for(j=i+1;j<N;j++)

　 {

　for(int m=i+1;m<N+1;m++)

　{

　 x[j][m]=x[j][m]-x[j][i]*x[i][m];

　}

　 }

　}

　for(i=N-1;i>=0;i--)

　{

　 result[i]=x[i][N];

　 for(int j=N-1;j>i;j--)

　 {

　result[i]=result[i]-result[j]*x[i][j];

　 }

　}

　for(i=0;i<N;i++)

　 cout<<result[i]<<＇ ＇;

　return 0; } 7、输出 实验结果如下：

　 8、分析体会

　消元过程中每次都选择绝对值最大者作为主元，这是高斯消去法很重要的一点。用高斯消去法计算高次较为快捷方便，这种算法可以用来解决所有线性方程组。即使一个方程组不能被化为一个三角形的格式，高斯消元法仍可找出它的解。学好高斯消去法对于解线性方程组是十分重要的。

　实验三 高斯-赛德尔迭代 1 实验目的 (5) 熟悉求解线性方程组的有关理论和方法； (6) 能编程实现雅可比及高斯-塞德尔迭代法、列主元高斯消去法、约当消去，追赶法 (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元的运输费. 在某个星期内，除了这三个企业间的彼此需求，煤矿得到50000元的订单，电厂得到25000元的电量供应要求，而地方铁路得到价值30000元的运输需求. 提问：这三个企业在这星期各应生产多少产值才能满足内外需求？ 提示：

　“投入”和“产出”之间的平衡关系：物资消耗和新创造的价值等于它生产的总产值.

　 输入：物资消耗和新创造的价值 输出：三个企业生产总值

　（3）一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立一个以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（1天文单位为地球到太阳的平均距离：9300万里）。他五个不同时间对小行星作五次观测，得到轨道上五个点的坐标分别为（5.764,0.648）（6.286,1.202）（6.759,1.823）（7.168,2.562）与（7.408，3.360）。由开普勒第一定律知小行星轨道为一椭圆，试建立它的方程。

　（4）用选主元高斯消去求行列式值 提示：

　D. E. 消元结果直接存储在系数矩阵中

　F. 当消元过程发生两行对调的情况为偶数次时，行列式值为对角线乘积，否则为对角线乘积的相反数

　（5）用选主元约当消去分别对矩阵求其逆矩阵，若不可逆输出奇异标志

　 和 提示：

　3

　算法基本原理 无论是三次样条还是拟合问题最终都归结为线性方程组，求解线性方程组在数值分析中非常重要，在工程计算中也不容忽视。

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

　 高斯-塞德尔迭代：

　 列主元高斯消去法：列主元 消元

　 回代 图4.2列主元的约当消去

　图4.1G-S迭代算法流程图

　约当消去

　 4

　算法描述 Gauss - 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> using namespace std; int main(){

　 double xishu[3][4];

　 double x[3]={0};

　 ifstream in(“D://1.txt“);

　 for(int i=0;i<3;i++)

　 {

　 for(int j=0;j<4;j++)

　 {

　 in>>xishu[i][j];

　 }

　 }

　 in.close();

　 //ofstream out(“d://2.txt“);

　 int t=0;

　 while(t++<7){

　 for(int i=0;i<3;i++){

　 x[i]=xishu[i][3];

　 for(int j=0;j<3;j++){

　 if(!(j==i))

　 {

　 x[i]-=xishu[i][j]*x[j];

　 }

　 }

　 x[i]/=xishu[i][i];

　 out<<x[i]<<＇ ＇;

　 }

　 out<<endl;

　 }

　 return 0; } 7、输入与输出

　8、分析体会

　本次试验通过理论与实例对线性方程组的解法、收敛性及误差分析进行了探讨.在对线性方程组数值解法的讨论下用到了高斯-赛德尔迭代法，让我进一步研究和总结了高斯-赛德尔迭代法的理论与应用，使我在分析问题与编辑程序时能更好的把握对高斯-赛德尔迭代法的应用。

　实验四 最小二乘法和三次样条插值 1

　实验目的 (1) 明确插值多项式和分段插值多项式各自的优缺点； (2) 编程实现拉格朗日插值算法，分析实验结果体会高次插值产生的龙格现象； (3) 理解最小二乘拟合，并编程实现线性拟合，掌握非线性拟合转化为线性拟合的方法 (4) 运用常用的插值和拟合方法解决实际问题。

　2

　实验内容 (1) 对于 要求选取11个等距插值节点，分别采用拉格朗日插值和分段线性插值，计算x为0.5, 4.5处的函数值并将结果与精确值进行比较。

　输入：区间长度，n(即n+1个节点)，预测点 输出：预测点的近似函数值，精确值，及误差 （2）对于给定的一元函数

　的n+1个节点值

　。

　数据如下：

　 0.4 0.55 0.65 0.80 0.95 1.05

　0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次拉格朗日多项式和分段三次插值多项式，计算 的值。

　输入：数据点集及个数，预测点。

　输出：预测点的近似函数值

　（3）在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间的拟合曲线。

　分） 0

　 5

　10

　 15

　20

　 25

　 30

　35

　 40

　45

　50

　55

　0

　1.27

　2.16

　2.86

　3.44

　3.87

　4.15

　4.37

　4.51

　4.58

　 4.02

　4.64 输入：

　数据点集及个数。

　输出：拟合函数，并打印出与的误差 提示：近似解析表达式为

　（4）用线性拟合编程实现给药方案 背景知识：一种新药用于临床之前，必须设计给药方案. 临床上，每种药物有一个最小有效浓度c1和一个最大有效浓度c2。

　设计给药方案时，要使血药浓度 保持在c1~c2之间。

　本题设c1=10，c2=25(ug/ml). 在实验方面，对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血药,测得血药浓度c(ug/ml)如下表: t (h)

　 0.25

　0.5

　 1

　 1.5

　 2

　3

　 4

　 6

　8 c (mg/ml)

　 19.21

　18.15

　15.36

　14.10

　12.89

　 9.32

　 7.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)在一些未知点处的函数值。最常用的插值函数是多项式插值。

　 所谓多项式插值即求 n 次多项式使得 基于基函数的拉格朗日插值是构造多项式插值最基本方法。也是推导数值微积分和微分方程数值解的公式的理论基础。

　拉格朗日插值公式为：

　 其中 当结点比较多，次数较高的插值多项式往往发生Runge现象，分段低次插值是避免Runge现象的重要手段。分段一次插值将整个区间分段，在每个小区间上，用一次多项式逼近 f (x)，直观上即用折线代替曲线，只要区间足够小就可以保证很好的逼近效果，但曲线缺乏光滑性。

　除了插值，逼近复杂函数f(x)的另一种方法是拟合。结点数据比较多，或者不够精确，往往采用拟合的方法。用来拟合的简单函数p(x)，不需严格通过给定结点，但在这些结点处总体误差达到极小，即满足最小二乘拟合的条件。线性拟合是最简单的拟合方法，许多非线性问题都可以转化为线性拟合。令线性函数为a0+a1x，要确定a0和a1，最终求解如下线性方程组：

　 图1.1拉格朗日插值算法 4

　算法描述

　 (2) 拉格朗日插值算法见流程图 (3) 分段线性插值算法要注意以下关键点：

　A. 定位：预测点x在第几个区间？ B. 假定x在第k个区间， 即

　 一次插值公式为 (3) 线性拟合算法 在没有讨论线性方程组数值解法的前提下，用最原始的消元法求解线性方程组 5

　计算用例的参考输出 （1）实验1参考输出如下 X

　y(精确)

　 y(拉格朗日) y(分段线性)

　误差(拉)

　误差(分) 0.500000

　0.800000

　0.843408

　0.750000

　 -0.043408

　 0.050000 4.500000

　0.047059

　1.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://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++)

　 a[i]=new float[2];

　// cout<<“输入各个点：“<<endl;

　 for(i=0;i<n;i++)

　for(int j=0;j<2;j++)

　in>>a[i][j];

　 smallest(a,n,&b,&k);

　// cout<<“输入自变量的值：“;

　 in>>x;

　 result=b+k*x;

　 out<<“拟合曲线方程:＼nY=“<<b<<“+“<<k<<“X“<<endl;

　 out<<“x=“<<x<<“时运算结果是：“<<result<<endl; } void smallest(float **a,int n,float*k,float *b) {

　 float can[2][3]={n,0,0,0,0,0};

　 for(int i=0;i<n;i++)

　 {

　 can[0][1]+=a[i][0];

　 can[0][2]+=a[i][1];

　 can[1][1]+=a[i][0]*a[i][0];

　 can[1][2]+=a[i][0]*a[i][1];

　 }

　 can[1][0]=can[0][1];

　 for(i=2;i>=0;i--)

　 can[1][i]=can[1][i]-can[0][i]*can[1][0]/can[0][0];

　 *b=can[1][2]/can[1][1];

　 *k=(can[0][2]-can[0][1]**b)/can[0][0]; } 三次样条插值：

　#include<iostream> #include<iomanip> #include<fstream> using namespace std;

　const int MAX = 50; float x[MAX], y[MAX], h[MAX]; float c[MAX], a[MAX], fxym[MAX];

　float f(int x1, int x2, int x3){

　 float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);

　 float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);

　 return (a - b)/(x[x3] - x[x1]); }

　//求差分 void cal_m(int n){

　//用追赶法求解出弯矩向量M……

　 float B[MAX];

　 B[0] = c[0] / 2;

　 for(int i = 1; i < n; i++)

　 B[i] = c[i] / (2 - a[i]*B[i-1]);

　fxym[0] = fxym[0] / 2;

　 for(i = 1; i <= n; i++)

　 fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);

　 for(i = n-1; i >= 0; i--)

　 fxym[i] = fxym[i] - B[i]*fxym[i+1]; } void printout(int n);

　int main(){

　int times=0;

　 int n,i; char ch;

　ifstream in(“d:＼＼1.txt“);

　 do{

　 //

　cout<<“Please put in the number of the dots:“;

　 in>>n;

　 for(i = 0; i <= n; i++){

　//

　 cout<<“Please put in X“<<i<<＇:＇;

　 in>>x[i];

　//cout<<endl;

　//

　 cout<<“Please put in Y“<<i<<＇:＇;

　 in>>y[i]; //cout<<endl;

　 }

　 for(i = 0; i < n; i++)

　//求 步长

　 h[i] = x[i+1] - x[i];

　// cout<<“Please 输入边界条件＼n 1: 已知两端的一阶导数＼n 2:两端的二阶导数已知＼n 默认:自然边界条件＼n“;

　 int t;

　 float f0, f1;

　 in>>t;

　 switch(t){

　 case 1://cout<<“Please put in Y0＼＇ Y“<<n<<“＼＇＼n“;

　 in>>f0>>f1;

　 c[0] = 1; a[n] = 1;

　 fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];

　 fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];

　 break;

　 case 2://cout<<“Please put in Y0＼“ Y“<<n<<“＼“＼n“;

　 in>>f0>>f1;

　 c[0] = a[n] = 0;

　 fxym[0] = 2*f0; fxym[n] = 2*f1;

　 break;

　 default:cout<<“不可用＼n“;//待定

　 };//switch

　 for(i = 1; i < n; i++)

　 fxym[i] = 6 * f(i-1, i, i+1);

　 for(i = 1; i < n; i++){

　 a[i] = h[i-1] / (h[i] + h[i-1]);

　 c[i] = 1 - a[i];

　 }

　 a[n] = h[n-1] / (h[n-1] + h[n]);

　 cal_m(n);

　 if(times==0)

　 cout<<“已知两端的一阶导数时：＼n“;

　 else

　cout<<“已知两端的二阶导数时：＼n“;

　 printout(n);

　 times++;

　 //

　cout<<“Do you to have anther try ? y/n :“;

　 in>>ch;

　 }while(ch == ＇y＇ || ch == ＇Y＇);

　 return 0; } void printout(int n){

　 cout<<setprecision(6);

　 for(int i = 0; i < n; i++){

　 cout<<i+1<<“: [“<<x[i]<<“ , “<<x[i+1]<<“]＼n“<<“＼t“;

　 /*

　 cout<<fxym[i]/(6*h[i])<<“ * (“<<x[i+1]<<“ - x)^3 + “<<<<“ * (x - “<<x[i]<<“)^3 + “

　 <<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<“ * (“<<x[i+1]<<“ - x) + “

　 <<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<“(x - “<<x[i]<<“)＼n“;

　 cout<<endl;*/

　 float t = fxym[i]/(6*h[i]);

　 if(t > 0)cout<<t<<“*(“<<x[i+1]<<“ - x)^3“;

　 else cout<<-t<<“*(x - “<<x[i+1]<<“)^3“;

　 t = fxym[i+1]/(6*h[i]);

　 if(t > 0)cout<<“ + “<<t<<“*(x - “<<x[i]<<“)^3“;

　 else cout<<“ - “<<-t<<“*(x - “<<x[i]<<“)^3“;

　 cout<<“＼n＼t“;

　 t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];

　 if(t > 0)cout<<“+ “<<t<<“*(“<<x[i+1]<<“ - x)“;

　 else cout<<“- “<<-t<<“*(“<<x[i+1]<<“ - x)“;

　 t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];

　 if(t > 0)cout<<“ + “<<t<<“*(x - “<<x[i]<<“)“;

　 else cout<<“ - “<<-t<<“*(x - “<<x[i]<<“)“;

　 cout<<endl<<endl;

　 }

　 cout<<endl; } 7、实验运行结果：

　最小二乘法：

　 三次样条差值：

　 8、分析体会

　在实际应用中，有时需解决的问题中运算很复杂且运算量也很大，这时我们就需要借助计算机的帮助解决实际问题，即利用最小二乘法和样条插值在C++语言中的编程实现，可以帮助我更好的理解计算方法，上级实践有益于我的进一步学习。因为在我们现实生活中我们需要通过已有的数据来发掘事物本身的内在规律，或者模拟出相应的数学模型来解决。所以这就需要我们用到上述相关知识来完成。

　 实验五、龙贝格算法 1 实验目的 (1) 熟悉复化梯形方法、复化Simpson方法、梯形递推算法、龙贝格算法； (2) 能编程实现复化梯形方法、复化Simpson方法、梯形递推算法、龙贝格算法； (3) 理解并掌握自适应算法和收敛加速算法的基本思想； (4) 分析实验结果体会各种方法的精确度，建立计算机求解定积分问题的感性认识 2

　实验内容 (1) 用龙贝格算法计算 输入：积分区间，误差限 输出：序列Tn，Sn,Cn,Rn及积分结果（参考书本P81的表2-4）

　（2）设计方案计算国土面积

　为了计算瑞士国土的面积,首先对地图作了如下测量:以西向东方向为x轴,由南向北方向为y轴,选择方便的原点,并将从最西边界到最东边界在x轴上的区间适当地划分为若干 段,在每个分点的y方向测出南边界点和北边界点的y坐标，数据如表(单位mm): x

　7.0

　10.5

　13.0

　17.5

　34.0

　40.5

　44.5

　48.0

　56.0

　y1

　44

　45

　47

　50

　50

　38

　30

　30

　34 y2

　44

　59

　70

　72

　93

　 100

　 110

　110

　 110 x

　 61.0

　68.5

　76.5

　80.5

　91.0

　96.0 101.0 104.0 106.5 y1

　36

　34

　41

　45

　46

　43

　 37

　 33

　 28 y2

　 117

　 118

　 116

　 118

　 118

　 121

　124

　 121

　 121 x

　111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 y1

　32

　65

　55

　 54

　 52

　50

　66

　 66

　 68 y2

　 121

　 122

　 116

　 83

　 81

　82

　86

　 85

　 68

　瑞士地图的外形如图(比例尺18mm:40km)

　问题：

　试由测量数据计算瑞士国土的近似面积，并与其精确值41288平方公里比较 提示：

　A. 将计算国土面积问题转化为求积分问题 B. 当被积函数为表格形式（仅提供一组数据点，函数形式未知），使用梯形公式最为简单 3

　算法基本原理 在许多实际问题中，常常需要计算定积分的值。根据微积分学基本定理，若被积函数f(x)在区间[a,b]上连续，只要能找到f(x)的一个原函数F(x)，便可利用牛顿-莱布尼兹公式求得积分值。

　但是在实际使用中，往往遇到如下困难，而不能使用牛顿-莱布尼兹公式。

　（1） 找不到用初等函数表示的原函数 （2） 虽然找到了原函数，但因表达式过于复杂而不便计算 （3） f(x)是由测量或计算得到的表格函数 由于以上种种困难，有必要研究积分的数值计算问题。

　利用插值多项式 则积分转化为，显然易算。称为插值型求积公式。最简单的插值型求积公式是梯形公式和Simpson公式，。当求积结点提供较多，可以分段使用少结点的梯形公式和Simpson公式，并称为复化梯形公式、复化Simpson公式。如步长未知，可以通过误差限的控制用区间逐次分半的策略自动选取步长的方法称自适应算法。梯形递推公式给出了区间分半前后的递推关系。由梯形递推公式求得梯形序列，相邻序列值作线性组合得Simpson序列, Simpson序列作线性组合得柯特斯序列, 柯特斯序列作线性组合的龙贝格序列。若|R2-R1|<e,则输出R2;否则…依此类推。如此加工数据的过程叫龙贝格算法，如下图所示：

　 图2.1龙贝格算法数据加工流程

　复化梯形公式 复化Simpson公式 梯形递推公式 加权平均公式：

　 龙贝格算法大大加快了误差收敛的速度，由梯形序列O(h2) 提高到龙贝格序列的O(h8) 4

　算法描述 (1) 梯形递推算法见流程图 图2.2梯形递推算法流程图

　图2.3龙贝格算法流程图

　(2) 龙贝格算法见流程图 5 计算用例的参考输出 实验1参考输出如下 图2.4 龙贝格参考输出

　6、代码 h=b-a;

　 T2=T1=(F(a)+F(b))*h/2;

　 //进行梯形公式数据装载

　for(int i=0;i<N+3;i++)//N为龙贝格数据个数

　{

　 T[i]=T2;

　 sum=0;

　 x=a+h/2;

　 while(x<b)

　 {

　 sum=F(x)+sum;

　 x+=h;

　 }

　 T2=(T1+h*sum)/2;

　 h=h/2;

　 T1=T2;

　 }

　 //进行梯形公式到辛甫生的转化

　for(int i=0;i<N+2;i++)

　 S[i]=(4*T[i+1]-T[i])/3;

　 //进行辛甫生到柯特斯的转化

　for(int i=0;i<N+1;i++)

　 C[i]=(16*S[i+1]-S[i])/15;

　 //进行柯特斯到龙贝格的转化

　for(int i=0;i<N;i++)

　 R[i]=(64*C[i+1]-C[i])/63;

　 7、实验运行如图：

　 8、分析体会

　用龙贝格积分法进行近似求积，其原理与埃特金插值有些类似，进行线性整合后使结果具有高精度的求积效果。在实际过程中，由于对于评判合理步长的困难，我们常采取变步长的办法进行计算，使结果满足精度要求。龙贝格运算过程：梯形公式装载数据→辛甫生化→柯特斯化→龙贝格高精度化。通过实验我牢记了龙贝格算法的演算过程，更体会到了这一方法的高效和高精度。

　实验六 欧拉算法、龙格库塔、亚当姆斯 1

　实验目的 (1) 熟悉数值微分中Euler法，改进Euler法，Rung-Kutta方法； (2) 能编程实现Euler法，改进Euler法，Rung-Kutta方法； (3) 通过实验结果分析各个算法的优缺点; (4) 明确步长对算法的影响并理解变步长的Rung-Kutta方法 2

　实验内容 （1）

　0 < x<1 取h=0.1时用Euler法，改进Euler法，Rung-Kutta方法求其数值解并与精确解进行比较。

　输入：求解区间，初值，数值解个数 输出：数值解

　（2）小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. 输入：相关数据 输出：引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． 提示：可得到火箭飞行过程的方程：

　初始条件为

　3

　算法基本原理 在许多科学技术问题中，建立的模型常常以常微分方程的形式表示。然而，除了少数特殊类型的常微分方程能用解析方法求其精确解外，要给出一般方程解析解的表达式是困难的。所以只能用近似方法求其数值解，在实际工作中常用计算机求常微分方程的数值解。所谓常微分方程的数值解即对于常微分方程初值问题计算出在一系列节点 a = x0< x1<…< xn= b 处的未知函数 y(x)近似值y0,y1,…yn，即找到一系列离散点（x0,y0）（x1,y1）（x2,y2）…（xn,yn）近似满足常微分方程。数值解法的基本思想用差商代替导数,实现连续问题离散化，选取不同的差商代替导数可以得到不同公式。

　改进欧拉公式是常用方法之一，包括预测和校正两步。先用欧拉公式进行预报，再将预报值代入梯形公式进行校正，从而达到二阶精度。通过龙格-库塔法我们可以获得更高精度。经典龙格-库塔法即在区间[xn,xn+1]取四点，并对这四点的斜率进行加权平均作为平均斜率，通过泰勒公式寻找使局部截断误差为O(h5)（即4阶精度）的参数满足条件。

　改进的欧拉公式：

　 预测 校正 四阶（经典）龙格-库塔公式

　4

　算法描述 (1) 改进的欧拉算法见流程图 (2)

　经典龙格-库塔算法见流程图 5 计算用例的参考输出 （1）实验1参考输出如下

　图5.3经典龙格库塔算法

　图5.2改进欧拉算法

　 （2）实验2参考输出如下 关闭引擎时（秒），此时火箭的高度h==8322.4米，速度v==254.1728米/秒，加速度为a== 4.0616米/秒，火箭到最高点的时间=51秒，高度=9186.4米。

　6 代码 改进的欧拉格式：

　/*取h=0.2，用改进欧拉方法求解下列初值问题。

　dy/dx=sqrt（2*x*x+3*y*y），y（0）=5,0<=x<=20. 程序：*/ #include <stdio.h> #include <math.h> #define MAXSIZE 50 double f(double x,double y); void main (void) {

　double a,b,h,x[MAXSIZE],y[MAXSIZE];

　long i,n;

　printf(“＼n请输入求解区间a,b:“);

　scanf(“%1f,%1f“,&a,&b);

　printf(“＼n请输入步长h:“);

　scanf(“%1f“,&h);

　n=(long)((b-a)/h);

　x[0]=a;

　printf(“＼n 请输入起点x[0]=%1f处得纵坐标y[0]:“,x[0]);

　scanf(“%1f“,&y[0]);

　for(i=0;i<n;i++)

　{

　x[i+1]=x[i]+h;

　y[i+1]=y[i]+h*f(x[i],y[i]);

　y[i+1]=y[i]+h*(f(x[i],y[i])+f(x[i+1],y[i+1]))/2;

　}

　printf(“＼n 计算结果为：“);

　for(i=0;i<=n;i++)

　printf(“＼nx[%1d]=%1f,y[%1d]=%1f“,i,x[i],i,y[i]); } double f(double x,double y) {

　return(sqrt(2*x*x+3*y*y)); } 龙哥库塔：

　#include <iostream> using namespace std;

　main() {

　int k;

　float X[k],Y[k],F[k];

　float h=0.001;

　double K1,K2,K3,K4;

　float f(float,float);

　X[0]=0.0;

　Y[0]=1.0;

　F[0]=f(X[0],Y[0]);

　for(k=0;k<100;k++)

　{ X[k+1]=X[k]+h;

　K1=h*F[k];

　K2=h*f(X[k]+0.5*h,Y[k]+0.5*K1);

　K3=h*f(X[k]+0.5*h,Y[k]+K2);

　K4=h*f(X[k]+h,Y[k]+K3) ;

　Y[k+1]=Y[k]+(K1+2.0*K2+2.0*K3+K4)/6.0;

　F[k+1]=f(X[k+1],Y[k+1]);

　cout<<X[k]<<endl;

　cout<<Y[k]<<endl;

　cout<<F[k]<<endl;

　}

　 system(“Pause“);

　 return 0; } float f(float a,float b) { float c;

　 c=b-a*2.0/b;

　 return (c);

　} 7、 欧拉算法

　龙格库塔

　亚当姆斯

　8、分析体会

　改进的欧拉算法、龙格库塔算法和亚当姆斯算法是计算方法中的难点，经过了上机实践了，掌握起来会更透彻一点。程序的编制过程中，我体会到了将一个抽象的数值计算方法变为一个具体计算机程序的不易。这个过程，让我对各种数值计算方法的数学公式有了更进一步的认识，因为要编制出计算机程序，就必须清楚数值计算方法的具体计算步骤。

　 学习了这门课，感觉实用性比较大。这门课程也是连接数学与计算机之间的桥梁，计算方法这门学科只有和计算机结合才能体现出其价值。借助计算机我们可以快速解决各种方程的求值，且结果具有高度精确性，这是单一的数学学习做不到的。