当前位置:首页 >> 政史地 >>

拟合与插值专题final(余绍权) (1)


2012数学建模培训 ----拟合与插值专题
余绍权 2012/10/28

? 在大量的应用领域中,人们经常面临用一个解析 函数描述数据(通常是测量值)的任务。对这个问 题有两种方法。
? 一种是插值法,数据假定是正确的,要求以某种方法描述数 据点之间所发生的情况。 ? 另一种方法是曲线拟合或回归。人们设法找出某条光滑曲线, 它最佳地拟

合数据,但不必要经过任何数据点。

? 本专题的主要目的是:了解插值和拟合的基本内容; 掌握用Matlab求解插值与拟合问题的基本命令。

内容提纲

?1.拟合问题引例及基本理论 ?2.Matlab求解拟合问题 ?3.应用实例 ?4.插值问题引例及基本理论 ?5.Maltab求解插值问题 ?6.应用实例

拟合问题

数据拟合在很多赛题中有应用,与图形
处理有关的问题很多与插值和拟合有关系,

例如98年美国赛A题,生物组织切片的三维插
值处理,94年A题逢山开路,山体海拔高度的

插值计算,2003年的“非典”问题也要用到
数据拟合算法,观察数据的走向进行处理,

2005年的雨量预报的评价的插值计算。2001
年的公交车调度拟合问题,2003年的饮酒驾

车拟合问题。

拟 合 问 题 引 例 1
温度t(0C) 20.5 32.7 51.0 73.0 95.7 已知热敏电阻数据: 电阻R(?) 765 求600C时的电阻R。
1100 1000 900 800 700 20

826

873

942 1032

设 R=at+b a,b为待定系数

40

60

80

100

拟 合 问 题 引 例 2 已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg) t (h) 0.25 0.5 1 1.5 2 3 4 6 8

