工程常用算法 牛顿 分段线性 三次样条 哈尔滨工程大学作业

哈尔滨工程大学机电工程学院综合实践大作业学年学期：2013-2014学年春季学期课程代号：0907322课程名称：工程常用算法授课教师：郭健副教授小组成员班级学号姓名成绩20110711 201107122 姚兴华20110711 201107124 张超凡2014年 4 月 9 日《工程常用算法》综合实践作业一非线性方程求根：(45组) 2014年3月9日班级 学号 姓名 主要工作说明 自评成绩07112011071122 姚兴华 编程&排版 优0711 2011071124 张超凡 流程图&总结已知方程x e x x y -+=232，请分别用二分法、牛顿法和割线法求此方程的根。

误 差限取：1210-=ε。

注意先确定出方程的有根区间。

二．程序流程图开始结束|x2-x1|<ex? 或|f(x2)|<ey?x2 = x1-f(x1)*(x1-x0)/(f(x1)-f(x0))x2值x0 = x1x1 = x2割线法2x0,x1,ex,eyny开始fd(x0)=0?结束|x1-x0|<ex?或|f(x1)|<ey?x1 = x0-f(x0)/fd(x0)x1值x0 <= x1牛顿迭代法：f(x)函数fd(x)一阶导函数无根输入：x0,ex,eynnnyy三．完整的程序及简要的注释 1.二分法# include <stdio.h> # include <math.h> #define ex 1e-12#define ey 1e-12 /* 宏定义ex 和ey 值 */ int n1=0,n2=0,n3=0; /* 计算次数*/ /* f(x)函数 */ double f(double x) {return x*x*x+2*x*x-exp(x); }/* 牛顿迭代法 */ double newton(double x1) {double x0;printf("牛顿法求根:\n"); if(f(x1)==0) printf("无根");else{do{n1++;x0=x1;x1=(2*pow(x0,3)+2*x0*x0+(1-x0)*exp(x0))/(3*x0*x0+ 4*x0-exp(x0));printf("n1=%d x1=%.6f f(x1)=%.6f\n",n1,x1,f(x1)); }while(fabs(x0-x1)>ex&fabs(f(x1))>ey);}printf("牛顿法计算次数n1=%d\n\n",n1);}/* 二分法*/double BiPartition(double a,double b){printf("二分法求根:\n");double x;if(f(a)*f(b)>0) printf("无根");else{do{x=(a+b)/2;if(f(a)*f(x)<0) b=x;else a=x;n2++;printf("n2=%d x=%.6f f(x)=%.6f\n",n2,x,f(x));}while(fabs(b-a)>ex&fabs(f(x))>ey);}printf("二分法计算次数n2=%d\n\n",n2);}/* 割线法*/double gexian(double x1,double x2){double x0;printf("割线法:\n");do{n3++;x0=x1;x1=x2;x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));printf("n3=%d x2=%.6ff(x2)=%.6f\n",n3,x2,f(x2));}while(fabs(x2-x1)>ex&fabs(f(x2))>ey);printf("割线法计算次数n3=%d\n\n",n3);}/* 将三种方法的计算次数按从小到大排序*/void sort(){printf("计算次数排序:\n");if(n1<n2&n2<n3) printf("牛顿法<二分法<割线法\n"); if(n1<n3&n3<n2) printf("牛顿法<割线法<二分法\n"); if(n2<n1&n1<n3) printf("二分法<牛顿法<割线法\n"); if(n3<n2&n2<n1) printf("割线法<二分法<牛顿法\n"); if(n2<n3&n3<n1) printf("二分法<割线法<牛顿法\n"); if(n3<n1&n1<n2) printf("割线法<牛顿法<二分法\n"); }/* 在主函数中运行各种求解方法*/int main(){newton (1);BiPartition(0,1);gexian(0,1);sort();}四．程序运行结果五、对不同实现方法的运行结果进行比较(10%)【从收敛速度，精度比较即可】解：牛顿法计算次数为4次，割线法计算次数为7次，二分法计算次数为39次，可知牛顿法收敛速度最快，割线法次之，二分法收敛次数最慢。

