Matlab随笔之插值与拟合(下)

2024-05-11 11:18
文章标签 matlab 随笔 插值 拟合

本文主要是介绍Matlab随笔之插值与拟合(下),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

1、二维插值之插值节点为网格节点

已知m x n个节点:(xi,yj,zij)(i=1…m,j=1…n),且xi,yi递增。求(x,y)处的插值z。

         Matlab可以直接调用interp2(x0,y0,z0,x,y,`method`)

        其中 x0,y0 分别为 维和 维向量,表示节点, z0 为 × 维矩阵,表示节点值, x,y 
为一维数组,表示插值点, x 与 y 应是方向不同的向量,即一个是行向量,另一个是列 
向量, z 为矩阵,它的行数为 的维数,列数为 的维数,表示得到的插值, 'method' 
的用法同上面的一维插值。

        如果是三次样条插值,可以使用命令 
          pp=csape({x0,y0},z0,conds,valconds),   z=fnval(pp,{x,y}) 
其中 x0,y0 分别为 维和 维向量, z0 为 × 维矩阵, 为矩阵,它的行数为 的维 
数,列数为 的维数,表示得到的插值,具体使用方法同一维插值。

eg:

image

(1)、用interp2函数插值:

复制代码
x=100:100:500; 
y=100:100:400; 
z=[636 697 624 478 450 698 712 630 478 420 680 674 598 412 400 662 626 552 334 310]; 
p=100:1:500; 
q=100:1:400; 
q=q';%须为列向量 
z0=interp2(x,y,z,p,q);%分段线性插值 
z1=interp2(x,y,z,p,q,'spline');%三次线条插值 
subplot(2,1,1); 
mesh(p,q,z0); 
title('分段线性插值'); 
subplot(2,1,2); 
mesh(p,q,z1); 
title('三次线条插值'); 
%可以观察出,三次线条插值的图像更平滑 
复制代码

运行结果:

image

 

(2)、用csape函数插值:

复制代码
x=100:100:500; 
y=100:100:400; 
z=[636 697 624 478 450 698 712 630 478 420 680 674 598 412 400 662 626 552 334 310]; 
p=100:1:500; 
q=100:1:400; 
q=q'; 
%三次线条插值 
pp=csape({x,y},z');%注意跟interp2的区别,有个转置 
z0=fnval(pp,{p,q}); 
mesh(p,q,z0');%注意跟interp2的区别,有个转置 
title('三次线条插值');
复制代码

运行结果:

image

2、二维插值之插值节点为散乱节点

已知 个节点: ( xi yi zi )(= 1,2,…, n) ,求点 (xy) 处的插值 。 
对上述问题, Matlab 中提供了插值函数 griddata,其格式为: 
ZI = GRIDDATA(X,Y,Z,XI,YI) 
其中 X、 Y、 Z 均为 n 维向量,指明所给数据点的横坐标、纵坐标和竖坐标。向量 XI、 
YI 是给定的网格点的横坐标和纵坐标,返回值 ZI 为网格( XI, YI)处的函数值。 XI 
与 YI 应是方向不同的向量,即一个是行向量,另一个是列向量。

eg:

复制代码
%散乱节点的二维插值 
x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5]; 
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; 
z=-[4 8 6 8 6 8 8 9 9 8 8  9 4 9]; 
x0=[75:1:200]; 
y0=[-85:1:145]'; 
z0=griddata(x,y,z,x0,y0,'cubic');%保凹凸性3次插值 
%[xx,yy]=meshgrid(x0,y0);无需采样,故不需要该函数 
mesh(x0,y0,z0);
复制代码

运行结果:

image

在上述问题中,补上寻找最大值的程序:

%max(z0)返回一个行向量,向量的第i个元素是矩阵A的第i列上的最大值 
%find(A) 寻找矩阵A非零元素下标,返回矩阵A中非零元素所在位置 
%[i,j,v]=find(A)返回矩阵A中非零元素所在的行i,列j,和元素的值v(按所在位置先后顺序输出) 
[p,q]=find(z0==max(max(z0))); 
zmax=z0(p,q)

3、最小二乘法实现曲线拟合

(1)用最小二乘法求一个形如 bx^ 2 的经验公式:

复制代码
%等价于[1,x^2][a;b]=y,转换成解超定方程问题,超定方程的解是根据最小二乘法得来的x=[19 25 31 38 44]'; 
y=[19.0 32.3 49.0 73.3 97.8]'; 
r=[ones(5,1),x.^2] 
ab=r\y 
x0=19:0.1:44; 
y0=ab(1)+ab(2)*x0.^2; 
plot(x,y,'o',x0,y0,'r')
复制代码

(2)多项式拟合

%a=polyfit(x,y,n)用多项式求过已知点的表达式,其中x为源数据点对应的横坐标,可为行向量、矩阵,y为源数据点对应的纵坐标,可为行向量、矩阵,n为你要拟合的阶数,一阶直线拟合,二阶抛物线拟合,并非阶次越高越好,看拟合情况而定,a为m+1的行向量。polyfit函数的数学基础是最小二乘法曲线拟合原理,所得到的函数值在基点处的值与原来点的坐标偏差最小,常用于数据拟合,polyfit 做出来的值从左到右表示从高次到低次的多项式系数。

如果要求拟合函数在x`点的函数值,可以调用polyval(a,x`)函数

