数学建模【数据处理方法(一维、二维插值方法;数据拟合方法;插值and拟合的MATLAB实现)】
目 录
3.1 一维插值方法
数据处理方法 数据处理 任务数据化 数学建模
数值算法、非数值算法
一维插值、二维插值、三维插值...
数据处理方法---比如:对未来做预测的数据处理方法。
数据查找方法、数据拟合、聚类分析、回值分析、数据逼近...
1、数据处理的应用背景
数据正确(准确)
函数通过数据点:插值法;
不要求目标函数经过所有的数据点(或者 所有的数据点都不经过!):拟合。
2、插值的基本原理
问题: 给定一批数据点,需确定满足特定要求的函数( 曲线或曲面 )。
解决方案: 若要求所求曲线(面)通过所给所有数据点,就是插值问题。
插值思想简介
yi:函数; 列表、作图
利用 已知数据(已知点及其函数值) 估算 未知数据(未知点及其函数值) :
函数f(x)的产生办法因观测数据与要求的不同而异 --> 2种方法:曲线拟合、插值
插值(精确):认为数据准确,没有误差,函数经过点;曲线拟合:存在误差,寻找误差小的函数。
3、引例
3.1、引例演示
由 下图可以看出:3条线的光滑程度 不同,
可以看出:最临近插值 光滑程度最差;线性插值 次之;样条插值 最光滑。
3.2、引例1、函数查表问题
最邻近插值(与哪个点接近,就用哪个点的函数值):【缺陷:误差大。】
3.2、引例2、绘制地图
3.3、引例3:机床加工
光滑性:导数阶数
4、一维插值定义
各个x*不同。 某函数【y*】:经过所有插值节点--->构造函数!!!
5、一维插值原理
求解插值问题的基本思路
先找到 近似函数f(x)
f(x*):插入到x0、x1、x2...xn之间的x*所对应的函数值。
构造不同的f(x),得到不同的y*。
构造什么样的函数:根据客户(题目)要求。满足客户要求。
插值:利用已知点的函数值 求 未知点的函数值。
6、一维插值方法
6.1、拉格朗日(Lagrange)插值
求一个经过n+1个插值节点的n次多项式函数。
各个Xi的值都不同,Li(X)存在,且唯一。
观察 公式:分子、分母---X、Xi。
一维:一元函数 y=f(x)
牛顿插值:略。
多项式:幂函数、常数 线性组合。【四则运算】
常数项、1次...n次:n+1个常数 --> n+1个方程
Li(xi) = 1
Pn(x) 是 n+1个多项式的和,是 n次多项式,并经过n+1个点。
function y = lagr1(x0, y0, x)
n = length(x0); m = length(x); % n=length() 长度
for k = 1:m
z = x(k); % 对x中第k个分量 进行求解,将其值赋给z
s = 0; % 初始化s
for i = 1:n
p = 1.0;
for j = 1:n
if j~=i % j!=i时
p = p*(z-x0(j))/(x0(i)-x0(j)); % 累乘 计算p函数(第i个拉格朗日函数)
end
end
s = p*y0(i) + s; % 累加
end
y(k) = s;
end
x0:n+1个点(x0、x1...xn); y0:n+1个点对应的纵坐标(y0、y1...yn); x:需要进行插值的点(未知点,可能不止一个,x1...xn;插值:根据实际要求进行处理)
输出结果(y1...ym):在已知点x0、y0(插值节点),n个点的情况下,对x进行插值 得到的输出结果。
x0 = [0 3 5 7 9 11 12 13 14 15];
y0 = [0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x = 0:0.05:15;
y = lagr1(x0, y0, x);
plot(x, y, x0, y0, ''r*)
axis([0, 15, -17, 3])
% 逼近函数 y = 1/(1+x^2)
n = 3; % 3:2次多项式 9 15
x0 = linspace(-5, 5, n);
y0 = 1./(1+x0.^2);
x = -7:0.4:7;
y_L = lagr1(x0, y0, x);
y_R = 1./(1+x.^2);
plot(x, y_R, 'r*', x, y_L);
axis([-6, 6, -1.5 2])
高次插值多项式:剧烈震荡(误差很大)
n=2:两点确定一条直线。
n=3:经过三点,存在一个二次函数。
随着 插值节点的增多,拉格朗日函数的次数会越来越高。--> 振荡现象(剧烈震荡)--> 误差增大 --> 分段插值方法
6.2、分段线性插值
n越大,插值函数越高。
一维分段线性插值:每两个相邻节点,确定一个一次函数(两点确定一条直线);存在一条折线,连接n+1个点。
每个点的左右导数不相等,光滑度低,不能满足实际的工程需要->寻找 光滑性好、次数低的、误差低的插值函数->三次样条插值
6.3、分段三次埃尔米特插值
比如:【每个区间:三次函数】
任意相邻点之间,找一个 三次多项式,使其满足要求。
6.4、三次样条插值
三次样条插值---插值函数
1)在每一个小区间
,都是一个三次函数;
一阶导数值、二阶导数值 相等;
端点处的二阶导数值 为 0。
满足的条件。每个小区间,都有一个三次,n个小区间,n+1个点,4n个常数。
三次样条插值精度 比较高
3.2 二维插值方法
1、高维插值的应用背景
在实际的应用领域中,人们经常面临寻求一个多元解析函数描述并通过高维空间中多维数据(通常是测量值)的任务。
如:天气预报风云高压图、山川河流的等高线测量 等。对此类问题多采用高维插值方法,数据假定是正确的,要求以某种方法描述数据点之间所发生的情况,以解决实际生产、生活中遇到的问题。其中,最常用的高维插值是二维插值。
2、二维插值的基本原理
二元函数 z=f(x, y)
2.1、规则网格下 二维插值
空间点 (x1, y1, z1)、(x2, y2, z2)、 ...、(xn, yn, zn) 在xOy平面投影 得来的散点图(定义域)。
投影方向 沿x、y轴方向 等距离地 画出,得到的网格,即为规则网格(容易寻找规律)。
xi:列; yj:行; 第i列 第j行
2.2、散乱节点网格下 二维插值
不规则网格查找:仅要求各个节点互不相同。散乱节点无规律,插值难度较大。
3、插值方法
二维插值方法 包括 最邻近插值法、分片线性插值、双线性插值、样条插值等一系列类似的方法,常用于图像处理、温度场计算、气压等压线和地图等高线的绘制、计算机模拟和数字信号处理等多种方面。
具体插值方法和应用参考数值逼近,本课采用matlab自带函数进行计算。
3.1、最邻近插值法
黑色点:插值节点。
最邻近插值法:按照 被插值节点 到 各插值节点的距离,最小的插值节点 所对应的函数值,就是 要求的插值结果。
离 被插值节点最近的插值节点 所对应的函数值,即为 插值。
采用 分片常值函数,一般不连续。
3.2、分片线性插值
分片线性插值:满足连续性。(整体连续)
将 方格的对角线 连到一块,分割 方格 为 2个三角形。
每个方格都这样处理 --> 可以看作:多个三角形拼接而成。
3个点 分别对应 空间中的3个点的投影。【这3个点不同;空间中的三个点也不同。】
空间中,不共线的3个点 可以唯一确定一个平面 --> 空间曲面可以看成多个三角形空间平面拼接而成。
每个平面对应一个二元线性函数 --> 所以 称之为分片线性插值。
3.3、双线性插值
4个边界,都是 直线!
3.3 数据拟合方法
数据插值:插值函数 经过 插值节点。插值节点过多 --> 数据插值误差可能较大 --> 数据拟合方法
1、拟合的应用背景
在大量的应用领域中,人们经常面临用一个解析函数描述数据(通常是测量值)的任务。对这个问题,人们设法找出某条光滑函数曲线,使它与实际数据误差尽可能地小,找出最佳拟合实际数据的函数,结果不要求目标函数经过任何数据点。这类方法统称为曲线拟合。
2、拟合的问题描述
问题:已知一组(二维)数据,即平面上 n个点(xi, yi) i=1,…,n, 寻求一个函数(曲线)y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近,即曲线拟合得最好。
解决方案:若不要求曲线(面)通过所有数据点,而是要求它反映 对象整体的变化趋势,这就是数据拟合,又称曲线拟合或曲面拟合(三维)
数据量多 --> 数据拟合 --> 反映对象整体(数据)变化趋势
3、拟合模型分类
线性(直线或平面)拟合 用一次函数或线性函数拟合已知现有数据。
非线性(曲线或曲面)拟合 若线性拟合效果不佳,为提高精度,可采用曲线(曲面)拟合数据 。常用的非线性拟合函数是二次函数,三次函数等高次多项式函数, 有时也会用到对数、指数、三角函数等。
分段拟合 在不同的数据段上,采用不同的低次多项式进行拟合。
常值---拟合-->最为简单。
线性拟合:二维:直线;三维:平面。
分段拟合:误差较大的情况下,类似于 分段插值的方法,将数据进行分段,在不同的数据段上 采用 不同的低次多项式进行拟合。
4、拟合模型常用解法
4.1、最小二乘法
最小二乘法是一个非常有用的数学方法,它的直接应用是求矛盾线性方程组的最小二乘解。在工程中,它可用来求经验公式,对实验数据进行曲线拟合;在统计学中,它可用来求最小二乘估计,多元回归等;另外,在数值分析领域,它可用来进行误差分析。 可以说,随着我们将来知识的增多,最小二乘法的应用领域将愈来愈大。
求解 拟合多项式的参数、系数
5、最小二乘法
a1,a2,a3,...,am :线性组合
f(xi) - yi :点到曲线的距离。
求得结果 越小越好 --> 拟合结果越好。
6、线性最小二乘法:预备知识
超出定解条件的方程组。
7、线性最小二乘法求解
转置矩阵R^T
7.1、化为二次型的形式
T:转置 列向量
7.2、正规方程的矩阵形式
7.3、算法说明
7.4、关于基函数的选择
3.4 数据插值的MATLAB实现
1、MATLAB作一维插值函数
x、y:数据点所构成的向量,插值节点对应的数值;
xi:被插值节点、插值向量
method:采用什么样的方法实现一维插值功能。
x,y,z:已知的插值节点数据。zi:插值结果。
method:采用什么样的方法实现二维插值功能。
常用的4种插值方法。
2、MATLAB作二维插值函数
注意:二维插值中已知数据点集(x0,y0)必须是栅格格式,一般用meshgrid函数产生。interp2要求(x0,y0)必须是严格单调的并且是等间距的,如果(x0,y0)不是等间距的,会将且变换为等间距形式,如果已知是等间距的,可在method参数前加星号,例如:'*cubic' 。
3、MATLAB作高维插值函数
例 气旋变化情况的可视化
axis([0 12 0 90 0 50]):坐标区间。
Xi,Yi,Zi:向量组。
3.5 数据拟合的MATLAB实现
1、MATLAB作线性最小二乘拟合
a:多项式系数。
x、y:拟合数据;n:拟合多项式的次数
s:误差预测结果,将数据代入后 与 实际的差值。
2、MATLAB作非线性最小二乘拟合
2.1、lsqcurvefit
lsqcurvefit、lsqnonlin :curve 刻画;non 非;lin 线性。
调用格式
2.2、lsqnonlin
3.6 数据插值的案例分析与实践
1、一维插值的应用
例1、推测24h内,每一秒的温度,绘制温度曲线
例1.1、推测13h的温度
x = 0:2:24; %时间:从0到24,增量为2
y = [12 9 9 10 18 24 28 27 25 20 18 15 13]; %温度
x1 = 13;
y1 = interp1(x, y, x1, 'spline')xi = 0:1/3600:24;
yi = interp1(x, y, xi, 'spline');
plot(x, y, 'o', xi, yi)
0、代码汇总
x = 0:2:24; %时间:从0到24,增量为2
y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 22, 24]; %温度
xi = 0:1/3600:24; %从0到24,增量为1/3600 24小时中的每一秒
yi = interp1(x, y, xi, 'nearest'); %interp1:一维多项式插值;x、y:插值数据(插值节点);nearest 最邻近插值
plot(xi, yi, 'r'); %plot():绘图命令;r:红色; 在xi时刻内所对应的插值结果记作yi
xlabel('时刻'), ylabel('温度'); %在x轴上标注‘时刻’
hold on %保留
x = 0:2:24;
y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 16, 13];
xi = 0:1/3600:24;
yi = interp1(x, y, xi, 'linear');
plot(xi, yi, 'g');
hold on
x = 0:2:24;
y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 16, 13];
xi = 0:1/3600:24;
yi = interp1(x, y, xi, 'spline');
plot(xi, yi, 'b');
hold on
x = 0:2:24;
y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 16, 13];
xi = 0:1/3600:24;
yi = interp1(x, y, xi, 'cubic');
plot(xi, yi, 'y');
hold on
1、nearest(最近插值)
2、linear(线性插值)
3、spline(三次样条插值)
4、cubic(三次插值)
例2、1-12h内,估计每隔1/10h的温度值
clc; %清除变量
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 Clesius')
引例3、机床加工
x0 = [0 3 5 7 9 11 12 13 14 15];
y0 = [0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x = 0:0.05:15;
y = lagr1(x0, y0, x);
y1 = interp1(x0, y0, x, 'spline')
subplot(2, 1, 1), plot(x, y, x0, y0, 'r*'), axis([0, 15, -17, 3])
subplot(2, 1, 2), plot(x, y1, x0, y0, 'r*')
2、二维插值的应用
例3 平板表面温度分布图
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 = 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) xi = 1:0.2:5; yi = 1:0.2:3; zi = interp2(x, y, temps, xi', yi, 'cubic'); mesh(xi, yi, zi)
例4 海域中哪些地方,船要避免进入?
1、海底曲面图(考察的海域范围)
问题分析(问题求解算法):
假设:该海域海底是平滑的。
1.输入插值基点数据;【由于测量点是散乱分布的,先作出测量点的分布图】
2.在矩形区域(75,200)x(-50,150)进行插值;【利用MATLAB中的二维插值方法求网格节点的水深】
3.作海底曲面图形和等高线图;
4.作出水深小于5的海域范围,即z=5的等高线。
%插值并作海底曲面图的程序 x = [129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5]; y = [7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -33.5]; z = [4 8 6 8 6 8 8 9 9 8 8 9 4 9]; x1 = 75:1:200; y1 = -50:1:150; [x1, y1] = meshgrid(x1, y1); %点阵 z1 = griddata(x, y, z, x1, y1, 'v4'); meshc(x1, y1, z1)
2、水深小于5的海域范围
%程序二:插值并作出水深小于5的海域范围。 x1 = 75:1:200; y1 = -50:1:150; [x1, y1] = meshgrid(x1, y1); z1 = griddata(x, y, z, x1, y1, 'v4'); %插值 z1(z1 >= 5) = nan; %将水深大于5的置为nan,这样绘图就不会显示出来 meshc(x1, y1, z1)
![]()
3.6节---代码汇总
% x = 0:2:24; %时间:从0到24,增量为2
% y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 22, 24]; %温度
% xi = 0:1/3600:24; %从0到24,增量为1/3600 24小时中的每一秒
% yi = interp1(x, y, xi, 'nearest'); %interp1:一维多项式插值;x、y:插值数据(插值节点);nearest 最邻近插值
% plot(xi, yi, 'r'); %plot():绘图命令;r:红色; 在xi时刻内所对应的插值结果记作yi
% xlabel('时刻'), ylabel('温度'); %在x轴上标注‘时刻’
% hold on %保留
% x = 0:2:24;
% y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 16, 13];
% xi = 0:1/3600:24;
% yi = interp1(x, y, xi, 'linear');
% plot(xi, yi, 'g');
% hold on
% x = 0:2:24;
% y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 16, 13];
% xi = 0:1/3600:24;
% yi = interp1(x, y, xi, 'spline');
% plot(xi, yi, 'b');
% hold on
% x = 0:2:24;
% y = [12, 9, 12, 15, 23, 30, 36, 35, 31, 27, 20, 16, 13];
% xi = 0:1/3600:24;
% yi = interp1(x, y, xi, 'cubic');
% plot(xi, yi, 'y');
% hold on
% clc; %清除变量
% 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 Clesius')
% 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)
% xi = 1:0.2:5;
% yi = 1:0.2:3;
% zi = interp2(x, y, temps, xi', yi, 'cubic');
% mesh(xi, yi, zi)
%插值并作海底曲面图的程序
% x = [129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5];
% y = [7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -33.5];
% z = [4 8 6 8 6 8 8 9 9 8 8 9 4 9];
% x1 = 75:1:200;
% y1 = -50:1:150;
% [x1, y1] = meshgrid(x1, y1); %点阵
% z1 = griddata(x, y, z, x1, y1, 'v4');
% meshc(x1, y1, z1)
%程序二:插值并作出水深小于5的海域范围。
x1 = 75:1:200;
y1 = -50:1:150;
[x1, y1] = meshgrid(x1, y1); %网格
z1 = griddata(x, y, z, x1, y1, 'v4'); %插值
z1(z1 >= 5) = nan; %将水深大于5的置为nan,这样绘图就不会显示出来
meshc(x1, y1, z1)
3.7 数据拟合的案例分析与实践
最小二乘法线性拟合的应用
例1 预测人口的增长情况
人口的增长是当前世界上引起普遍关注的问题,并且我们会发现在不同的刊物预报 同一时间的人口数字不相同,这显然是由于用了不同的人口模型计算的结果。我国是世界第一人口大国,基本上地球每九个人中就有一个中国人。有效地控制我国人口的增长是使我国全面进入小康社会、到21世纪中叶建成富强民主文明的社会主义国家的需要。而有效控制人口增长的前提是要认识人口数量的变化规律,建立人口模型,作出较准确的预报。
1949年——1994年,我国人口数据资料如下:
年份xi 1949 1954 1959 1964 1969 1974 1979 1984 1989 1994
人口数yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8[ 人口数 单位:亿!!! ]
建模分析我国人口增长的规律,预报1999年我国人口数。
1、模型(线性增长模型、指数增长模型)
2、Matlab程序and计算结果
% Matlab计算程序如下: x=[1949 1954 1959 1964 1969 1974 1979 1984 1989 1994]; y=[5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 ]; a=polyfit(x,y,1); x1=[1949:10:1994]; y1=a(2)+a(1)*x1; b=polyfit(x,log(y),1); y2=exp(b(2))*exp(b(1)*x1); plot(x,y,'*') hold on plot(x1,y1,'--r') hold on plot(x1,y2,'-k') legend('原曲线','模型一曲线','模型二曲线')
3、拟合结果与实际结果比较
4、结果比较与分析
例2 假设有一组实测数据...
假设有一组实测数据
xi 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 yi 2.3201 2.647 2.9707 3.2885 3.6008 3.909 4.2147 4.5191 4.8232 5.1275
模型假设:该数据满足的原型函数
试求出满足下面数据的最小二乘解的a,b,c,d的值。
采用最小二乘拟合:lsqcurvefit()
建立原型函数:
function y=f1(a,x)
y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);
在命令窗口中输入:
function y=f1(a,x)
y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);x=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0];
y=[2.3201 2.647 2.9707 3.2885 3.6008 3.909 4.2147 4.5191 4.8232 5.1275];
a=lsqcurvefit('f1',[1;2;2;2],x,y)
建模实例 估计水塔的水流量
1、问题的提出
2、问题分析
自动加水:自动压力器
水位升高:气压增大,停止加水。【8.2米,开始加水;10.8米,停止加水。】
高度变化 ——》水流变化、水流速度
3、模型假设
假设1:忽略因素---水位越高,流速越快(重力影响)
假设4:一直有人在用水
水流速:水流量的导数
向前差分 or 向后差分,求导数。 时刻变化 / 体积变化
利用 二阶差商公式 近似计算。。。
4、模型建立与求解
x0:时间;y0:流速;x:被插结点;y:这些点所对应的函数值(插值结果)
误差较小,结果可以接受。
5、模型的稳定性分析与检验
可能:上午、下午、晚上,用水量不同。。。所以,进行分段。
误差较小,结果可以接受。
模型评价
插值:利用 已知数据 求得 未知数据 所对应函数的方法(必须经过所有已知点)。
更多推荐
所有评论(0)