数学建模数值算法

数值算法2014.8主要内容插值和拟合方程f(x)=0求根神经网络、遗传算法、模拟退火、模糊控制算法选用算法应遵循的原则1.尽量简化计算步骤，减少乘除运算的次数.

例如，计算多项式通常运算的乘法次数为若采用递推算法,

则乘法次数仅为n.又如2.防止大数“吃掉”小数当|a|>>|b|时，尽量避免a+b。例如，假设计算机只能存放10位尾数的十进制数，则3.尽量避免相近数相减例如，当x很大时，应

，

当x接近于0时，应4.避免绝对值很小的数做分母当|b|<<|a|时，应尽量避免。5.选用数值稳定性好的算法，以控制舍入误差高速增长例如若（误差）则计算

时误差扩大了倍,而

（n=1,2,…）

是稳定的。

构造一个(相对简单的)函数通过全部节点,即再用计算插值，即插值拟合算法一维插值的定义已知n+1个节点其中互不相同，不妨设求任一插值点处的插值节点可视为由产生,g(x)表达式复杂,或无封闭形或未知。

称为拉格朗日插值基函数。

已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn

。求一n次多项式函数Pn(x)，使其满足：

Pn(xi)=yi,i=0,1,…,n.

解决此问题的拉格朗日插值多项式公式如下其中Li(x)为n次多项式：拉格朗日(Lagrange)插值例：选取n+1个不同插值节点，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出下列函数拉格朗日插值插值多项式图形.

分段线性插值计算量与n无关;n越大，误差越小.xjxj-1xj+1x0xnxoy比分段线性插值更光滑。xyxi-1xiab

在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值

三次样条插值g(x)为被插值函数。二维插值的定义xyO第一种（网格节点）：例：样条插值的应用

用三次样条插值描绘如图所示图形的上部轮廓线。

已知mn个节点其中互不相同，不妨设

构造一个二元函数通过全部已知节点,即再用计算插值，即第二种（散乱节点）：yx0已知n个节点其中互不相同，

构造一个二元函数通过全部已知节点,即再用计算插值，即

注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。

将四个插值点（矩形的四个顶点）处的函数值依次简记为：