c (?g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01
求血药浓度随时间的变化规律c(t). 作半对数坐标系(semilogy)下的图形
10
2

MATLAB(aa1)
? kt

c (t ) ? c0 e
10
1

c0 , k为待定系数
0 2 4 6 8

10

0

曲 线 拟 合 问 题 的 提 法
已知一组(二维)数据,即平面上 n个点(xi,yi) i=1,…n, 寻求一个函数(曲线)y=f(x), 使 f(x) 在某种准则下与所 有数据点最为接近,即曲线拟合得最好。 y + + + +

+ ?i (x+ i) i,y

+

+

+

y=f(x)

x

?i 为点(xi,yi) 与曲线 y=f(x) 的距离

线性最小二乘拟合 f(x)=a1r1(x)+ …+amrm(x)中 函数{r1(x), …rm(x)}的选取 1. 通过机理分析建立数学模型来确定 f(x);
2. 将数据 (xi,yi) i=1, …n 作图,通过直观判断确定 f(x):

f=a1+a2x + + +
+ + f=a1+a2/x + + +

f=a1+a2x+a3x2

f=a1+a2x+a3x2 + + + + f=ae-bx + + +

+
+ + + +

+

f=aebx +

+

+ +

+ + +

+ +

曲线拟合问题最常用的解法——线性最小二乘法的基本思路
第一步:先选定一组函数 r1(x), r2(x), …rm(x), m<n, 令 f(x)=a1r1(x)+a2r2(x)+ …+amrm(x) 其中 a1,a2, …am 为待定系数。 第二步: 确定a1,a2, …am 的准则(最小二乘准则): (1)

使n个点(xi,yi) 与曲线 y=f(x) 的距离?i 的平方和最小 。


J ( a1 , a2 , ? am ) ? ? ? i2 ? ? [ f ( xi ) ? yi ]2
i ?1 n i ?1

n

n

? ? [? ak rk ( xi ) ? yi ]
i ?1 k ?1

m

2

( 2)

问题归结为,求 a1,a2, …am 使 J(a1,a2, …am) 最小。

线性最小二乘法的求解:预备知识 超定方程组:方程个数大于未知量个数的方程组

? r11a1 ? r12 a2 ? ? ? r1m am ? y1 ? ?? ( n ? m) ? ?r a ? r a ? ? ? r a ? y nm m n ? n1 1 n 2 2
? r11 R ? ?? 其中 ? ? rn1 ? r12 ? rn 2
n

即 Ra=y

? ? ?

r1m ? ? a1 ? ? y1 ? ? ?, a ? ? ? ?, y ? ? ? ? ? ? ? ? ? ? am ? ? yn ? rnm ? ? ? ? ? ?

(ri1a1 ? ri 2 a2 ? ? ? rim am ? yi ) 2 达到最小, 如果有向量a使得 ?
i ?1

则称a为上述超定方程的最小二乘解。

线性最小二乘法的求解

所以,曲线拟合的最小二乘法要解决的问题,实际上就是求以 下超定方程组的最小二乘解的问题。 Ra=y (3) ? r1 ( x1 ) ? rm ( x1 ) ? ? a1 ? ? y1 ? ?, a ? ? ? ?, y ? ? ? ? 其中 R ? ??? ? ? ? ? ? ? ? r1 ( xn ) ? rm ( xn ) ? ? am ? ? yn ? ? ? ? ? ? ?
定理:当RTR可逆时,超定方程组(3)存在最小二乘解, 且即为方程组 RTRa=RTy 的解:a=(RTR)-1RTy ------正则(正规)方程组

用MATLAB解拟合问题

1、线性最小二乘拟合

2、非线性最小二乘拟合

用MATLAB作线性最小二乘拟合 1. 作多项式f(x)=a1xm+ …+amx+am+1拟合,可利用已有命令: a=polyfit(x,y,m) 输出拟合多项式系数 a=[a1,…,am,am+1]’ (数组) 输入同长度 数组X,Y

拟合多项式 次数

2.多项式在x处的值y的计算命令:y=polyval(a,x) 3.对超定方程组

Rn?m am?1 ? yn?1 (m ? n) ,用 a ? R \ y

可得最小二乘意义下的解。

例 对下面一组数据作二次多项式拟合
xi yi 0 0.1 0.2 3.28 0.3 6.16 0.4 7.08 0.5 7.34 0.6 7.66 0.7 9.58 0.8 9.48 0.9 1

-0.447 1.978

9.30 11.2

2 即要求 出二次多项式: f ( x) ? a1x ? a2 x ? a3

中 的 A ? (a1 , a2 , a3 )? 使得:

[ f ( xi ) ? yi ]2 最小 ?
i ?1

11

? x , x1 , 1 ? ? y1 ? ? 2 ? ? a1 ? ? ? ? x2 , x2 , 1 ? ? a ? ? ? y2 ? ???? ? ? 2 ? ?? ? ? ? ? a3 ? ? ? ? ? 2 ? x11 , x11 , 1? ? y11 ? ? ?
2 1

解法1.解超定方程的方法
此时 ? x12 ? R ? ?? ? x2 ? 11 1? ? ? ?? x11 1 ? ? x1
RA ? y?

1)输入命令: x=0:0.1:1;

y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];

R=[(x.^2)', x', ones(11,1)];
A=R\y'
MATLAB(zxec1)

2)计算结果: A = [-9.8108, 20.1293, -0.0317]

f ( x) ? ?9.8108x 2 ? 20.1293x ? 0.0317

12

解法2.用多项式拟合的命令
1)输入命令: x=0:0.1:1;

10 8 6 4 2 0 -2

0

0.2

0.4

0.6

0.8

1

y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2]; A=polyfit(x,y,2)

z=polyval(A,x);
plot(x,y,'k+',x,z,'r') %作出数据点和拟合曲线的图形
MATLAB(zxec2)

2)计算结果:A = [-9.8108, 20.1293, -0.0317]

