MATLAB基础运算【6】：插值拟合与数值微积分

HelloFeynman

已于 2024-08-31 11:02:37 修改

阅读量1.1k

点赞数 32

分类专栏： MATLAB操作说明系列文章标签： matlab 开发语言矩阵算法

于 2024-08-29 15:39:36 首次发布

本文链接：https://blog.csdn.net/janezhang123011/article/details/141652479

版权

MATLAB操作说明系列专栏收录该内容

9 篇文章 0 订阅

订阅专栏

系列更新MATLAB各种运算和操作的说明。

本篇主要包括插值拟合/数值微积分相关内容（可利用侧边目录直接定位）：

拉格朗日插值拟合：“自定义lagrange”
线性插值拟合：“interp1”
三次样条插值拟合：“spline”
多项式拟合曲线：“polyval”、“polyfit”
干扰下的多项式拟合
非线性插值拟合：“lsqcurvefit”
数值微积分：“gradient”、“quad”、“dblquad”、“integral2”、“int”

前文：

MATLAB基础运算【5】：图像绘制学科专题

例：对"sinx"(0-2pi),在n个节点上(n不要太大,如5~11)用拉格朗日﹑分段线性﹑三次样条三种插值方法,计算m个插值点的函数值(m要适中,如50～100).通过数值和图形输出,将三种插值结果与精确值进行比较。适当增加n,再作比较，由此作初步分析.

1- 拉格朗日插值拟合：“自定义lagrange”

取m=61，n=9，代码：

function y = lagrange(x0, y0, x)
    n = length(x0);
    m = length(x);
    for i = 1 : m % 遍历每一个需要计算插值的点
        z = x(i); % 将当前点的 x 坐标赋值给变量 z
        s = 0.0; % 初始化插值结果 s 为 0
        for k = 1 : n
            p = 1.0; % 初始化当前项的系数 p 为 1
            % 计算当前项的系数 p，排除当前已知点 x0(k)
            for j = 1 : n
                if j ~= k
                    p = p * ((z - x0(j)) / (x0(k) - x0(j))); % 根据拉格朗日插值公式计算系数 p
                end
            end
            s = s + p * y0(k); % 将当前项的值累加到插值结果 s 中
        end
        y(i) = s; % 将计算得到的插值结果赋值给 y 向量的第 i 个元素
    end
end

命令行输入：

x = [1:0.5:5];
y = [0.841471, 0.997495, 0.909297, 0.598472, 0.14112, -0.350783, -0.756802, -0.97753, -0.958924]; % 定义与 x 坐标点对应的 y 坐标值
t = 0:0.1:6; % 定义 t 坐标点，步长为 0.1，从 0 到 6
plot(t,lagrange(x,y,t),t,sin(t))

得到图像：

2- 线性插值拟合：“interp1”

取m=61，n=13，命令行输入：

x = 0:0.5:6;
y = sin(x);
t = 0:0.1:6;
yi = interp1(x,y,t);
plot(t,yi,t,sin(t))

得到图像：

增加n使得n=31：

3- 三次样条插值拟合：“spline”

取m=61，n=13，命令行输入：

x = [1:0.5:5];
y = sin(x);
xx = 0:0.1:6;
yy = spline(x,y,xx);
plot(x,y,'o',xx,yy)

得到图像：

增大n至n=16：

4- 多项式拟合曲线：“polyval”、“polyfit”

例1：x,y数据如下：

x

y

x

y

x

y

0

0.1

14

2.022

28

0.4308

2

1.884

16

2.65

30

0.203

4

2.732

18

1.5838

32

0.1652

6

3.388

20

1.35

34

-0.073

8

3.346

22

1.0082

36

-0.002

10

3

24

0.718

38

-0.1122

12

2.644

26

0.689

40

0.106

以空心圆点做出x,y的散点图（红色空心圆点），并分别采用1次、2次、3次、5次多项式拟合x，y，作出拟合曲线（蓝色实线），用subplot函数展示出2x2的四个字图，并分别用text函数将拟合结果展示在各自的子图中。