分片线性插值xy(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足

双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O拟合与插值的关系

函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题；常用的是最小二乘拟合(略)最临近插值、线性插值、样条插值与曲线拟合结果：二、方程f(x)=0求根1.区间迭代法：二分法

2.迭代法

(1)将方程f(x)=0等价变形为x=(x)，建立迭代关系式：

xk+1=(xk)，k=0,1,2,…

若(x)连续，且迭代序列{xk}的极限存在，则该极限就是x*。

(2)局部收敛与全局收敛(3)收敛定理5．牛顿迭代法

(2).牛顿法的收敛性：平方收敛(3)牛顿迭代法的优缺点：收敛速度快，需要求导数，局部收敛。(4)牛顿下山法:取适当参数(0,1]，使成立k=0,1,2,…。

(1).牛顿迭代公式6．弦截法

(1)弦截法的公式k=0,1,2,…

(2)快速弦截法

k=0,1,2,…

(3)弦截法的收敛性例.快速求解空间任意螺旋线与任意平面的全部交点。

圆柱螺旋线的向量式：(x,y,z)T=(cos(t-t0),sin(t-t0),z)T

经坐标轴旋转-平移得空间一般螺旋线方程平面的一般方程

ax+by+cz=d把螺旋线方程代入平面方程整理，用t/代替原来的t，cost、sint、t的系数分别为A、B、C常数项为D,得：Acost+Bsint+Ct+D=0

令：f(t)=Acost+Bsint+Ct+D问题归为解非线性方程f(t)=0。

分析方程左边函数f(t)，其具有如下的特点：

(1)f(t)极值的出现具有周期2(2)f(t)=0的大部分根在t轴上下方两相邻极值点之间，个别在极值点处。(3)令f

(t)=0，由得f(t)极值点位于

(4)若相邻两极值f(t1)f(t2)<0，则必有根t

(t1，t2)；过点A1(t1，f(t1))、A2(t2，f(t2))作直线A1A2，以A1A2与t轴的交点t0为初值，一般地可用牛顿法求解，特别地当根t位于极值点附近时，牛顿法收敛速度太慢，可取区间[t1，t2]，用二分法求解。

(5)特殊情况

1.当C=0时，f(t)是周期函数，振幅不超过，当时，f(t)=0无解。当时，其解周期性出现f(t)=0有唯一解或无解。如果有解一定在弯曲点处。

算法描述输入平面方程一般式，螺旋线方程一般参数式，用3个欧拉角和映像z轴到螺旋线中心轴的变换向量把输入平面方程一般式，螺旋线方程参数式判断方程是否奇异求极值点有解否？输出无解标志用二分法求唯一解求周期解解在极值附近否？用二分法求解用牛顿法求解是唯一解?yyyNNNN

二分法有根区间的确定牛顿法初值的确定：(牛顿法局部收敛，而且会出现加收敛)以连接相邻两个极值点的直线与t轴交点横坐标为初始点。数值积分和数值微分

1．黎曼和：

2．牛顿—柯特斯(Newton—Cotes)求积公式

将积分区间[a,b]等分n份,记,分点为,k=0,1,...,n,用拉格朗日插值多项式代替被积函数f(x)，得到n阶Newton—Cotes公式3.常用的Newton—Cotes公式n=1时，得梯形公式:In=2时，得辛普森公式:I4复化求积公式随着n的增加可减少积分误差，但高阶Newton—Cotes公式又会造成数值不稳定，因而采用复化公式。将积分区间[a,b]等分n等份,在每个小区间[xk,xk+1]上使用低阶Newton—Cotes公式,(k=0,1,...,n).1).复化梯形公式

在每个子区间上用梯形公式,将它们累加得复化梯形公式：2).复化辛浦生公式记xk+1/2是区间的中点,在每个子区间上用辛浦生公式,并将它们累加得复化辛浦生公式：5.变步长梯形法(逐次分半梯形公式)以复化梯形公式为例n等分区间2n等分区间近似有：类似，复化Simpsom公式6.先看看事后误差估计（不同的误差表达式，事后误差估计式是不同的）自适应计算记为复化一次，2次的Simpson公式控制求7．自适应步长法

函数变化有急有缓，为了照顾变化剧烈部分的误差，需要加密格点。对于变化缓慢的部分，加密格点会造成计算的浪费。以此我们介绍一种算法，可以自动在变化剧烈的地方加密格点计算，而变化缓慢的地方，则取稀疏的格点。是8．龙贝格(Romberg)算法

受事后误差估计式启发，用低阶的公式线性组合后成为一个高阶的公式。

Romberg

算法：<?<?<?………………