f ( x) ? ?9.8108x 2 ? 20.1293x ? 0.0317

用MATLAB作非线性最小二乘拟合
两个求非线性最小二乘拟合的函数: lsqcurvefit、lsqnonlin。 相同点和不同点:两个命令都要先建立M-文件fun.m,定义函 数f(x),但定义f(x)的方式不同,请参考例题。 1. lsqcurvefit 已知数据点:xdata=(xdata1,xdata2,…,xdatan) ydata=(ydata1,ydata2,…,ydatan) lsqcurvefit用以求向量值函数 F(x,xdata)=(F(x,xdata1),…,F(x,xdatan))T

中的参变量x(向量),使得

1 n ? ( F ( x, xdatai ) ? ydatai ) 2 i ?1

2

最小

输入格式: (1) x = lsqcurvefit (‘fun’,x0,xdata,ydata); (2) x =lsqcurvefit (‘fun’,x0,xdata,ydata,lb, ub); (3) x =lsqcurvefit (‘fun’,x0,xdata,ydata, lb, ub, options); (4) [x, resnorm] = lsqcurvefit (‘fun’,x0,xdata,ydata,…); (5) [x, resnorm, residual] = lsqcurvefit (‘fun’,x0,xdata,ydata,…);

说明:x = lsqcurvefit (‘fun’,x0,xdata,ydata,options); fun是一个事先建立的 定义函数F(x,xdata) 的 M-文件, 自变量为x和 xdata 属性 迭代初值 已知数据点

2. lsqnonlin
已知数据点: xdata=(xdata1,xdata2,…,xdatan) ydata=(ydata1,ydata2,…,ydatan) lsqnonlin用以求向量值函数 f(x)=(f1(x),f2(x),…,fn(x))T
T 2

中的参量x,使得
2 2

f ( x) f ( x) ? f1 ( x) ? f 2 ( x) ? ? ? f n ( x)
最小。 其中 fi(x)=f(x,xdatai,ydatai) =F(x,xdatai)-ydatai

输入格式:
1) x=lsqnonlin(‘fun’,x0); 2) x= lsqnonlin (‘fun’,x0,lb,ub); 3) x= lsqnonlin (‘fun’,x0, ,lb,ub,options); 4) [x, resnorm]= lsqnonlin (‘fun’,x0,…); 5) [x, resnorm , residual]= lsqnonlin (‘fun’,x0,…); 说明:x= lsqnonlin (‘fun’,x0,options); fun是一个事先建立的 定义函数f(x)的M-文件, 自变量为x

属性
迭代初值

例2 用下面一组数据拟合 c(t ) ? a ? be0.0.2 kt
中的参数a,b,k
tj
100 200 300 400 500 600 700 800 900 1000

c j ? 10?3 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59

该问题即解最优化问题:
1 1 10 ?0.02 kt j min F (a, b, k ) ? ? [a ? be ? c j ]2 2 2 j ?1

解法1. 用命令lsqcurvefit

F(x,tdata)= (a ? be?0.02 kt 1 ,?, a ? be ?0.02 kt 10 )T ,x=(a,b,k) 1)编写M-文件 curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中 x(1)=a; x(2)=b;x(3)=k; 2)输入命令 tdata=100:100:1000 cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39, 6.50,6.59]; x0=[0.2,0.05,0.05]; x=lsqcurvefit ('curvefun1',x0,tdata,cdata) f= curvefun1(x,tdata)
MATLAB(fzxec1)

3)运算结果:
f =0.0043 0.0062 x =0.0063 0.0051 0.0062 -0.0034 0.0056 0.0063 0.2542 0.0059 0.0063 0.0061 0.0063

4)结论:a=0.0063, b=-0.0034, k=0.2542

c ? t ? ? 0.0063 ? 0.0034e

?0.02?0.2542 t

? 0.0063 ? 0.0034e ?0.005084t

解法 2 用命令lsqnonlin