第12章三次样条众所周知，使用高阶多项式的插值常常产生病态的结果。

目前，有多种消除病态的方法。

在这些方法中，三次样条是最常用的一种。

在MATLAB中，实现基本的三次样条插值的函数有spline，ppval，mkpp和unmkpp。

在这些函数中，仅spline在《MATLAB参考指南》中有说明。

下面几节，将展示在M文件函数中实现三次样条的基本特征。

12.1 基本特征在三次样条中，要寻找三次多项式，以逼近每对数据点间的曲线。

在样条术语中，这些数据点称之为断点。

因为，两点只能决定一条直线，而在两点间的曲线可用无限多的三次多项式近似。

因此，为使结果具有唯一性。

在三次样条中，增加了三次多项式的约束条件。

通过限定每个三次多项式的一阶和二阶导数，使其在断点处相等，就可以较好地确定所有内部三次多项式。

此外，近似多项式通过这些断点的斜率和曲率是连续的。

然而，第一个和最后一个三次多项式在第一个和最后一个断点以外，没有伴随多项式。

因此必须通过其它方法确定其余的约束。

最常用的方法，也是函数spline所采用的方法，就是采用非扭结(not-a-knot)条件。

这个条件强迫第一个和第二个三次多项式的三阶导数相等。

对最后一个和倒数第二个三次多项式也做同样地处理。

基于上述描述，人们可能猜想到，寻找三次样条多项式需要求解大量的线性方程。

实际上，给定N个断点，就要寻找N-1个三次多项式，每个多项式有4个未知系数。

这样，所求解的方程组包含有4*(N-1)个未知数。

把每个三次多项式列成特殊形式，并且运用各种约束，通过求解N个具有N个未知系数的方程组，就能确定三次多项式。

这样，如果有50个断点，就有50个具有50个未知系数的方程组。

幸好，用稀疏矩阵，这些方程式能够简明地列出并求解，这就是函数spline所使用的计算未知系数的方法。

12.2 分段多项式在最简单的用法中，spline获取数据x和y以及期望值xi，寻找拟合x和y的三次样条内插多项式，然后，计算这些多项式，对每个xi的值，寻找相应的yi。

B班大作业要求：1. 使用统一封皮；2. 上交大作业内容包含：一摘要二数学原理三程序设计（必须对输入变量、输出变量进行说明；编程无语言要求，但程序要求通过）四结果分析和讨论五完成题目的体会与收获3. 提交大作业的时间：本学期最后一次课，或考前答疑；过期不计入成绩；4. 提交方式：打印版一份；或手写大作业，但必须使用A4纸。

5. 撰写的程序需打印出来作为附录。

课程设计课程名称：设计题目：学号：姓名：完成时间：题目一：非线性方程求根 一 摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。

本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目，分析并总结不同方法处理问题的优缺点。

观察迭代次数，收敛速度及初值选取对迭代的影响。

用Newton 法计算下列方程（1） 310x x --= ， 初值分别为01x =，00.45x =，00.65x =；（2） 32943892940x x x +-+= 其三个根分别为1,3,98-。

当选择初值02x =时给出结果并分析现象,当6510ε-=⨯，迭代停止。

二 数学原理对于方程f(x)=0，如果f(x)是线性函数，则它的求根是很容易的。

牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。