T1=)0(0T

T8=)3(0TT4=)2(0T

T2=)1(0T

S1=)0(1T

R1=)0(3T

S2=)1(1T

C1=)0(2T

C2=)1(2T

S4=)2(1T

直到|Tm(0)-Tm-1(0)|<ε或

二重积分的复化梯形公式

系数，在积分区域的四个角点为1/4，4个边界为1/2，内部节点为1定义：若积分的数值积分公式对于任意一个次数不高于m次的多项式都精确成立，而存在一个m+1次多项式使之不精确成立，则称该数值积分公式的代数精确度为m。

9．积分公式的代数精确度

引理:n阶Newton—Cotes公式的代数精确度至少是n。

当n为奇数时，n阶Newton—Cotes公式的代数精确度为n；当n为偶数时，n阶Newton—Cotes公式的代数精确度是n+1

10．高斯型求积公式

适当选取其中2n+2个参数Ak,xk(k=0,1,…,n),使得数值积分公式的代数精确度达到2n+1,这类求积公式称为高斯型求积公式。

对于求积公式:(2)求出pn(x)的n个零点x1,x2,…xn即为Gsuss点.

(1)求出区间[a,b]上权函数为W(x)的正交多项式pn(x).(3)计算积分系数

Gauss型求积公式的构造方法(3)复化高斯求积公式(4)高精度Lobatto数值积分Lobatto积分公式是高斯型求积公式的修正，始终将上下端点作为节点，n+1阶Lobatto积分公式的代数精度达到2n+1

(1)Gauss-Laguerre求积公式(2)Gauss-Hermite求积公式

将积分区间[a,b]分成n等分，在每个子区间上用两点高斯求积公式，再把他们累加得区间[a,b]上复化Gauss-Legendre求积公式。

11.利用样条函数计算列表函数的数值积分：先构造样条函数，然后再求积分。

1、插值型求导公式已知f(x)在n+1个点上的函数值，构造插值多项式Pn(x)近似代替原函数f(x)，利用P’n(x)来代替f’(x)。这就称为插值型求导公式

常微分方程的数值解法考虑一阶常微分方程初值问题:

常微分方程的数值解：设微分方程的解y(x)的存在区间是[a,b]，在[a,b]内取一系列节点a=x0<x1<…<xn=b,其中hk=xk+1-xk

；(一般采用:h=(b-a)/n)。在每个节点xk求解函数y(x)的近似值:yk≈y(xk)，这样y0,y1,...,yn称为微分方程的数值解。用数值方法，求得f(xk)的近似值yk，再用插值或拟合方法就求得y(x)的近似函数。(二)、主要算法1.显式欧拉格式：k=1,2,n-1.

隐式欧拉格式：2.两步欧拉格式：k=1,2,n-1.3.梯形方法:；3.改进的欧拉方法：预测值:

校正值:

2．局部截断误差局部截断误差：当y(xk)是精确解时，由y(xk)按照数值方法计算出来的的误差y(xk+1)-称为局部截断误差。注意：yk+1和的区别。因而局部截断误差与误差ek+1=y(xk+1)-yk+1不同。如果局部截断误差是O(hp+1),我们就说该数值方法具有p阶精度。

5、龙格-库塔法龙格-库塔法的思想：选取[xk,xk+1]上若干个点的导数，将它们作线性组合作为平均斜率，使数值方法达到预期的阶。

1).二阶龙格-库塔法：当:时，得一簇龙格-库塔公式，其局部截断误差均为O(h3)都是二阶精度。特别取

就是改进欧拉公式。2)、经典四阶龙格－库塔格式(也称为标准龙格－库塔方法)：

四阶龙格－库塔方法的截断误差为O(h5)，具有四阶精度。一般一阶常微分方程初值问题用四阶龙格－库塔方法计算，其精度均满足实际问题精度要求。

3).变步长龙格－库塔方法：从节点xk出发，以步长h据四阶龙格－库塔方法求出一个近似值,然后以步长h/2求出一个近似值，若

>,继续将步长折半，若

,将步长加倍，直到>，这时再步长折半一次，即得结果()。

4).RKF格式：变步长龙格－库塔方法，因频繁加倍或折半步长会浪费计算量。Felhberg改进了传统龙格－库塔方法，得到RKF格式，较好解决了步长的确定，而且提高了精度与稳定性，为Matlab等许多数值计算软件采用。4/5阶RKF格式：由4阶龙格－库塔方法与5阶龙格－库塔方法结合而成。

6.亚当斯方法(Adams):使用前面已算出的数值解一边减少计算量。（记：；）1).显式Adams方法：二阶显式Adams方法：;

三阶显式Adams方法：;

四阶显式Adams方法：.2).隐式Adams方法二阶隐式Adams方法：;

三阶隐式Adams方法：;

四阶隐式Adams方法：

3）.Adams预报-校正系统：先用显式格