f ? x ? ? F ? x, tdata , ctada ? ? ( a ? be ?0.02 kt1 ? c1 ,? , a ? be ?0.02 kt10 ? c10 )T x=(a,b,k) 函数curvefun2的自变量是x,
cdata和tdata是已知参数,故应 1)编写M-文件 curvefun2.m 将cdata tdata的值写在 function f=curvefun2(x) curvefun2.m中

tdata=100:100:1000; cdata=1e-03*[4.54,4.99,5.35,5.65,5.90, 6.10,6.26,6.39,6.50,6.59]; f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata

2)输入命令: x0=[0.2,0.05,0.05]; x=lsqnonlin('curvefun2',x0) f= curvefun2(x)

MATLAB(fzxec2)

3)运算结果为 f =1.0e-003 *(0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792) x =0.0063 -0.0034 0.2542 4)结论:即拟合得a=0.0063 b=-0.0034 k=0.2542

可以看出,两个命令的计算结果是相同的。

说明:拟合与统计回归
? 区别与联系
? 统计回归
– 线性回归 (regress命令) – 非线性回归

非线性回 归 1、回归:

是事先用m-文件定 义的非线性函数

(1)确定回归系数的命令: [beta,r,J]=nlinfit(x,y,’model’, beta0) 估计出的 回归系数 残差 输入数据x、y分别为 n? m矩阵和n维列向 量,对一元非线性回 归,x为n维列向量。 回归系数 的初值

Jacobian矩阵

(2)非线性回归命令:nlintool(x,y,’model’, beta0,alpha) 2、预测和预测误差估计: [Y,DELTA]=nlpredci(’model’, x,beta,r,J) 求nlinfit 或nlintool所得的回归函数在x处的预测值Y及预测值的显 著性为1-alpha的置信区间Y ? DELTA.

例题的求解:
1、对将要拟合的非线性模型 y=a e b / x ,建立 m-文件 volum.m 如下: function yhat=volum(beta,x) yhat=beta(1)*exp(beta(2)./x);

2、输入数据: x=2:16; y=[6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76]; beta0=[8 2]'; To MATLAB(liti41) 3、求回归系数: [beta,r ,J]=nlinfit(x',y','volum',beta0); beta 即得回归模型为: 得结果:beta = 1.10641 11.6036 ? y ? 11.6036 e x -1.0641