eg:

复制代码
x0=[1990 1991 1992 1993 1994 1995 1996]; 
y0=[70 122 144 152 174 196 202]; 
%画出散点图 
plot(x0,y0,'ro'); 
hold on 
%用线性拟合 
p=polyfit(x0,y0,1); 
z0=polyval(p,x0); 
plot(x0,z0);
复制代码

运行结果:

image

 

4、最小二乘优化 (最小二乘:least square)

image

1、lsqlin函数

image

eg:

复制代码
%拟合形如y=a+bx^2的函数 
%采样点 
x=[19 25 31 38 44]'; 
y=[19 32.3 49 73.3 97.8]'; 
r=[ones(5,1),x.^2]; 
ab=lsqlin(r,y) 
x0=19:0.1:44; 
y0=ab(1)+ab(2)*x0.^2; 
plot(x,y,'o',x0,y0,'r')
复制代码

运行结果:

image

 

5、曲线拟合与函数逼近

image

eg:

求 f(x) =cos x(-pi/2<=x<=pi/2) 在 = Span{1, x^2 , x^4} 中的最佳平方逼近多项式。

程序如下:

复制代码
syms x%定义符号数值 
base=[1,x^2,x^4]; 
y1=base.'*base 
y2=cos(x)*base.' 
r1=int(y1,-pi/2,pi/2) 
r2=int(y2,-pi/2,pi/2) 
a=r1\r2%a为符号数值 
xishu1=double(a)%化简符号数值 
digits(8)%设置符号数值的精度 
xishu2=vpa(a)%任意精度(符号类)数值 
复制代码

运行结果:

y1 =
   
[   1, x^2, x^4]
[ x^2, x^4, x^6]
[ x^4, x^6, x^8]
   
   
y2 =
   
      cos(x)
x^2*cos(x)
x^4*cos(x)
   
   
r1 =
   
[      pi,  pi^3/12,   pi^5/80]
[ pi^3/12,  pi^5/80,  pi^7/448]
[ pi^5/80, pi^7/448, pi^9/2304]
   
   
r2 =
   
                     2
            pi^2/2 - 4
pi^4/8 - 6*pi^2 + 48
   
   
a =
   
(15*(pi^4 - 308*pi^2 + 3024))/(4*pi^5)
    -(210*(pi^4 - 228*pi^2 + 2160))/pi^7
    (1260*(pi^4 - 180*pi^2 + 1680))/pi^9
  
xishu1 =
     0.9996
    -0.4964
     0.0372
   
xishu2 =
   
   0.99957952
-0.49639233
0.037209327
   
>>

  

所以y的最佳平方逼近多项式为y=0.9996-0.4964x^2+0.0372x^4

这篇关于Matlab随笔之插值与拟合(下)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/979385

