Transcript 插值与拟合实验
数学实验 Experiments in Mathematics 插 值 与 拟 合 重庆邮电学院基础数学教学部 实验目的 1、掌握用Matlab计算拉格朗日、分段线性、 三次样条三种插值的方法，改变节点的数目， 对三种插值结果进行初步分析。. 2、掌握用Matlab作线性最小二乘的方法. 3、通过实例学习如何使用插值方法与拟合方 法解决实际问题，注意二者的区别和联系 实验内容 1、插 2、拟 值. 合. 3、 数学建模实例 实验软件 MATLAB 插 值 （一） 插值问题的提法 （二）解决插值问题的基本方法 数学建模实例 1、船在该海域会搁浅吗 2、薄膜渗透率的测定 插值问题的提法： 已知n 1个节点( x j , y j ), j 0,1,, n ，其中x j 互不 相同，不妨设 a x0 x1 xn b ，求任一插 x ( x ) 值点 j 处的插值y ，( x j , y j ) 可以看成由某 个函数 y g (x) 产生的，g 的解析表达式可能十分 复杂，或不存在封闭形式，也可以未知。 求解的基本思路： 构造一个相对简单的函数 y f (x) ，使 f 通过全部节点， 即 f ( x j ) y j ( j 0,1,, n) ，再用f (x) 计算插值，即 y f (x ) 拉格朗日多项式插值 从理论和计算角度看，多项式是最简单的函数，设 f(x) 是 n 次多项式，记作 Ln ( x ) an x n an 1x n 1 a1x a0 对于节点( x j , y j ) 应有 Ln ( x j ) y j (1) j 0,1,, n (2) 为了确定插值多项式Ln (x) 中的系数an , an 1,, a0 将（1）代入（2） ，有 an x0n an 1x0n 1 a1x0 a0 y0 n 1 a x n a x a1xn a0 y0 n 1 n n n (3) , 记 x0n x0n 1 1 X xnn xnn 1 1 A ( an , an 1 ,..., a0 )T ,Y ( y0 , y1 ,..., y n )T ， 方程组（3）简写作 XA Y （4） 其中 detX 是 Vandermonde 行列式，利用行列式性质可得 det X ( xk x j ) 0 j k n det X 0 因 x j 互不相同，故 ,于是方程（4）中有唯一解， 即根据 n+1 个节点可以确定唯一得 n 次插值多项式。 实际上比较方便的作法不是解方程（4）求 A，而是先构 造一组基函数： li ( x ) ( x x0 )...(x xi 1 )( x xi 1 )...(x xn ) , ( xi x0 )...(xi xi 1 )( xi xi 1 )...(xi xn ) i 0,1,, n （5） li (x) 是 n 次多项式，满足 i j 1, li ( x j ) 0 , i j, i, j 0,1,, n （6） 令 n ln ( x ) yi li ( x ) i 0 （7） 显然ln (x) 是满足（2）的 n 次多项式，由方程（4）解的唯 一性， （7）式表示ln (x) 的与（1）式相同。 （5） 、 （7）称拉 格朗日插值多项式，ln (x) 计算插值称拉格朗日多项式插值。 分段线性插值 将两个相邻的节点用直线连起来，如此形成的一条折线就是 分段线性插值函数，记作 I n (x) ，它满足I n ( x j ) y j 且I n (x) 在每个小区间[ x j , x j 1 ] ( j 0,1,, n) 上是线性函数。 I n (x) 可以表示为 n I n ( x) y j l j ( x) j 0 x x j 1 x x , x j 1 x x j ( j 0时舍去） j 1 j x x j 1 l j ( x) , x j x x j 1 ( j n时舍去） x j x j 1 0, 其它 注 ： I n (x) 有 良 好 的 收 敛 性 ， 即 对 于 x [a, b] 有 x 点的插值时，只用到 x 左 lim I n ( x ) g ( x ) 。用 I n (x) 计算 n n 无关。但 n 越大，分段越 右的两个节点，计算量与节点个数 多，插值误差越小。实际上用函数表作插值计算时，分段线 性插值就够了，如数学、物理中用的特殊函数表，数理统计 中用的概率分布表等。 三次样条插值 三次样条函数记作 S（x）,a x b ,要求它满足以下条件： a) 在每个小区间xi 1, xi ,i 1,2,, n 上是 3 次多项式； b) 在 a x b 上二阶导数连续； c) S ( xi ) yi ，i 1,2,, n (14) 由条件 a，不妨将 S（x）记为 S ( x ) si ( x ), x xi 1, xi , i 1,...n S i ( x ) ai x 3 bi x 2 ci x d i (15) 其中ai ，bi ，ci ，d i 为待定系数，共 4n 个。 由条件 b: Si ( xi ) Si 1 ( xi ) Si ( xi ) Si1 ( xi ) S ( x ) S ( x ) i i i 1 i i 1,2,, n 1 (16) 容易看出， （14） 、 （16）式共含有 4n-2 个方程，为确定 S（x） 的 4n 个待定参数，尚需再给出 2 个条件。最常用的是所谓 自然边界条件： S ( x0 ) S ( x0 ) （17） 可以证明，4n 阶线性方程组（14） 、 （16）、 （17）有唯一解， 即 S（x）被唯一确定。但是，这种解法的工作量太大，方 程组又常呈病态，所以实际上要设计简便的解法。 曲线拟合问题的提法: 已知一组（二维）数据，即平面上的 n 个点( xi , yi ) ， i 1,2,, n, xi 互不相同，寻求一个函数（曲线）y f (x) ， 使 f (x) 在某中准则下与所有数据点最为接近，即曲线拟合 得最好，如图: y +( xi , yi ) i + + + + + + O y f (x) + + x 线性最小二乘法是解决曲线拟合最常用的方法， 基本思路是令： f ( x) a1r1 ( x) a2r2 ( x) amrm ( x) 其 中 rk (x) 是 事 先 选 定 的 一 组 函 数 ， ak 是 待 定 系 数 (k 1,2,, m, m n). 拟 合 准 则 是 使 n 个 点( xi , yi ) ， i 1,2,, n ，与 y f ( xi ) 的距离 i 的平方和最小，称最小 二乘准则。 一、系数的确定 记 n n i 1 i 1 J (a1 ,, am ) i2 [ f ( xi ) yi ]2 为求 a1,, am 使 J 达到最小，只需利用极值的必要条件 J 0(k 1,, m) 得到关于a1,, am 的线性方程组： ak m n r1 ( xi )[ ak rk ( xi ) yi ] 0 k 1 i 1 m n rm ( xi )[ ak rk ( xi ) yi ] 0 k 1 i 1 二、常用的曲线函数：rk ( x) 1、直线： y a1x a2 m y a x am x am1 2、多项式： 1 a1 3、双曲线（一支） ： y a2 x a x 4、指数曲线： y a1e 2 船在该海域会搁浅吗？ 在某海域测得一些点( x, y ) 处的水深 z（单位：英尺） 由下表给出，水深数据是在低潮时测得的。船的吃水 深度为 5 英尺，问在矩形(75,200) (50,150) 里的 哪些地方船要避免进入。 水道水深测量数据（单位：英尺） x 129.0 140.0 103.5 88.0 185.5 195.0 105.5 Y 7.5 Z 4 X 141.5 23.0 147.0 22.5 137.5 85.5 8 6 157.5 107.5 77.0 8 6 8 8 81.0 162.0 162.0 117.5 Y -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5 Z 9 9 8 8 9 4 9 一、问题分析： 假设：该海域海底是平滑的。由于测量点是散乱分布的，先 在平面上作出测量点的分布图，在利用二维插值方法补充一 些点的水深，然后作出海底曲面图和等高线图，并求出水深 小于5的海域范围。 二、问题求解： 1、作出测量点 的分布图： 150 100 50 0 -50 -100 60 80 100 120 140 160 180 200 2、作出海底地貌图 3、危险区域海底地貌图 4、危险区域平面图 150 100 50 0 -50 80 100 120 140 160 180 200 薄膜渗透率的测定 某种医用薄膜有允许一种物质的分子穿 透它，从高浓度的溶液向低浓度的溶液 扩散的功能，在试制时需测定薄膜被这VB 种分子穿透的能力。测定方法如下： 用面积 S 的薄膜将容器分成体积分别为V A ,VB 的两部分，在两 部分中分别注满该物质的两种不同浓度的溶液。此时该物质分 子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。通过单位 面积膜分子扩散的速度与膜两侧溶液的浓度差成正比，比例系 数 K 表证了薄膜被该物质分子穿透的能力，称为渗透率。定时 测量容器中薄膜某一侧的溶液浓度值，以此确定 K 的值。 VA S 一、假设 1、薄膜两侧的溶液始终是均匀的，即在任何时刻膜两侧的每 一处溶液的浓度都是相等的 2、当两溶液的浓度不一致时，物质的分子穿透薄膜总是从高 浓度溶液向低浓度溶液扩散 3、通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成 正比 4、薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的 性能是相同的 二、符号说明 1、C A (t ), CB (t ) 表示 t 时刻膜两侧的溶液浓度； 2、 A , B 表示初始时刻膜两侧的溶液浓度（单位： 毫克/立方厘米） ； 3、K 表示渗透率； 4、VA ,VB 表示由薄膜阻隔的容器两侧的体积； 三、建模 考察时段[t,t+Δt]薄膜两侧容器中该物质质量的变化。以容 器A为例，在该时段物质质量的增加量为： VAC A (t t ) VAC A (t ) 另一方面从B侧渗透至A侧的该物质质量为： SK(CB C A )t 由质量守恒定律有： VACA (t t ) VACA (t ) SK(CB CA )t 由此得： dCA SK (CB C A ) dt VA 又整个容器 中含有该物质的质量应该不变，所以有下式： VACA (t ) VBCB (t ) VA A VB B 即 所以 VB VB C A (t ) A B CB (t ) VA VA dCB (t ) 1 1 A B SK( )CB SK( ) dt VA VB VB VA 在利用初始条件 CB (t ) CB (0) B 得 AVA BVB VA VB VA ( B A ) e VA VB SK ( 1 1 )t V A VB 至 此 ， 问 题 归 结 为 利 用 CB 在 时 刻 t j 的 测 量 数 据 C j ( j 1,2,, N ) 来辩识参数 K 和 A , B ，对应的数学 模型为求函数： N SK ( E ( K , a, b) [a be j 1 1 1 )t j V A VB Cj] 的最小值点(K,a, b)，其中： a AVA BVB VA VB VA ( B A ) ,b VA VB 2 四、模型求解 设VA VB 1000立方厘米，S=10 平方厘米，求容器的 B 部分溶液浓度的测试结果如下表（其中C j 的单位为 毫克/立方厘米） t j (秒) 100 200 300 400 500 5 4.54 c c jj(10 ) 4.99 5.35 5.65 5.90 t j (秒) 600 700 800 900 1000 c j (105 ) 6.10 6.26 6.39 6.50 6.59 此时极小化的函数为： 10 E ( K , a, b) [a be j 1 用Matlab软件进行计算 20 K t j Cj ] 2