4、预测及作图: [YY,delta]=nlpredci('volum',x',beta,r ,J); plot(x,y,'k+',x,YY,'r')

插值问题

在实际中,常常要处理由实验或测量所 得到的一些离散数据。插值与拟合方法就是 要通过这些数据去确定某一类已知函数的参 数或寻求某个近似函数,使所得到的近似函 数与已知数据有较高的拟合精度。 如果要求这个近似函数(曲线或曲面) 经过所已知的所有数据点,则称此类问题为 插值问题。 (不需要函数表达式)

如果不要求近似函数通过所有数据点, 而是要求它能较好地反映数据变化规律的近 似函数的方法称为数据拟合。(必须有函数 表达式) 近似函数不一定(曲线或曲面)通过所

有的数据点。

1、联系

都是根据实际中一组已知数据来构造一个能够
反映数据变化规律的近似函数的方法。 2、区别 插值问题不一定得到近似函数的表达形式,仅 通过插值方法找到未知点对应的值。数据拟合 要求得到一个具体的近似函数的表达式。

插值的使用及求解
当数据量不够,需要补充,且认定已有数 据可信时, 通常利用函数插值方法。 实际问题当中碰到的函数 f (x) 是各种各 样的,有的表达式很复杂,有的甚至给不出数 学的式子,只提供了一些离散数据,警如,某 些点上的函数值和导数值。

拟合与插值
问题:给定一批数据点,需确定满足特定要求的曲线或曲面 解决方案: ?若要求所求曲线(面)通过所给所有数据点,就是插值问题;

?若不要求曲线(面)通过所有数据点,而是要求它反映对象 整体的变化趋势,就是数据拟合,又称曲线拟合或曲面拟合。 说明:
函数插值与曲线拟合都是要根据一组数据构造一个函数作为近 似,由于近似的要求不同,二者在数学方法上完全不同。 实例:下面数据是某次实验所得,希望得到x和 f之间的关系?
x f 1 1.5 2 3.9 4 6.6 7 11.7 9 15.6 12 13 18.8 19.6 15 20.6 17 21.1

MATLAB(cn)

最临近插值、线性插值、样条插值与曲线拟合结果:
25

25

??????????
20

20

15

?? ? ? ??? ?? à????

nearest

15

linest ?????à????????

???? ???? ?? 10

10

5
5

0
0 0 2 4 6 8 10 12 14 16 18

0

2

4

6

8

10

12

14

16

18

25

?????????? 20

15 spline

10 ?????à???????? 5

0

0

2

4

6

8

10

12

14

16

18

一 一、插值的定义 二、插值的方法







拉格朗日插值

分段线性插值

三次样条插值 三、用Matlab解插值问题

二维插值
一、二维插值定义 二、网格节点插值法

最邻近插值 分片线性插值 双线性插值

三、用Matlab解插值问题

网格节点数据的插值
散点数据的插值
返回

一维插值的定义
已知 n+1个节点 ( x j , y j ) ( j ? 0,1, ? n, 其中

xj 互不相同,不妨设 a ? x0 ? x1 ? ? ? xn ? b),

求任一插值点

x (? x j ) 处的插值 y * .
*

节点可视为由

y
y0

*

y ? g (x)
?
?

产生,

y1

?

?

?

g 表达式复杂,
或无封闭形式, 或未知。

?

x0 x1 x

*

xn

构造一个(相对简单的)函数 y ? f (x), 通过全部节点, 即

f ( x j ) ? y j ( j ? 0,1, ? n)
再用

f (x) 计算插值,即 y ? f ( x ).
* *

y
y0

*

y1

? ?

?

?

?

?

x0 x1 x *

xn
返回

拉格朗日(Lagrange)插值
已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为 y0,y1,…,yn 。求一个n次多项式函数Pn(x),使其满足: Pn(xi)=yi,i=0,1,…,n. 解决此问题的拉格朗日插值多项式公式如下

Pn ( x ) ? ? L i ( x ) ? y i
i?0

n

其中Li(x) 为n次多项式:
( x ? x 0 )( x ? x 1 ) ? ( x ? x i ?1 )( x ? x i ?1 ) ? ( x ? x n ) L i (x) ? ( x i ? x 0 )( x i ? x 1 ) ? ( x i ? x i ?1 )( x i ? x i ?1 ) ? ( x i ? x n )

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

拉格朗日(Lagrange)插值
特别地: 两点一次(线性)插值多项式:

x ? x0 x ? x1 L1 ?x ? ? y0 ? y1 x0 ? x1 x1 ? x0
三点二次(抛物)插值多项式:

?x ? x1 ?? ?x ? x2 ? y ? ?x ? x0 ?? ?x ? x2 ? y ? ?x ? x0 ?? ?x ? x1 ? y L2 ?x ? ? 0 1 2 ?x0 ? x1 ?? ?x0 ? x2 ? ?x1 ? x0 ?? ?x1 ? x2 ? ?x2 ? x0 ?? ?x2 ? x1 ?
直接验证可知, Ln ?x ?满足插值条件.



1 g ( x) ? , ?5? x ? 5 2 1? x

采用拉格朗日多项式插值:选取不同插值 节点个数n+1,其中n为插值多项式的次数,当n 分别取2,4,6,8,10时,绘出插值结果图形.
To Matlab lch(larg1)

拉格朗日多项式插值的 这种振荡现象叫 Runge(龙格)
现象
返回

高次插值的龙格现象
对于代数插值来说,插值多项式的次数很 高时,逼近效果往往很不理想。例如,考察 5 p 5 函数 f ? x? ? 1/ ?1 ? x ? , ?5 ? x ? ,设将区间 ?-5,?分为 n 等份,n ? x ? 表取n ? 1 个等分点作节点的插值多项式,如下 pn 图所示,当 n 增大时,? x ? 在两端会发出激烈 的振荡,这就是所谓龙格现象。
2

龙格现象
2 y=1/(1+x 2) y=p4(x) y=p10(x) 1.5

1

0.5

0

-0.5 -5

-4

-3

-2

-1

0 x

1

2

3

4

5

分段线性插值
y ? ? ? o
Ln ( x ) ?

?

?

? xj-1 xj xj+1 xn x

x0

?y
j ?0

n

j

l j ( x)

? x ? x j ?1 , x j ?1 ? x ? x j ? n越大,误差越小. ? x j ? x j ?1 ? x ? x j ?1 ? l j ( x) ? ? , x j ? x ? x j ?1 lim Ln ( x) ? g ( x), x0 ? x ? x j ? x j ?1 n ?? ?0, 其它 ? ? ?

计算量与n无关;

? xn

样条函数与样条插值
样条函数的概念
定义:设给定区间 [a, b] 的一个分划 ? : a ? x0 ? x1 ? ? ? ? ? xn ? b , 如果函数 s (x ) 满足条件:
1)在每个子区间 [ xi ?1 , xi ](i ? 1,2,? ? ?, n) 上是