相关文章

matlab读取NC文件(含group)

matlab读取NC文件(含group): NC文件数据结构: 代码: % 打开 NetCDF 文件filename = 'your_file.nc'; % 替换为你的文件名% 使用 netcdf.open 函数打开文件ncid = netcdf.open(filename, 'NC_NOWRITE');% 查看文件中的组% 假设我们想读取名为 "group1" 的组groupName

利用matlab bar函数绘制较为复杂的柱状图,并在图中进行适当标注

示例代码和结果如下:小疑问:如何自动选择合适的坐标位置对柱状图的数值大小进行标注?😂 clear; close all;x = 1:3;aa=[28.6321521955954 26.2453660695847 21.69102348512086.93747104431360 6.25442246899816 3.342835958564245.51365061796319 4.87

OpenCV结构分析与形状描述符(11)椭圆拟合函数fitEllipse()的使用

操作系统:ubuntu22.04 OpenCV版本:OpenCV4.9 IDE:Visual Studio Code 编程语言:C++11 算法描述 围绕一组2D点拟合一个椭圆。 该函数计算出一个椭圆,该椭圆在最小二乘意义上最好地拟合一组2D点。它返回一个内切椭圆的旋转矩形。使用了由[90]描述的第一个算法。开发者应该注意,由于数据点靠近包含的 Mat 元素的边界,返回的椭圆/旋转矩形数据

C# double[] 和Matlab数组MWArray[]转换

C# double[] 转换成MWArray[], 直接赋值就行             MWNumericArray[] ma = new MWNumericArray[4];             double[] dT = new double[] { 0 };             double[] dT1 = new double[] { 0,2 };

libsvm在matlab中的使用方法

原文地址:libsvm在matlab中的使用方法 作者: lwenqu_8lbsk 前段时间,gyp326曾在论坛里问libsvm如何在matlab中使用,我还奇怪,认为libsvm是C的程序,应该不能。没想到今天又有人问道,难道matlab真的能运行libsvm。我到官方网站看了下,原来,真的提供了matlab的使用接口。 接口下载在: http://www.csie.ntu.edu.

Matlab/Simulink中PMSM模型的反电动势系数和转矩系数

Matlab/Simulink中PMSM模型的反电动势系数和转矩系数_matlab pmsm-CSDN博客

MATLAB层次聚类分析法

转自:http://blog.163.com/lxg_1123@126/blog/static/74841406201022774051963/ 层次聚类是基于距离的聚类方法,MATLAB中通过pdist、linkage、dendrogram、cluster等函数来完成。层次聚类的过程可以分这么几步: (1) 确定对象(实际上就是数据集中的每个数据点)之间的相似性,实际上就是定义一个表征

MATLAB的fix(),floor()和ceil()函数的区别与联系

fix(x),floor(x)和ceil(x)函数都是对x取整,只不过取整方向不同而已。 这里的方向是以x轴作为横坐标来看的,向右就是朝着正轴方向,向左就是朝着负轴方向。 fix(x):向0取整(也可以理解为向中间取整) floor(x):向左取整 ceil(x):向右取整 举例: 4个数:a=3.3、b=3.7、c=-3.3、d=-3.7 fix(a)=3 fl

MATLAB中的eig函数

在MATLAB中,计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有5种: E=eig(A):求矩阵A的全部特征值,构成向量E。 [V,D]=eig(A):求矩阵A的全部特征值,构成对角阵D,并求A的特征向量构成V的列向量。 [V,D]=eig(A,'nobalance'):与第2种格式类似,但第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量,而格式3直接求矩阵A的特

MATLAB中的diag函数

diag函数功能:矩阵对角元素的提取和创建对角阵 设以下X为方阵,v为向量 1、X = diag(v,k)当v是一个含有n个元素的向量时,返回一个n+abs(k)阶方阵X,向量v在矩阵X中的第k个对角线上,k=0表示主对角线,k>0表示在主对角线上方,k<0表示在主对角线下方。例1: v=[1 2 3]; diag(v, 3) ans =      0     0     0