33、MATLAB 数值技术：拟合、微分与积分的综合应用

MATLAB 数值技术：拟合、微分与积分的综合应用

1. 结果验证与分析

在进行数值计算后，我们需要对结果进行验证。将我们计算得到的结果与文献中报道的结果进行比较，发现二者接近但并不完全相同。例如，对于热容 (C_p) 的计算，我们拟合的结果为 (C_p = 2.737 \times 10^{-10}T^3 - 10.63 \times 10^{-7}T^2 + 1.552 \times 10^{-3}T + 4.683 \times 10^{-1})，而文献中的结果是 (C_p = 1.698 \times 10^{-10}T^3 - 7.957 \times 10^{-7}T^2 + 1.359 \times 10^{-3}T + 5.059 \times 10^{-1})。这并不令人意外，因为我们建模时使用的数据点数量有限，而文献中的模型使用了更多的数据，所以可能更准确。

2. 使用交互式拟合工具

MATLAB 提供了交互式绘图工具，无需使用命令窗口即可对绘图进行注释，还包含基本曲线拟合、复杂曲线拟合和统计工具。

2.1 基本拟合工具

要使用基本拟合工具，首先需要创建一个图形。以下是具体步骤：

x = 0:5;
y = [0,20,60,68,77,110];
plot(x,y,'o');
axis([−1,7,−20,120]);

这些命令会生成一个包含样本数据的图形。要激活曲线拟合工具，在图形的菜单栏中选择“工具” -> “基本拟合”，基本拟合窗口将出现在绘图上方。通过勾选“线性”、“三次”和“显示方程”，可以生成相应的拟合曲线。勾选“绘制残差”框会生成第二个绘图，显示每个数据点与计算线的距离。

在基本拟合窗口的右下角有一个箭头按钮，点击两次可打开窗口的其余部分。窗口的中间面板显示曲线拟合的结果，并提供将结果保存到工作区的选项；右侧面板允许根据中间面板显示的方程选择 (x) 值并计算 (y) 值。

除了基本拟合窗口，还可以从图形菜单栏访问数据统计窗口。选择“工具” -> “数据统计”，该窗口允许根据图形中的数据交互式计算统计函数，如均值和标准差，并将结果保存到工作区。

2.2 曲线拟合工具箱

除了基本拟合工具，MATLAB 还包含用于执行专业统计和数据拟合操作的工具箱。特别是曲线拟合工具箱提供了一个应用程序，允许拟合不仅仅是多项式的曲线。不过，该工具箱需要单独购买。

3. 差分与数值微分

3.1 diff 函数

函数 (y = f(x)) 的导数衡量了 (y) 随 (x) 的变化情况。如果能定义 (x) 和 (y) 的关系方程，可以使用符号工具箱中的函数求导数方程；如果只有数据，则可以通过 (y) 的变化量除以 (x) 的变化量来近似导数：
(\frac{dy}{dx} \approx \frac{\Delta y}{\Delta x} \approx \frac{y_2 - y_1}{x_2 - x_1})

例如，对于数据 (x = 0:5) 和 (y = [15, 10, 9, 6, 2, 0])，这种导数近似对应于连接数据的每个线段的斜率。MATLAB 有一个内置函数 diff ，可用于查找向量中元素值之间的差异，并可用于计算有序数据对的斜率。

delta_x = diff(x);
delta_y = diff(y);
slope = delta_y./delta_x;

使用 diff 函数返回的向量比输入向量短一个元素，因为是在计算差值。如果要绘制这些斜率与 (x) 的关系，最好创建一个条形图，因为变化率是不连续的。

x = x(:,1:5)+diff(x)/2;
bar(x,slope);

diff 函数也可用于在已知 (x) 和 (y) 关系的情况下数值近似导数。例如，对于函数 (y = x^2)，使用更多的 (x) 和 (y) 值可以使绘图更平滑。

x = −2:2;
y = x.^2;
big_x = −2:0.1:2;
big_y = big_x.^2;
plot(big_x,big_y,x,y,'−o');

通过不同数量的点来近似函数的斜率，使用的点越多，近似越准确。

3.2 前向、后向和中心差分技术

如果要在某一点近似导数，可以使用相邻点之间的斜率作为该点导数的近似值。

