系列更新MATLAB各种运算和操作的说明。
本篇主要包括插值拟合/数值微积分相关内容(可利用侧边目录直接定位):
- 拉格朗日插值拟合:“自定义lagrange”
- 线性插值拟合:“interp1”
- 三次样条插值拟合:“spline”
- 多项式拟合曲线:“polyval”、“polyfit”
- 干扰下的多项式拟合
- 非线性插值拟合:“lsqcurvefit”
- 数值微积分:“gradient”、“quad”、“dblquad”、“integral2”、“int”
前文:
例:对"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- 干扰下的多项式拟合
例:用给定的多项式,如
,产生一组数据(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:已知参数方程
, 0<t<1.5, 试取t的步长0.01, 求
和
的数值解。
命令行输入:
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:用积分法计算下列椭圆的周长
![]()
命令行输入:
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:求函数
( -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(误差极小忽略不计),这说明他们的计算路径不同。
【本篇完】