命令行输入：

x = 0:2:40;
y = [0.1, 1.884, 2.732, 3.388, 3.346, 3, 2.644, 2.022, 2.65, 1.5838, 1.35,
     1.0082, 0.718, 0.689, 0.4308, 0.203, 0.1652, -0.073, -0.002, -1.1122, 0.106]; % 定义与 x 坐标点对应的 y 坐标值
n = [1, 2, 3, 6]; % 定义多项式的阶数
hold on;
for j = 1:2 % 绘制散点图和多项式拟合曲线的循环
    subplot(2, 1, j);  % 创建子图布局，这里是两行一列，选择第 j 个子图
    scatter(x, y, 'MarkerEdgeColor', 'r'); % 绘制散点图
    xlim([0, 40]); 
    ylim([-2, 6]);  
    hold on;
    for i = 1:length(n) % 对每个多项式的阶数进行循环
        p = polyfit(x, y, n(i)); % 使用 polyfit 函数进行多项式拟合
        p1 = round(p * 10^8) / 10^8; % 将拟合得到的多项式系数四舍五入到小数点后八位
        xi = 0:0.05:40;
        yi = polyval(p, xi); % 计算拟合曲线的 y 坐标点
        plot(xi, yi, '-b');
    end
end
hold off;

图像结果：

5- 干扰下的多项式拟合

例：用给定的多项式,如 $y=x^3-6x^2+5x-3$ ,产生一组数据(xi,yi,i=1,2,…,n),再在yi上添加随机干扰(可用rand产生(0,1)均匀分布随机数,或用randn产生N(0,1)分布随机数),然后用xi和添加了随机干扰的yi作3次多项式拟合,与原系数比较.如果作2或4次多项式拟合,结果如何?

三次多项式拟合，取41个点，命令行输入：

y = @(x) x.^3 - 6*x.^2 + 5*x -3;
xi = -10:0.5:10;
y1 = y(xi);
y2 = y(xi) + rand(1,41);
p = polyfit(xi,y2,3)

得到系数p（保留六位有效数字）：

1.00035, -5.99948, 4.97683, -2.54198

改为取11个点时得到系数p：

0.999004, -6.00018, 5.06891, -2.48958

取41个点和11个点时，拟合得到三次多项式系数在非常数项中与原系数几乎一致，误差极小，只有常数项大约比原常数项大接近0.5，这是因为生成了0到1间的随机数添加在yi上，其平均数为0.5，常数项减去0.5即接近原多项式。

二次多项式拟合时：

y = @(x) x.^3 - 6*x.^2 + 5*x -3;
xi = -10:0.5:10; % 2.改为“xi = -30:1.5:30”
y1 = y(xi);
y2 = y(xi) + rand(1,41);
p = polyfit(xi,y2,2)

绘图可得到

2次多项式拟合时，2次系数和常数项系数与3次多项式拟合得到结果类似，即2次系数与原多项式2次系数类似，常数项约为原常数项加0.5，而一次项不准确，随着所取原函数局部特性变化而变化，有较大差异。

四次多项式：

y = @(x) x.^3 - 6*x.^2 + 5*x -3;
xi = -10:0.5:10;
y1 = y(xi);
y2 = y(xi) + rand(1,41);
p = polyfit(xi,y2,4)

得到系数p（保留六位有效数字）：

0.000018, 1.00041, -6.00119, 4.98226, -2.56249

四次系数接近于0，剩下的部分与三次拟合的结果相似。

即大于或等于三次多项式拟合时都能得到较好的结果。

6- 非线性插值拟合：“lsqcurvefit”

（上一篇的例题2）

命令行输入：


