第四章 MATLAB的数学运算

4.1 多项式与插值

数学运算时MATLAB的核心，其他领域的应用都以这些数学运算为基础，如矩阵运算、代数运算、解方程等等。MATLAB的强大之处正是在于其对几乎所有数学问题的解决方案。本章将介绍MATLAB的部分基础数学功能，包括多项式，线性插值，傅里叶变换和微分方程等。这些数学知识在今后的学习中会经常用到，非常重要。

4.1.1 多项式的表示

多项式在数学中有着极为重要的作用，而多项式的运算也是工程和应用中经常会遇到的问题。
在MATLAB中多项式用一个行向量表示，向量中的元素为该多项式的系数，按照降序排列。例如，多项式$4x^3+3x^2+6x^1+9$可以表示为向量$P=[4,3,6,9]$。用户可以用创建向量的方式创建多项式，再利用poly2sym函数将其显示为多项式。

例子：

clear all; clc;

P = [4, 3, 6, 9];
% 默认变量为x
y = poly2sym(P)

4.1.2 多项式的四则运算

使用向量来表示多项式使得多项式的四则运算变得简单，因为可以直接转换为向量的运算。
多项式的加减即为对应系数项的加减，因此可以通过向量的加减来完成。进行加减运算的两个向量需要有相同的长度，因此对低阶多项式需要补零
多项式的乘法和除法实际上是多项式系数之间的卷积和反卷积运算。MATLAB有相应的卷积函数conv和deconv。
例子：

clear all; clc;
% polynomial operation
p1=[1 2 1];                         		%定义多项式
p2=[1 1];
length_of_p1=length(p1);
length_of_p2=length(p2);
if length_of_p1 == length_of_p2     		%判断两个多项式的长度是否相等
    p1_plus_p2 =p1+p2;              		%多项式相加
    p1_minus_p2=p1-p2;              	%多项式相减
elseif length_of_p1 < length_of_p2
        temp_p1=[zeros(length_of_p2-length_of_p1) p1];
        p1_plus_p2 =temp_p1+p2;
        p1_minus_p2=temp_p1-p2;
else
    temp_p2=[zeros(length_of_p1-length_of_p2) p2];
    p1_plus_p2 =p1+temp_p2;
    p1_minus_p2=p1-temp_p2;
end
p1_multiply_p2=conv(p1,p2);         		%多项式相乘
p1_divide_p2  =deconv(p1,p2);       		%多项式除法
p1=poly2sym(p1)                    		%显示多项式 p1
p2=poly2sym(p2)                    		%显示多项式 p2
p1_plus_p2 =poly2sym(p1_plus_p2)
p1_minus_p2=poly2sym(p1_minus_p2)
p1_multiply_p2=poly2sym(p1_multiply_p2)
p1_divide_p2  =poly2sym(p1_divide_p2)

4.1.3 多项式的其他运算

表4-1 多项式运算函数

这些都是多项式运算中常见的函数。下面分别简要举例介绍一下。
roots() 和 poly()

clear all; clc;
% roots求多项式的根--注意与root的区别
A = [1 2 1];
yA = poly2sym(A)
a = roots(A)

%%
B = [1, 4, 5, 3, 6, 7];
yB = poly2sym(B)
b = roots(B)

%% 
% poly求具有特定根的多项式
A2 = poly(a)

B2 = poly(b)

%%
% poly求矩阵的特征多项式
P = magic(4);
p = poly(P);
y = poly2sym(p)
x = roots(p)

polyvar多项式求值

clear all; clc;

% 求多项式的值
p = [4 3 2 1];
polyval(p, 4)

x = roots(p)
polyval(p, x)

%%
% 尝试用另一种方法？
syms x
y = poly2sym(p);
disp(y)
r = root(y, x);
vr = vpa(r);
disp(vr)
n = polyval(p, double(vr))

polyder多项式求导

clear all; clc;

a = [1, 2, 1];
% x^2 + 2x + 1
ya = poly2sym(a)

b = [1, -1];
% x-1
yb = poly2sym(b)

% da
p1 = polyder(a);
dya = poly2sym(p1)