k 次多项式;

2) s (x ) 及直到 k -1 阶的导数在 [a, b] 上连续。

则 称 s (x ) 是 关 于 分 划 △ 的 一 个 k 次 多 项 式 样 条 函 数 ,

x0 , x1 ,? ? ?, xn 称为样条节点, x1 , x2 ,? ? ?, xn?1 称为内节点, x0 , x n
称为边界节点,全体记作 S P (?, k ) ,称为 k 次样条函数空间。
48 2013年8月13日

三次样条插值
比分段线性插值更光滑。
y
? ? ? ? ? ?

? ?

?

a

xi-1

xi

b

x

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

三次样条插值

S ( x) ? {si ( x), x ? [ xi ?1 , xi ], i ? 1,? n}
1) si ( x ) ? ai x 3 ? bi x 2 ? ci x ? d i (i ? 1, ? n) 2) S ( xi ) ? yi (i ? 0,1, ? n) 3) S ( x ) ? C 2 [ x0 , xn ]

? si ( xi ) ? si ?1 ( xi ), si?( xi ) ? si??1 ( xi ), si??( xi ) ? si???1 ( xi ) (i ? 1, ?, n ? 1)
4) S ??( x0 ) ? S ??( xn ) ? 0 ( 自然边界条件) 2) 3) 4) ? ai , bi , ci , d i ? S ( x)
lim S ( x) ? g ( x )
n ??

g(x)为被插值函数。

用MATLAB作插值计算
一维插值函数:

yi=interp1(x,y,xi,'method')
xi处的插 值结果 插值节点 被插值点 插值方法

‘nearest’ :最邻近插值 ‘linear’ : 线性插值; ‘spline’ : 三次样条插值; ‘cubic’ : 立方插值。 缺省时: 分段线性插值。

注意:所有的插值方法都要求x是单调的,并且xi不 能够超过x的范围。