设已知方程f(x)=0有近似根x k (假定k f'(x )0≠) ，将函数f(x)在点x k 进行泰勒展开，有k k k f(x)f(x )+f'(x )(x-x )+≈⋅⋅⋅于是方程f(x)=0可近似的表示为k k k f(x )+f'(x )(x-x )=0这是个线性方程，记其根为x k+1，则x k+1的计算公式为k+1k ()x =x -'()k k f x f x ，k=0,1,2,… 这就是牛顿迭代法或简称牛顿法。

三 程序设计（本程序由Fortran 语言编制）（1）对于310x x --=，按照上述数学原理，编制的程序如下program newton implicit nonereal :: x(0:50),fx(0:50),f1x(0:50)!分别为自变量x ，函数f(x)和一阶导数f1(x) integer :: kwrite(*,*) "x(0)="read(*,*) x(0) !输入变量：初始值x(0)open(10,file='1.txt') do k=1,50,1fx(k)=x(k-1)**3-x(k-1)-1 f1x(k)=3*x(k-1)**2-1x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k 及x 的值 write(10,'(I3,1x,f11.6)') k,x(k)if(abs(x(k)-x(k-1))<1e-6) exit !终止迭代条件 end do stop end（2）对于32943892940x x x +-+=，按照上述数学原理，编制的程序如下program newton implicit nonereal :: x(0:50),fx(0:50),f1x(0:50)!分别为自变量x ，函数f(x)和一阶导数f1(x) integer :: kwrite(*,*) "x(0)="read(*,*) x(0) !输入变量：初始值x(0)open(10,file='1.txt') do k=1,50,1fx(k)=x(k-1)**3+94*x(k-1)**2-389*x(k-1)+294 f1x(k)=3*x(k-1)**2+188*x(k-1)-389x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k 及x 的值 write(10,'(I3,1x,f11.6)') k,x(k)if(abs(x(k)-x(k-1))<5e-6) exit !终止迭代条件 end do stop end四 结果分析和讨论（1）对于方程310x x --=，当初值分别为01x =，00.45x =，00.65x =时，所得结果如下结果分析与讨论：从计算结果可以看到，当取的初始值不同时，虽然均得到了近似解x *=1.324718，但收敛速度明显不同。

5.7 三次样条函数在制造船体和汽车外形等工艺中传统的设计方法是，首先由设计人员按外形要求，给出外形曲线的一组离散点值，施工人员准备好有弹性的样条（一般用竹条或有弹性的钢条）和压铁，将压铁放在点的位置上，调整竹条的形状，使其自然光滑，这时竹条表示一条插值曲线，我们称为样条函数。

从数学上看，这一条近似于分段的三次多项式，在节点处具有一阶和二阶连续微商。

样条函数的主要优点是它的光滑程度较高，保证了插值函数二阶导数的连续性，对于三阶导数的间断，人类的眼睛已难以辨认了。

样条函数是一种隐式格式，最后需要解一个方程组，它的工作量大于多项式拉格朗日型式或牛顿型式等显式插值方法。

定义给定区间上个节点和这些点上的函数值，。

若满足；在每个小区间上至多是一个三次多项式；在上有连续的二阶导数，则称为关于剖分的三次样条插值函数，称为样条节点。

要在每个子区间上构造三次多项式，共需要个条件，由插值条件，提供了个条件；用每个内点的关系建立条件又得到个条件；再附加两个边界条件，即可惟一确定样条函数了。

用待定系数法确定了构造样条函数的存在性和惟一性。

在具体构造样条函数时一般都不使用计算量大的待定系数法。

下面给出构造三次样插值的关系式和关系式的方法。

5.7.1 三次样条插值的M关系式引入记号。

用节点处二阶导数表示样条插值函数时称为大关系式，用一阶导数表示样条插值函数时称为小关系式。

问题5.8给定插值点，怎样构造用二阶导数表示的样条插值函数，即怎样构造关系式？假设。

由于在上为线性函数，故在上做的分段线性插值函数：令，得到（5.29）对积分两次有（5.30）将代入式（5.27）可解出故在上有（5.31）在每个小区间上具有不同的表达式，但由于在整个区间上是二阶光滑的，故有列出每一个关系式，再经计算得：（5.32）其中：由式（5.32）得到个未知数的个方程组。

现补充两个边界条件，使方程组只有惟一解。

下面分三种情况讨论边界条件。

（1）给定的值时，称为自然边界条件），此时阶方程组有个未知量，即（5.33）（2）给定的值，它们分别代入在中的表达式，得到另外两个方程：于是需要解阶的方程组：（5.34）（3）被插函数以为基本周期时，即，即；即。

计算方法大作业选题（任选一题）1.插值方法：编制用牛顿插值、哈密特插值、分段插值、样条插值的计算机程序，并就计算结果分 析它们的特点。

2. 数值积分方法：编制用牛顿一科特斯、复化求积、龙贝格公式计算积分的计算机程序，并就计算结果 分析它们的特点。

3. 求解线性方程组：编制用直接法（消去法、三角分解法）与间接法（迭代法）解线性方程组的计算机 程序，并就计算结果分析它们的特点。

4. 用高斯-勒让德公式计算定积分设计用高斯-勒让德求积公式计算定积分的程序，并用数值例子计算。

5. 求大区间积分的数值方法2000 2 设计一种方法求解积分I 二 e 」dx 的数值方法，并分析它的可行性6. 椭圆数值积分2 -IT I ------------- -- ------- -- ---已知椭圆的周长可以表示成 s=a,「｛4 + P cos Q d （0 v P v 1），取a=1，（1）针对「从0.1到0.9 （步长h=0.1）分别求出周长s ；（用Romberg 积分方法）（2）对于以上数据，求出s 关于「的插值多项式;对于（1）中数据，试用最小二乘的思想求作拟合多项式（要求是偶次）优劣进行比较（1）试采用Romberg 求积算法计算当x = 0.5、1、1.5、2、2.5、3、3.5时的门（x ）的值;（2） 针对（1）中求出的叮」（x ）值，选择适当的曲线做曲线拟合。

备选曲线：①多项式曲线 ②y x③y ■夫（选择其中一种即可） Ax+B 1+Ce8. 曲线拟合1601年，德国天文学家开普勒发表了行星运行第三定律： T 二Cx 3/2，其中，T 为行星绕太阳 旋转一周的时间（单位：天），x 表示行星到太阳的平均距离（单位：百万公里），并测得水 星、金星、地球、火星的数据（x ，T ）分别为（58，88）、（108，225）、（150，365）、（228,687）。

（1） 用最小二乘法估计C 的值；（2） 分别作出上述数据点的直线、抛物线、三次、四次多项式拟合，求出残差平方和 Q ，，并对这些多项式的7•正态曲线的拟合标准正态分布的分布函数 G丄2 e 2dt并比较优劣；(3)用函数y =ae x bx c 来对数据点进行曲线拟合，并求出残差平方和Q 9. 求方程实根用二分法和牛顿迭代法(包括弦截法)编程求方程2 xs i nx 02 的实根，要求误差不超过10・。

《工程常用算法》综合实践作业三151.6飞机下轮廓线形状大致如下图所示：要求分别用拉格朗日插值法、Newton插值法、分段线性插值法和三次样条插值法计算x每改变0.5时y的值，即x 取 0.5, 1, 1.5, … , 14.5 时对应的y值。

比较采用不同方法的计算工作量、计算结果和优缺点。

二、程序流程图三、完整的程序及简要的注释拉格朗日#include<stdio.h>#include<math.h>float x[10]={0,3,5,7,9,11,12,13,14,15},y[10]={0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1. 6},n[29]; /*定义29个隔值 */void X() /*定义隔值函数*/{ int i;for(i=0;i<29;i++) /*隔值循环*/n[i]=(i+1)*0.5;}float l(float xx,int i) /*求Lagrange基函数*/{ int j;float l=1;for(j=0;j<10;j++)if(i!=j) l=l*(xx-x[j])/(x[i]-x[j]);return l;}float Lagrange(float x) /*Lagrange插值法*/{ int i;float f=0;for(i=0;i<10;i++)f=f+y[i]*l(x,i);return f;}/*调用Lagrange插值函数*/int main(){ int i;float Y[29];X();printf("Lagrange插值法:\n");for(i=0;i<29;i++){ float Y[29];Y[i]=Lagrange(n[i]);printf("n[%d]=%.2fY[%d]=%8.4f\n",i,n[i],i,Y[i]);}}Newton #include<stdio.h>#include<math.h>float x[10]={0,3,5,7,9,11,12,13,14,15},y[10]={0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1. 6},n[29]; /*定义29个隔值 */floatf[10][10]={0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1. 0,1.6};void X() /*定义隔值函数*/{ int i;for(i=0;i<29;i++) /*隔值循环*/n[i]=(i+1)*0.5;}void CS() /*定义差商函数*/{ int i,j;for(i=1;i<10;i++)for(j=i;j<10;j++)f[i][j]=(f[i-1][j]-f[i-1][j-1])/(x[j]-x[j-i ]);}float Newton(float xx) /*Newton插值法*/ { int i,j;float t=1;float N=f[0][0];for(i=1;i<10;i++){ t=t*(xx-x[i-1]);N=N+f[i][i]*t;}return N;}/*调用Newton插值函数*/int main(){ int i;float Y[29];X();CS();printf("Newton插值法:\n");for(i=0;i<29;i++){ float Y[29];Y[i]=Newton(n[i]);printf("n[%d]=%.2fY[%d]=%8.4f\n",i,n[i],i,Y[i]);}}三次样条插值法#include<stdio.h>#include<math.h>float x[10]={0,3,5,7,9,11,12,13,14,15},y[10]={0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1. 6},G[29]; /*定义29个隔值 */void X() /*定义隔值函数*/{ int i;for(i=0;i<29;i++) /*隔值循环*/G[i]=(i+1)*0.5;}float S(float xx) /*三次样条插值法*/{ float h[9],q[9],w[9],C[9],m[10],b[9],c[8],a[9],n[ 9],p[8],r[9],yp[9],s[9];int i;for(i=0;i<9;i++)h[i]=x[i+1]-x[i]; /*计算插值点之间的距离*/for(i=1;i<9;i++){ q[i]=h[i-1]/(h[i]+h[i-1]);w[i]=1-q[i]; /*计算方程组中的两个系数*/C[i]=6/(h[i]+h[i-1])*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]); /*计算方程组的结果*/}m[0]=m[9]=0; /*求线性方程组得m值*/ for(i=1;i<9;i++){ b[i]=2;if(i<8)c[i]=w[i];if(i>1)a[i]=q[i];if(i==1){n[1]=b[1];p[1]=c[1]/n[1];}if(i>1){ r[i]=a[i];n[i]=b[i]-r[i]*p[i-1];if(i<8) p[i]=c[i]/n[i];}}for(i=1;i<9;i++){ if(i==1) yp[1]=C[1]/n[1];else yp[i]=(C[i]-r[i]*yp[i-1])/n[i];}for(i=8;i>0;i--){ if(i==8) m[8]=yp[8];else m[i]=yp[i]-p[i]*m[i+1];}for(i=0;i<9;i++)if(xx>=x[i]&&xx<x[i+1]) { s[i]=pow(xx-x[i+1],3)*m[i]/6/h[i]+pow(xx -x[i],3)*m[i+1]/6/h[i]-(y[i]-h[i]*h[i]*m[i]/6)*(xx-x[i+1])/h[i]+( y[i+1]-h[i]*h[i]*m[i+1]/6)*(xx-x[i])/h[i]; return s[i];}}int main(){ int i;float Y[29];X();printf("三次样条插值法:\n");for(i=0;i<29;i++){ float Y[29];Y[i]=S(G[i]);printf("G[%d]=%.2fY[%d]=%8.4f\n",i,G[i],i,Y[i]);}}分段线性插值#include<stdio.h>#include<math.h>float x[10]={0,3,5,7,9,11,12,13,14,15},y[10]={0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1. 6},G[29]; /*定义29个隔值 */void X() /*定义隔值函数*/{ int i;for(i=0;i<29;i++) /*隔值循环*/G[i]=(i+1)*0.5;}float FD(float xx) /*定义分段线性插值函数*/{ int i;float p;for(i=0;i<10;i++)if((xx>=x[i])&&(xx<x[i+1])){ p=(xx-x[i+1])*y[i]/(x[i]-x[i+1])+(xx-x[i] )*y[i+1]/(x[i+1]-x[i]);return p;}}int main() /*调用分段线性插值函数*/{ int i;X();printf("分段线性插值法:\n");for(i=0;i<29;i++){ float Y[29];G[i]=(i+1)*0.5;Y[i]=FD(G[i]);printf("G[%d]=%.2fY[%d]=%8.4f\n",i,G[i],i,Y[i]);}}四、程序运行结果五、对不同实现方法的运行结果进行比较六、问题与总结一、计算公式及计算方法。