Re = [2, 4, 8, 16, 32, 64, 128, 256, 512];
K = [56.3615,28.2406,14.2390,7.34264,4.0626,2.6171,2.0271,1.7572,1.5907];
g = @(y, Re) (y(1) + y(2) ./ Re); % 定义函数g，它是一个匿名函数，接受y和Re作为输入
f = @(x, Re) ((x(1).^20) + (x(2) ./ Re).^20).^(1 ./ 20); % 定义函数f，它是一个匿名函数，接受x和Re作为输入
x = [10, 10];
y = [10, 10]; % 定义初始猜测值x和y
b = lsqcurvefit(g, y, Re, K); % 使用lsqcurvefit函数拟合g函数，以找到最佳拟合参数b
a = lsqcurvefit(f, x, Re, K); % 使用lsqcurvefit函数拟合f函数，以找到最佳拟合参数a
scatter(Re, K, 'r'); % 绘制K关于Re的散点图
x1 = 2:0.1:512;
x2 = 2:10:512;
hold on
z1 = g(b, x1);% 计算g函数在b参数和x1上的值
plot(x1, z1, '-g'); % 绘制z1关于x1的曲线，使用绿色实线
z2 = f(a, x2); % 计算f函数在a参数和x2上的值
plot(x2, z2, '.k'); % 绘制z2关于x2的点图，使用黑色点
xlabel('Re');
ylabel('K');

输出图像：

7- 数值微积分：“gradient”、“quad”、“dblquad”、“integral2”、“int”

例1：已知参数方程 $\left\{\begin{matrix} x=ln(cost )\\ y=cost-tsint \end{matrix}\right.$ , 0<t<1.5, 试取t的步长0.01, 求 $\frac{dy}{dx}$ 和 $\frac{dy}{dx}|_{x=-1}$ 的数值解。

命令行输入：

t = 0:0.01:1.5;
x = log(cos(t));
y = cos(t) - t.*sin(t);
f = gradient(y,x)

可以绘制图像：

橙色线为y(x)，蓝色线为dy/dx(x)，横坐标x，纵坐标y

例2：用积分法计算下列椭圆的周长 $x^2/4+y^2/9=1$

命令行输入：

f = @(theta) sqrt(4*cos(theta).^2) + 9*(sin(theta).^2));
quad(f,0,pi*2)

输出：

15.865439474091897

x=2sin(theta),y=3cos(theta)，变量代换后，利用ds=sqrt(dx^2+dy^2)并且对ds积分可以得到椭圆的周长。

例3：求函数 $z=xe^{-x^2-y^2}$ ( -1<x<1, 0<y<2)构成曲面的面积。

三种方法，命令行输入：

% method 1: integral2
f = @(x,y) (x.*exp(-x.^2 - y.^2));
Q1 = integral2(f,-1,1,0,2)

% method 2: dblquad
f = @(x,y) (x.*exp(-x.^2 - y.^2));
Q2 = dblquad(f,-1,1,0,2)

% method 3: int*2
syms x y
Q3 = int(int(x*exp(-x^2-y^2),x,-1,1),y,0,2)

输出：

% method 1: integral2
Q1 = 1.6263e-17

% method 2: dblquad
Q2 = -6.7720e-18

% method 3: int*2
Q3 = 0

可以使用dblquad二重积分、使用integral2二重积分，或者利用int嵌套实现。MATLAB使用文档中不推荐使用dblquad函数。

从结果上来看，定义符号对象x和y，int嵌套后得到的结果就是0，而dblquad和integral2得到的结果也是0，但其实是近似0（误差极小忽略不计），这说明他们的计算路径不同。

【本篇完】

HelloFeynman

关注

32
点赞
踩
14

收藏

觉得还不错? 一键收藏
0
评论
MATLAB基础运算【6】：插值拟合与数值微积分

拉格朗日插值拟合：“自定义lagrange”；线性插值拟合：“interp1”；三次样条插值拟合：“spline”；多项式拟合曲线：“polyval”、“polyfit”；干扰下的多项式拟合；非线性插值拟合：“lsqcurvefit”；数值微积分：“gradient”、“quad”、“dblquad”、“integral2”、“int”
复制链接

扫一扫