例:在1-12的11小时内,每隔1小时测量一次温 度,测得的温度依次为:5,8,9,15,25,29,31, 30,22,25,27,24。试估计每隔1/10小时的温度 值。
hours=1:12; temps=[5 8 9 15 25 29 31 30 22 25 27 24]; h=1:0.1:12; t=interp1(hours,temps,h,'spline'); plot(hours,temps,'+',h,t,hours,temps,'r:') %作图 xlabel('Hour'),ylabel('Degrees Celsius’) To MATLAB (temp)

估计结果:
35 original spline linear

30

25

Degrees Celsius

20

15

10

5

0

2

4

6 Hour

8

10

12



已知飞机下轮廓线上数据如下,求x每改变0.1时的y值。

X Y

0 0

3 1.2

5 1.7

7 2.0

9 2.1
y

11 2.0

12 1.8

13 1.2

14 1.0

15 1.6

机翼下 轮廓线

? ?

?

?

?

?
? ?

?

x

To MATLAB(plane) 返回

lagrange 20

机翼下 轮廓线

0

-20 0 5 piecewise linear 4 10 15

2

0 0 5 spline 4 10 15

2

0 0 5 10 15

二维插值的定义
第一种(网格节点):
y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

O

x

已知 m?n个节点
其中 互不相同,不妨设

构造一个二元函数

通过全部已知节点,即

再用

计算插值,即

第二种(散乱节点):
y ? ? ? ? ? ? ? ?

?
? ? ? ?

? ?

0

x

已知n个节点
其中 互不相同,

构造一个二元函数

通过全部已知节点,即

再用

计算插值,即
返回

最邻近插值
y ?

(x1, y2) ?

(x2, y2) ?

? ? ?

?

?
?

? ? (x1, y1) (x2, y1)
? ?

?
?

O

x

二维或高维情形的最邻近插值,与被插值点最邻近的 节点的函数值即为所求。 注意:最邻近插值一般不连续。具有连续性的最简单 的插值是分片线性插值。
返回

分片线性插值
y
? ? (xi, yj+1) (xi+1, yj+1) ? ? ? ? ? (xi, yj) (xi+1, yj) ?

?
?

?

?

?

?

?
x

O

将四个插值点(矩形的四个顶点)处的函数值依次 简记为: f (xi, yj)=f1,f (xi+1, yj)=f2,f (xi+1, yj+1)=f3,f (xi, yj+1)=f4

分两片的函数表达式如下:

第一片(下三角形区域): (x, y)满足

插值函数为: (y ? yj) ( x ? xi ) f ( x , y ) ? f1 ? ( f 2 ? f 1 ) ? ( f3 ? f 2 ) ( x i ?1 ? x i ) ( y i ?1 ? y j ) 第二片(上三角形区域):(x, y)满足 y j?1 ? y j y? (x ? x i ) ? y i x i ?1 ? x i 插值函数为: (y ? yj) ( x ? xi ) f ( x , y ) ? f1 ? ( f 4 ? f 1 ) ? ( f3 ? f 4) ( y i ?1 ? y j ) ( x i ?1 ? x i ) 注意:(x, y)当然应该是在插值节点所形成的矩形区 返回 域内。显然,分片线性插值函数是连续的;

双线性插值
y
? ? ?

(x1, y2) ?

(x2, y2) ?

? ? ?

? ? ? x

? ? (x1, y1) (x2, y1) ? ?

O

双线性插值是一片一片的空间二次曲面构成。 双线性插值函数的形式如下: f ( x , y ) ? ( ax ? b )( cy ? d ) 其中有四个待定系数,利用该函数在矩形的四个顶 点(插值节点)的函数值,得到四个代数方程,正 好确定四个系数。 返回

用MATLAB作网格节点数据的插值
z=interp2(x0,y0,z0,x,y,’method’)

被插值点 的函数值

插值 节点

被插值点

插值方法

‘nearest’ 最邻近插值 ‘linear’ 双线性插值 ‘cubic’ 双三次插值 缺省时, 双线性插值
要求x0,y0单调;x,y可取为矩阵,或x取 行向量,y取为列向量,x,y的值分别不能超出 x0,y0的范围。

例:测得平板表面3*5网格点处的温度分别为: 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)的图形。
1.先在三维坐标画出原始数据,画出粗糙的温度分布曲面图. 输入以下命令: x=1:5; y=1:3; temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86]; mesh(x,y,temps) 2.以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.