% d(ab)
p2 = polyder(a,b);
dy2 = poly2sym(p2)

% d(a/b)
[q, d] = polyder(a, b);
dy3 = poly2sym(q) / poly2sym(d)

polyfit多项式拟合
曲线拟合是工程中经常要用到的技术之一。MATLAB提供了曲线拟合工具箱和多项式拟合函数。函数polyfit给出了使用最小二乘法拟合的多项式。函数调用格式为p = polyfit(x, y, n)，其中x, y分别为要拟合数据的坐标，应该是相同长度的向量，n则是用于拟合的多项式的次数。

clear all; clc;
x = 0 : pi/10 : 2*pi;
y = sin(x);

z = polyfit(x, y, 3);

% 画图
plot(x, y, 'r*')
hold on
f = poly2sym(z);
ezplot(f, [0, 2*pi])

4.1.4 数据插值

另一种工程上常见的问题是数据插值，即根据已知数据推断未知数据。
插值运算时根据已有数据的分布规律，找到一个可以连接起已知各点的函数表达式，并用这一函数表达式来预测处于已有数据两点之间任意位置的数据。
MATLAB提供了对数组的任意一维数据进行插值的工具，这些工具大多需要用到多维数组的操作。本节主要介绍对一维数据插值。
一维数据插值在曲线拟合和数据分析中具有重要的地位。在MATLAB中，一维数据插值主要由函数interp1来实现（即interpolation1的缩写）。该函数的调用格式为yi = interp1(x, y, xi, method)，其中x,y 为所采用的数据的坐标，均是一维向量；xi为待插值的位置，method为采用的插值方法，返回插值结果。而该函数有四种method，不同method得到的结果会有些许差异。

表4-2 插值函数中可选的方法

例子：

% Interpolation using the four methods
x=[0 3 5 7 9 11 12 13 14 15];
y=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.4 1.6];
length_of_x=length(x);
scalar_x=[x(1):0.5:x(length_of_x)];
length_of_sx=length(scalar_x);
for i=1:length_of_sx
    y_nearest(i)=interp1(x,y,scalar_x(i),'nearest');
    y_linear(i) =interp1(x,y,scalar_x(i),'linear');
    y_spline(i) =interp1(x,y,scalar_x(i),'spline');
    y_cubic(i) =interp1(x,y,scalar_x(i),'pchip');
end
subplot(2,2,1),plot(x,y,'*'),hold on,plot(scalar_x,y_nearest),title('method=nearest');
subplot(2,2,2),plot(x,y,'*'),hold on,plot(scalar_x,y_linear),title('method=linear');
subplot(2,2,3),plot(x,y,'*'),hold on,plot(scalar_x,y_spline),title('method=spline');
subplot(2,2,4),plot(x,y,'*'),hold on,plot(scalar_x,y_cubic),title('method= pchip ');

$\ast$最近邻插值将插值点xi的值设置为举例最近的点的对应值。
$\ast$线性插值用分段线性函数拟合已有数据，返回拟合函数在xi处的值。
$\ast$三次样条插值采用样条函数对数据拟合，并且任意两点之间的函数为三次函数。
$\ast$三次插值为一组方法，通过pchip函数对数据进行三次Hermite插值，这种方法可以保持数据的一致性和数据曲线形状。

可以看出，最近邻插值效果最差，其他三种方法差别不是很大，实际上三次插值的效果相对较好，但是运行效率更差。在实际应用中，需要根据实际情况综合考虑来选取合适的插值方法。

傅里叶插值法

MATLAB针对一维数据还有另一种方法，基于傅里叶变换的方法。这种方法将输入数据视为周期函数的采样数据，对数据进行傅里叶变换，然后对更多的点进行傅里叶变换。函数名为interpft（interpolation fourier transform），调用方式为y = interpft(x, n)，其中x为周期函数的均匀采样数据，n为待返回的数据个数。当我们对周期数据或近似周期数据进行插值时，该方法更加合适。

例子：

clear all;clc;

x = 0:pi/5:2*pi; % 10个数据
y = sin(x);
plot(x, y);

hold on
y1 = interpft(y, 20);

x1 = linspace