前向差分：
(\left(\frac{dy}{dx}\right) i = \frac{y {i + 1} - y_i}{x_{i + 1} - x_i})
可以通过 diff 函数实现：

x = linspace(0,pi/2,10);
y = sin(x);
dydx_analytical = cos(x);
dydx_approx = diff(y)./diff(x);
dydx_approx(length(x)) = NaN;
error_percentage = (dydx_analytical - dydx_approx)./dydx_analytical*100;
table = [x; dydx_analytical; dydx_approx; error_percentage];
disp('Forward Difference Approximation of the derivative of sin(x)');
disp(' x dy/dx dy/dx %error');
disp(' cos(x) forward approx.');
fprintf('%8.4f\t%8.4f\t%8.4f\t%8.4f\n',table);

后向差分与前向差分类似，只是将导数近似值分配给范围中的最后一个值。

中心差分：
(\left(\frac{dy}{dx}\right) i = \frac{y {i + 1} - y_{i - 1}}{x_{i + 1} - x_{i - 1}})
MATLAB 中的 gradient 函数使用前向差分技术处理数组中的第一个点，后向差分处理最后一个点，中心差分处理其余点来近似导数。

g = gradient(y,x);

gradient 函数也可用于二维数组来近似偏导数。以下是一个示例：

X = −2:.1:2;
Y = −2:.1:2;
[x,y] = meshgrid(X,Y);
z = 3*(1−x).^2.*exp(−(x.^2) − (y+1).^2) ...
    − 10*(x/5 − x.^3 − y.^5).*exp(−x.^2−y.^2) ...
    − 1/3*exp(−(x+1).^2 − y.^2);
surf(x,y,z);
[dzdx, dzdy] = gradient(z,X,Y);
contour(x,y,z);
hold on;
quiver(x,y,dzdx,dzdy);

4. 数值积分

积分通常被认为是曲线下的面积。可以通过将面积划分为矩形，然后对所有矩形的贡献求和来计算曲线下的面积，这就是梯形法则。

avg_y = y(1:5)+diff(y)/2;
sum(diff(x).*avg_y);

MATLAB 有一个内置函数 trapz 可以实现相同的结果：

trapz(x,y);

要近似由函数定义的曲线下的面积，可以创建一组有序的 (x - y) 对。增加 (x) 和 (y) 向量中的元素数量可以得到更好的近似。例如，对于函数 (y = x^2)，从 0 到 1 的积分可以这样计算：

x = 0:0.1:1;
y = x.^2;
trapz(x,y);

MATLAB 还提供了 integral 函数，可用于计算函数的积分，无需用户指定矩形的定义方式。以下是一个示例：

fun_handle = @(x) −x.^3 + 20*x.^2 −5;
fplot(fun_handle,[−5,25]);
integral(fun_handle,0,25);

5. 计算移动边界功

以活塞 - 气缸装置为例，使用 MATLAB 的数值积分技术 integral 来计算功。假设 (PV = nRT)，其中 (P) 是压力，(V) 是体积，(n) 是物质的量，(R) 是通用气体常数，(T) 是温度。

5.1 问题陈述

求活塞 - 气缸装置产生的功。

5.2 输入和输出

输入：
- (T = 300K)
- (n = 1kmol)
- (R = 8.314kJ/kmolK)
- (V_1 = 1m^3)
- (V_2 = 5m^3)

输出：活塞 - 气缸装置所做的功。

5.3 手动计算示例

由 (PV = nRT) 可得 (P = \frac{nRT}{V})，则功 (W = \int PdV = nRT \int \frac{dV}{V} = nRT \ln(\frac{V_2}{V_1}))。代入数值可得 (W = 4014kJ)，由于功为正，说明是系统对外做功。

5.4 MATLAB 解决方案

%Example 13.5
%Calculating boundary work, using MATLAB's quadrature function
clear, clc;
%Define constants
n = 1;      % number of moles of gas
R = 8.314;  % universal gas constant
T = 300;    % Temperature, in K
%Define an anonymous function for P
P = @(V) n*R*T./V;
%Use integral to evaluate the integral
integral(P,1,5);

运行结果与手动计算结果相同。

5.5 结果验证