再输入以下命令:
xi=1:0.2:5; yi=1:0.2:3;

zi=interp2(x,y,temps,xi',yi,'cubic');
mesh(xi,yi,zi) 画出插值后的温度分布曲面图. To MATLAB (wendu)

实验结果对比:插值前(左),插值后(右)

例 山区地貌:
在某山区测得一些地点的高程如下表。平面区域为 1200<=x<=4000,1200<=y<=3600) 试作出该山区的地貌图和等高线图,并对几种插值方法进行比较。
X Y 1200 1600 2000 2400 2800 3200 3600 1200 1130 1320 1390 1500 1500 1500 1480 1600 1250 1450 1500 1200 1200 1550 1500 2000 1280 1420 1500 1100 1100 1600 1550 2400 1230 1400 1400 1350 1550 1550 1510 2800 1040 1300 900 1450 1600 1600 1430 3200 900 700 1100 1200 1550 1600 1300 3600 500 900 1060 1150 1380 1600 1200 4000 700 850 950 1010 1070 1550 980

通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插 值效果进行比较。

To MATLAB (moutain) 返回

用MATLAB作散点数据的插值计算
插值函数griddata格式为:

cz =griddata(x,y,z,cx,cy,‘method’) 被插值点 的函数值
插值 节点 被插值点 插值方法

‘nearest’ 最邻近插值 ‘linear’ 双线性插值 ‘cubic’ 双三次插值 'v4'- Matlab提供的插值方法 缺省时, 双线性插值

要求cx取行向量,cy取为列向量。

作业与练习

练习1 用给定的多项式,如y=x3-6x2+5x-3,产生一
组数据(xi,yi,i=1,2,…,n),再在yi上添加随机干扰(可用

rand产生(0,1)均匀分布随机数,或用rands产生N(0,1)分
布随机数),然后用xi和添加了随机干扰的yi作的3次多 项式拟合,与原系数比较。
分别作1、2、4、6次多项式拟合,比较结果,体会欠拟合、 过拟合现象。

练习2 用电压V=10伏的电池给电容器充电,

电容器上t时刻的电压为 v(t ) ? V ? (V ? V0 )e ? ,其中
V0是电容器的初始电压, 是充电常数。试由

?

t

下面一组t,V数据确定V0, 。 ?
分别应用非线性最小二乘拟合以及非线性回归命令求解,并 作比较,体会统计回归与拟合方法的区别。

?

t (秒) V (伏)

0.5 1 2 3 4 5 7 9 6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63

练习3

在某海域测得一些点(x,y)处的水深z由下 表给出,船的吃水深度为5英尺,估计在矩形区域 (75,200)*(-50,150)里的哪些地方船要避免 进入。
用插值方法作海底曲面图.作出水深小于5的海域范围,即z=5 的等高线.

x y z x y z

129 140 103.5 88 185.5 195 7.5 141.5 23 147 22.5 137.5 4 8 6 8 6 8 157.5 -6.5 9 107.5 -81 9 77 3 8

105 85.5 8

81 162 162 117.5 56.5 -66.5 84 -33.5 8 9 4 9

Thanks!


相关文章:
2012年秋数学建模课程表
周日上午教 205 11 月 1 日,周四晚上基础楼 101 数据的插值与拟合 初等模型 余绍权 杨迪威 10 月 21 日,周日上午教一 205 10 月 25 日, 周四晚上基础...
数模培训日历2014定稿
(教三-111) :讲解 教材第三章插值拟合,SURFER ...(建模实例,要 求完整答案) ,余绍权负责 上机 上机...(要求有建模论文的格式,每次的主要内容与一专题有...
更多相关标签:
matlab插值拟合 | 插值拟合 | 插值与拟合 | 插值和拟合的区别 | 插值与拟合matlab实现 | 插值与拟合的区别 | 插值拟合 建模案例 | 插值和拟合 |