本题主要目的是建立相关模型解决在修建水渠过程中的诸多问题,从而实现工程量最优化。 针对问题一,为了求得开掘水渠的土石方量,本文通过对比分段三次Hermite 插值与三次样条插值,最终采用分段三次 Hermite 插值的方法对已知数 据点进行插值拟合,得到关于水渠的曲线方程 y f x ,对水渠曲线方程积分即得到水渠长度 14550 7650 21 dx yL ,利用 MATLAB 求解得到水渠长度为: m 。因此最终解得开掘水渠的总土石方量为: 3 m135405LSV 。 针对问题二,在问题一的基础上,本文建立积分上限函数模型:令 7650 a ,1x 满足 1 2 16xaVS y dt ,求出 i x 后, 1 ix 满足 1 2 16iixxVS y dt ,从而将总土石方量的六等分,得到 7 .8736x 1 , 2 .9862x 2 , 10956 x3 , 12116 x4 ,13353x 5 ,进而确定了六等分点的坐标 y,x 。 针对问题三,设在沿水渠的公路上有三个变量,分别为 k ji x ,x,x ,为使得运输工作量最小,本文建立了无约束规划模型,利用 MATLAB 求解得到最小运输量为 4 m 。并给出了修建两条公路时水渠上的位置坐标 和 , 。 关键词:Hermite 插值 MATLAB 积分上限函数 无约束规划 一、问题重述 在某地区开掘水渠,已知该水渠经过的若干点。 问题一,求解水渠施工的总石方量; 问题二,如果将水渠的分成 6 截,每截土石方量相同,分段点应该取在何位置; 问题三,设平行于水渠修一条路。河道中挖出的土石方要运往 A(9500,4000)处为了方便运输,计划在沿水渠的公路上选择两点修建通往 A 处的临时公路,使得总的土石方运输工作量最小。 二、问题的分析 针对问题一,本题要求开掘水渠的总土石方量,已知水渠截面积,则主要目的在于求得水渠长度。已知水渠经过的若干点的位置,要得到水渠的长度,本文想到用插值拟合可以得到水渠曲线,对曲线积分则得到水渠长度。插值与拟合的方法有多种,样条插值会较光滑,但不一定能保持原有形状,考虑到要更好的保持水渠的形状,于是,本文选用 Hermite 方法进行插值拟合。 针对问题二, 要将水渠六等分且每段的土石方量相同,此问题为函数的反解问题,因此,在已知水渠曲线函数的情况下,本文可以考虑到用积分上限函数求解,从而确定 x 点,进而得到 y 点。 针对问题三,要修建公路以运输土石方,从而使运输量工作量最小。此问题为规划问题,在问题二中,本文已知 x 与土石方量 V 存在关系,又因为运输工作量等于土石方量与距离的乘积,因此,本文使用无约束规划模型,求工作量最小值即可。 三、模型假设 1、修建的两条临时公路为直线。 2、沿水渠的公路函数曲线近似与水渠的曲线函数相同。 四、符号说明 xf 水渠曲线方程 V 土石方量 S 水渠截面积 L 水渠长度 ix 水渠上点的横坐标 iy 水渠上点的纵坐标 iW 土石方运输工作量 1L 临时公路 2L 临时公路 五、模型的建立与求解 问题一 插值与拟合 由已知水渠经过的点,做出散点图(图 1) 1 x 10 420002500300035004000450050005500600065007000X/mY/m水渠散点图 图 1.水渠散点图 方法 1、利用 Hermite 方法对已知数据点进行插值。 【3】 设 已 知 函 数 xfy 在 1 n 个互异节点 n 10 x ,L,x,x 上 的 函 数 值 ii xfy n,L,1,0i 和导数值 i ' i ' x fy ,要求一个至多 2 n +1 次的多项 式 xH ,使得 i i yxH i ' i ' y xH n,1,0i Hermite 插值多项式为: 2 ' i i i i i i H x h x x a y y y 其中,2nij 0j j ij i x x xx h , n ij 0j j i i x x 1 a 。 利用 MATLAB 进行插值,得到插值曲线(图 2)。 1 x 10 420002500300035004000450050005500600065007000Hermite插值曲线与原始数据点X/mY/m Hermite插值曲线 原始数据点 图 插值曲线与原始数据点 方法 2、利用样条差值对已知数据点进行插值。 【3】 定义样条函数: 数学上将具有一定光滑性的分段多项式称为样条函数。具体的说,给定区间 a,b 的一个划分 0 1 1nn :a x x x x b 如果函数 () sx满足: 1. 在每个小区间 1 , ( 0,1, , 1) ii x x i n 上 () sx是k 次多项式; 2. () sx在 a,b 上具有 1 k 阶连续导数。 则称 () sx为关于划分的k 次样条函数,其图形称为k 次样条曲线。 01 , , , n x x x 称为样条节点, 1 2 1 , , , n x x x 称为内节点, 0, n xx称为边界点,这样样条函数的全体记作 ( , ) p Sk ,称为k 次样条函数空间。 显然,折线是一次样条曲线。 若 ( ) ( , ) p s x S k ,则 () sx是关于分划的k 次多项式样条函数。k 次多项式样条函数的一般形式为 101 ( ) ( ) !! i kn j ki kj ij x s x x x ik 其中 ( 0,1, , ) i ik 和 ( 1,2, , 1) j jn 均为任意常数,而 ( ) , ( ) , ( 1,2, , 1) 0, k jjk j j x x x x x x j n xx 本文使用 3 k 的情况:即为三次样条函数。 三次样条函数:对于 a,b 上的划分 0 1 1nn :a x x x x b ,则 1 2 3 3 323 0 11 ( ) ( ) ( ,3) 2! 3! 3! n j jp j aa s x x x x x x S 其中3 3 ( ) , ( ) , ( 1,2, , 1) 0, jj j j x x x x x x j n xx 三次样条函数差值: 由于 3( ) ( ,3) ps x S中含有 3 n 个待定系数,故应需要 3 n 个插值条件,已知插值节点 i x 和相应的函数值 ( ) ( 0,1,2, , ) ii f x y i n ,这里提供了 1 n 个条件,还需要 2 个边界条件。 常用的三次样条函数的边界条件有 3 中类型: (1) 3 0 3 ( ) , ( ) n s a y s b y 。由这中边界条件建立的样条插值函数称为 () fx的 完备三次样条插值函数。 特别的, 0'0 n yy 时,样条曲线在端点处呈水平状态。 如果 () fx 不知道,可以要求 3() sx 与 () fx 在端点处近似相等。这时以0 1 2 3 , , , x x x x 为节点作一个三次 Newton 插值多项式 () a Nx,以 1 2 3 , , , n n n n x x x x 作一个三次 Newton 插值多项式 () b Nx,要求 ( ) ( ), ( ) ( ) ab s a N a s b N b 由这种边界条件建立的三次样条称为 () fx的 Lagrange 三次样条插值函数。 (2) 3 0 3 3 ( ) , ( ) s a y s b y 。特别的 0 nn yy 时,称为自然边界条件。 (3) 3 3 3 3 ( 0) ( 0), ( 0) ( 0) s a s b s a s b ,(这里要求 33 ( 0) ( 0) s a s b )此条件称为周期条件。 利用 MATLAB 进行三次样条插值,得到插值曲线(图 3)。 1 x 10 420002500300035004000450050005500600065007000X/mY/m三次样条插值曲线与原始数据点 三次样条插值曲线 原始数据点 图 3.三次样条插值曲线与原始数据点 Hermite 插值与三次样条插值的对比【5】: SPLINE 提供的函数 s(x)的构建方法和 PCHIP 里面的函数 p(x)完全相同,只 O x y 0 AM 1M2M1nM n BM 图 4 不过在 X(j)处的斜率的选择方法不一样, SPLINE 函数的 s(x)在 X(j)的二阶导数 D^2s(x)也是连续的,这导致了如下结果: (1) SPLINE 更加光滑,即,D^2s(x)是连续的。 (2) 如果数据是一个光滑函数的值,则 SPLINE 更加精确。 (3) 如果数据不是光滑的,则 PCHIP 没有 overshoots,也不太震荡。 (4) PCHIP 建立的难度较小。 (5) 这两种函数估计的难度是一样的。 三次样条比 Hermite 插值光滑,样条的两阶导数连续,而 Hermite 插值一阶导数连续。不连续的两阶导数隐含着不连续的曲率。人的眼睛可以检测出图形上曲率的不连续。另一方面,Hermite 插值是保形状的,而样条插值不一定保形状。 通过对比 Hermite 插值与三次样条插值,针对本题并无明显差异。为了更好的保证图形形状,减小误差,本文采用 Hermite 插值。 求解水渠长度 圆的周长可以利用圆的内接正多边形的周长当边数无限增多时的极限确定。类似的方法,可以用来建立平面连续曲线的弧长,应用定积分来计算弧长。 设 AB 、 是曲线弧的两个端点。在弧 AB 上以此取分点: 0 1 2 1 1 , , , , , , , , i i n n A M M M M M M M B ,并以此连接相邻分点得一折线(图4)。 当分点的数目无限增加且每小段 1ii MM 都缩向一点时,如果此折线的长11niiiMM 的极限存在,则称此极限为曲线弧 AB 的弧长,并称此曲线弧 AB 是 可求长的。 由于光滑曲线弧是可求长的,故可应用定积分来计算弧长。 设曲线弧由参数方程: () , () xt t yt 给出,其中 ( ), ( ) tt 在 , 上具有连续导数,且 ( ) ( ) tt 、 不同时为零,现计算该曲线弧的长度。 取参数t为积分变量,它的变化区间为 , 。相应于 , 上任一小区间 , t t dt 的小弧段的长度 s 近似等于对应的弦的长度 22 ( ) ( ) xy ,因为 ( ) ( ) ( ) x t dt t dx t dt ( ) ( ) ( ) y t dt t dy t dt 所以, s 的近似值(弧微分)即弧长元素为 2 2 2 2 2 2 2 2 ( ) ( ) ( )( ) ( )( ) ( ) ( ) ds dx dy t dt t dt t t dt 于是所求弧长为 22 ( ) ( ) s t t dt 当曲线弧由直角坐标方程 ( ) ( ) y f x a x b 给出,其中 () fx在 , ab上具有一阶连续导数,这时曲线弧由参数方程 ()() xxa x by f x 从而所求的弧长为 21b a s y dx 利用插值后得到的水渠的曲线函数,对其进行积分,则为水渠长度。 14550 7650 21 dx yL 用 MATLAB 求解得到 m L 求解土石方量 已知,水渠长度,水渠截面积。 则: 2 m1822810S 3m135405LSV 问题二 设函数 () fx在区间 , ab上连续,并且设x为 , ab上的一点。观察 () fx在部 分区间 , ax上的定积分 ()xa f x dx 首先,由于 () fx在 , ax上依旧连续,因此该定积分存在。这里,x即表示定积分的上限,又表示积分变量。因为定积分与积分变量的记号无关,所以,为了明确起见,可以吧积分变量改用其他符号,例如用t表示,则上面的定积分可以写成 ()xa f t dt 如果上限x在区间 , ab上任意变动,则对于每一个取定的x值,定积分有一 个对应值,所以它在 , ab上定义了一个函数,记作 () x : ( ) ( ) ( ) x a x f t dt a x b () x 便为积分上限函数。 本文针对问题二建立积分上限函数模型: 216xaVS y dt 通过起点a作为积分下限,求得第一个积分上限,即第一个等分点,第一个等分点为积分下限,求得第二个积分上限,即第二个等分点,以此类推。改变积分上下限,确定5个等分点,将水渠六等分,且每段土石方量相同。 利用 MATLAB 求解(见附录 ),得到等分点坐标为: , , , , , 3540,13353 且每段的土石方量为: 3 问题三 由问题二知土石方量 V 与水渠曲线函数存在关系。首先建立 xFV 模型。设在沿水渠的公路上有三个变量为 k ji x ,x,x ,修建的临时公路需要保证运输工作量最小,因此,在 D 点左边开掘水渠的土石方均运到 B 处,在 D 右边开掘水渠的土石方都运往 C 处。最终将土石方由 B、C 两处运往 A 处(示意图见图4) y = f(x)(xk)(xj)D(xi)L2L1ACB 图4.水渠临时公路修建示意图 运输工作量等于土石方量乘以距离,因此对于水渠曲线上的运输工作量本文建立的模型为: 以 i 0 x ~x 段为例,设:该段水渠长度为n L L,L i 0 i0 ,该段土石方量为n SLn V V,V i 0i0 i0 , 则: )LnL(V)L2L(V)LL(VW i0i0i0i0 )n21(LVVLn i0 n2 1nLV 当 i 0 i 0 x x2xx2 i0 dx y1dxy1S 2 1 LV 2 1 W,n 同理可得 1 kjkij W,W,W 则: iji01kjk1 WWWWW 14550 14550 2 2 2 2 11 1 1 1 1 22 kk j j k k xx x x x x S y dx y dx S y dx y dx j i j i ii x x 2x x 2x 7650 2x 7650 2 dx y1dxy1S 2 1 dxy1dxy1S 2 1 对于由 B、C 两点运往 A 处的运输工作量本文建立的模型为: 2 0k 2 0k 14550 x 22 0i 2 0ix76502 2 y yxxdxy1Syyxxdxy1SWjj 要使运输工作量最小,即 1 W 、 2 W 之和达到最小,因此,本文建立无约束规划模型: kji2kji1 x,x,xWx,x,xWMin 即: 2 0i 2 0i x 7650 2 2x x 2 2x 7650 2 y yxxdxy1dxy1 2 1 dxy1 2 1 Min j j i i 2 0k 2 0k 14550 x 2 214550 x 2 2 x x 2 y yxxdxy1dxy1 2 1 dxy1 2 1 jk k j运用 MATLAB 求解得:(程序代码见附录 ) 最小运输工作量为: 4 m B、C 两点的坐标为: 和 , 。 六、模型的优缺点 优点: 1、 通过对比 Hermite 插值与三次样条插值,发现求得的水渠长度分别为 和 ,对本题无明显差异。 2、 对于运输量的规划问题,准确的反应了最优解。 缺点: 1、 对于插值函数的曲线积分,近似了曲线的导数,存在一定误差。 2、 规划问题的运算量较大。利用 MATLAB 算法优势不明显。