将 MATLAB 计算结果与手动计算结果进行比较，结果一致。同时，使用符号工具箱求解也有助于验证结果。之所以需要两种 MATLAB 解决方案，是因为有些问题无法使用符号工具解决，而有些问题（存在奇点的问题）不适合使用数值方法。

以下是一个总结流程图：

graph LR
    A[结果验证] --> B[交互式拟合工具]
    B --> B1[基本拟合工具]
    B --> B2[曲线拟合工具箱]
    B --> C[差分与数值微分]
    C --> C1[diff函数]
    C --> C2[前向、后向和中心差分技术]
    C --> D[数值积分]
    D --> D1[梯形法则]
    D --> D2[trapz函数]
    D --> D3[integral函数]
    D --> E[计算移动边界功]
    E --> E1[问题陈述]
    E --> E2[输入和输出]
    E --> E3[手动计算示例]
    E --> E4[MATLAB解决方案]
    E --> E5[结果验证]

通过以上内容，我们详细介绍了 MATLAB 中的数值技术，包括拟合、微分和积分，以及如何应用这些技术解决实际问题。希望这些内容对大家在 MATLAB 编程和数值计算方面有所帮助。

MATLAB 数值技术：拟合、微分与积分的综合应用

6. 实践练习巩固知识

为了更好地掌握上述数值技术，下面提供一些实践练习，通过实际操作加深对这些方法的理解和运用。

6.1 数值微分练习

• 问题 1 ：考虑方程 (y = x^3 + 2x^2 - x + 3)，定义 (x) 向量从 (-5) 到 (+5)，使用 diff 函数以向前差分的方法近似 (y) 关于 (x) 的导数。解析导数为 (\frac{dy}{dx} = y’ = 3x^2 + 4x - 1)。
  matlab x = -5:5; y = x.^3 + 2*x.^2 - x + 3; dydx_approx = diff(y)./diff(x); dydx_analytical = 3*x.^2 + 4*x - 1; % 由于向前差分少一个值，添加 NaN 使长度一致 dydx_approx = [dydx_approx, NaN]; % 比较结果 table = [x; dydx_analytical; dydx_approx]; disp('Forward Difference Approximation of the derivative of y = x^3 + 2x^2 - x + 3'); disp(' x dy/dx_analytical dy/dx_approx'); fprintf('%8.4f\t%8.4f\t%8.4f\n', table);
• 问题 2 ：对以下函数及其导数重复上述操作：
  | 函数 | 导数 |
  | ---- | ---- |
  | (y = \sin(x)) | (\frac{dy}{dx} = \cos(x)) |
  | (y = x^5 - 1) | (\frac{dy}{dx} = 5x^4) |
  | (y = 5xe^x) | (\frac{dy}{dx} = 5e^x + 5xe^x) |

以下是 (y = \sin(x)) 的示例代码：

x = -5:5;
y = sin(x);
dydx_approx = diff(y)./diff(x);
dydx_analytical = cos(x);
dydx_approx = [dydx_approx, NaN];
table = [x; dydx_analytical; dydx_approx];
disp('Forward Difference Approximation of the derivative of y = sin(x)');
disp(' x dy/dx_analytical dy/dx_approx');
fprintf('%8.4f\t%8.4f\t%8.4f\n', table);

• 问题 3 ：使用 gradient 函数求上述问题中的导数值。以 (y = x^3 + 2x^2 - x + 3) 为例：

x = -5:5;
y = x.^3 + 2*x.^2 - x + 3;
dydx_gradient = gradient(y, x);
disp('Derivative of y = x^3 + 2x^2 - x + 3 using gradient function');
disp(' x dydx_gradient');
fprintf('%8.4f\t%8.4f\n', [x; dydx_gradient]);

• 问题 4 ：绘制结果并比较两种方法。以 (y = x^3 + 2x^2 - x + 3) 为例：

x = -5:5;
y = x.^3 + 2*x.^2 - x + 3;
dydx_approx = diff(y)./diff(x);
dydx_approx = [dydx_approx, NaN];
dydx_gradient = gradient(y, x);
dydx_analytical = 3*x.^2 + 4*x - 1;

figure;
subplot(3,1,1);
plot(x, dydx_analytical, 'b', 'DisplayName', 'Analytical');
title('Analytical Derivative');
xlabel('x');
ylabel('dy/dx');
legend;

subplot(3,1,2);
plot(x, dydx_approx, 'r', 'DisplayName', 'Forward Difference');
title('Forward Difference Approximation');
xlabel('x');
ylabel('dy/dx');
legend;

subplot(3,1,3);
plot(x, dydx_gradient, 'g', 'DisplayName', 'Gradient Function');
title('Gradient Function Approximation');
xlabel('x');
ylabel('dy/dx');
legend;

6.2 数值积分练习

• 问题 1 ：考虑方程 (y = x^3 + 2x^2 - x + 3)
  • (a) 使用 trapz 函数估计 (y) 关于 (x) 从 (-1) 到 (1) 的积分，使用 11 个 (x) 值。

x = linspace(-1, 1, 11);
y = x.^3 + 2*x.^2 - x + 3;
integral_trapz = trapz(x, y);
disp(['Integral using trapz: ', num2str(integral_trapz)]);

- **(b)** 使用 `integral` 函数求 \(y\) 关于 \(x\) 从 \(-1\) 到 \(1\) 的积分。

fun_handle = @(x) x.^3 + 2*x.^2 - x + 3;
integral_integral = integral(fun_handle, -1, 1);
disp(['Integral using integral: ', num2str(integral_integral)]);

- **(c)** 与使用符号工具箱函数 `int` 的解析解进行比较。解析解为 \(\int_{a}^{b} (x^3 + 2x^2 - x + 3) dx = \frac{1}{4} (b^4 - a^4) + \frac{2}{3} (b^3 - a^3) - \frac{1}{2} (b^2 - a^2) + 3(b - a)\)，这里 \(a = -1\)，\(b = 1\)。

syms x;
a = -1;
b = 1;
integral_analytical = int(x^3 + 2*x^2 - x + 3, x, a, b);
disp(['Integral using symbolic toolbox: ', char(integral_analytical)]);

• 问题 2 ：对以下函数重复上述操作：
  | 函数 | 积分 |
  | ---- | ---- |
  | (y = \sin(x)) | (\int_{a}^{b} \sin(x)dx = \cos(a) - \cos(b)) |
  | (y = x^5 - 1) | (\int_{a}^{b} (x^5 - 1)dx = \frac{b^6 - a^6}{6} - (b - a)) |
  | (y = 5xe^x) | (\int_{a}^{b} (5xe^x)dx = (-5e^b + 5be^b) - (-5e^a + 5ae^a)) |

以 (y = \sin(x)) 为例，使用 trapz 函数的代码如下：

x = linspace(0, pi, 11);
y = sin(x);
integral_trapz = trapz(x, y);
disp(['Integral of sin(x) using trapz: ', num2str(integral_trapz)]);

7. 总结与应用拓展

通过上述内容，我们全面了解了 MATLAB 中的数值技术，包括结果验证、交互式拟合工具的使用、差分与数值微分以及数值积分等方面。这些技术在科学研究、工程计算、数据分析等领域都有广泛的应用。

例如，在物理实验中，我们可以使用数值微分来计算物理量的变化率，如速度、加速度等；在化学工程中，数值积分可用于计算反应热、物质的量等。在实际应用中，我们可以根据具体问题选择合适的方法和函数，同时要注意结果的验证和误差分析。

以下是一个应用拓展的流程图，展示了如何将这些数值技术应用到实际问题中：

graph LR
    A[实际问题] --> B[数据采集]
    B --> C[选择合适的数值技术]
    C --> C1{拟合?}
    C1 -- 是 --> D[使用拟合工具]
    C1 -- 否 --> C2{微分?}
    C2 -- 是 --> E[使用差分或 gradient 函数]
    C2 -- 否 --> C3{积分?}
    C3 -- 是 --> F[使用 trapz 或 integral 函数]
    D --> G[结果分析与验证]
    E --> G
    F --> G
    G --> H[应用结果解决问题]

总之，MATLAB 的数值技术为我们解决各种复杂的数值问题提供了强大的工具。通过不断的实践和应用，我们可以更好地掌握这些技术，提高解决实际问题的能力。希望大家在今后的学习和工作中能够灵活运用这些方法，取得更好的成果。