“MATLAB程序设计（王赫然）”阅读笔记（二）——数学计算

代码下载

1.初等数学

1.1离散数学：

%分解质因数
primeFactor = factor(20)
%最大公约数
greatestCommonDivisor = gcd(20,15)
%最小公倍数
leastCommonMultiple = lcm(20,15)
%确定哪些数组元素为质数
isprime = isprime([2 3 0 6 10])
%小于等于输入值的质数
PrimeNumbers = primes(20)
%所有可能的排列
permulations = perms([1 2 3])
%阶乘
inputFactorial = factorial(5)
%有理分式近似值
rationalFraction = rat(pi,1e-7)

1.2多项式

%创建上述多项式
p = [1 -10 35 -50 24]
%计算x=2时多项式的值
polyValue = polyval(p,2)
%同时计算x为1到5的值
polyValue2 = polyval(p,1:5)
%多项式作图
x = 0.5:0.1:4.5;
plot(x,polyval(p,x));yline(0);
%求多项式的根
r = roots(p)
%由根反求多项式
pCalc = poly(r)
%显示多项式,poly2sym本质就是将多项式向量转换为符号矩阵
poly2sym(p)

2.线性代数

关于矩阵的基本运算，对于符号阵和数值阵都是通用的

syms a b c d
x = [a b;c d]
diagX = diag(x)
diagMatrixX = diag(diag(x))
traceX = trace(x)
rankX = rank(x)
detX = det(x)

syms a b c
x = [a b;c 0]
y = [1 2;3 0]
% 三角分解
% A=LU
% A=P'LU
[Lx,Ux] = lu(x)
[Lu,Uy] = lu(y)
[Ly,Uy,Py] = lu(y)
% 特征分解 A=VDV-1
[Vx,Dx] = eig(x)
[Vy,Dy] = eig(y)

3.微积分

1.极限

syms x a
f = (1+a/x)^x;
limF1 = limit(f,x,1)
limF2 = limit(f,x,0,'right')
limF3 = limit(f,x,-inf)
% 作a=1时的图像
fplot(@(x)(1+1/x)^x,[-50 50]);
yline(exp(1));%渐近线

对多元函数求极限，仍然是用limit()函数，

limit(limit(f,x,x0),y,y0)
%等价于
limit(limit(f,y,y0),x,x0)

2.导数
用diff函数，diff(f,x)表示f对x求导，diff(f,x,2)表示f对x求二阶导
3.积分

syms x a 
f = a*x^2
intf1 = int(f,x)
intf2 = int(int(f,x),x)
intf3 = int(f,x,0,1)

int对于不定积分与定积分通用，多次积分需要嵌套实现，不定积分结果省略了C。
4.泰勒展开

syms x
f = sin(x)
f1 = taylor(f,x,0,'order',3)
f2 = taylor(f,x,0) % 等效于6阶
% 作图
fplot(f,'LineWidth',2); hold on;
fplot(f1);
fplot(f2);

5.傅里叶展开

syms x
f = abs(x)/x
p = 10;
L = pi;
fs = int(f,x,-L,L)/2/L;
for n=1:p
    a = int(f*cos(n*pi*x/L),x,-L,L)/L;
    b = int(f*sin(n*pi*x/L),x,-L,L)/L;
    fs = fs+a*cos(n*pi*x/L)+b*sin(n*pi*x/L);
end
fs %显示傅里叶展开结果
fplot(f);hold on;
fplot(fs);

4.插值与拟合

1.一维插值

x = 0:pi/4:2*pi;
v = sin(x);
xq = 0:pi/16:2*pi;
ax(1) = subplot(121);
vq1 = interp1(x,v,xq);
plot(x,v,'o',xq,vq1,':.');
title('(Default) Linear Interpolation');
% 三次样条插值
ax(2) = subplot(122);
vq2 = interp1(x,v,xq,'spline');
plot(x,v,'o',xq,vq2,':.');
title('Spline Interpolation');
axis(ax,"square")

2.二维网格数据插值
调用格式：interp2(x0,y0,z0,x1,y1,'method')

[X,Y] = meshgrid(-3:3);
V = peaks(X,Y);
subplot(2,2,1);
surf(X,Y,V)
title('Original Sampling');
[Xq,Yq] = meshgrid(-3:0.2:3);
%
Vq = interp2(X,Y,V,Xq,Yq);
subplot(2,2,2);
surf(Xq,Yq,Vq);
title('Linear');
%
Vq = interp2(X,Y,V,Xq,Yq,'cubic');
subplot(2,2,3);
surf(Xq,Yq,Vq);
title('Cubic');
%
Vq = interp2(X,Y,V,Xq,Yq,'spline');
subplot(2,2,4);
surf(Xq,Yq,Vq);
title('Spline');

对于非网格数据，MATLAB提供了一个面向更为一般数据的插值函数griddata()函数：
griddata(x0,y0,z0,x,y,'method')

x = -3+6*rand(30,1);
y = -3+6*rand(30,1);
v = sin(x).^4.*cos(y);%构造已知数据
[xq,yq] = meshgrid(-3:0.1:3);%构造插值网
subplot(2,2,1);
z1 = griddata(x,y,v,xq,yq,'nearest');
plot3(x,y,v,'o');hold on
mesh(xq,yq,z1)
title('Nearest Neighbor')
subplot(2,2,2);
z2 = griddata(x,y,v,xq,yq,'linear');
plot3(x,y,v,'o');mesh(xq,yq,z2)
title('Linear')
subplot(2,2,3);
z3 = griddata(x,y,v,xq,yq,'natural');
plot3(x,y,v,'o');mesh(xq,yq,z3)
title('Natural Neighbor')
subplot(2,2,4);
z4 = griddata(x,y,v,xq,yq,'cubic');
plot3(x,y,v,'o');mesh(xq,yq,z4)
title('Cubic')

3.多项式拟合
拟合行为需要样本点数据作为输入

% 准备样本数据并作图
x = linspace(0,1,5);
y = 1./(1+x);
plot(x,y,'o');hold on
fplot(@(x) 1./(1+x),[-0.5 2],'--')
p = polyfit(x,y,4)
fitFunc = vpa(poly2sym(p),3)%显示便于观看的多项式
fplot(poly2sym(p),[-0.5 2])
legend('Sample','Origin','Fit')

5.代数方程与优化

1.代数方程

syms x
eqn = x^2+x-1 % 方程定义
solveSym1 = solve(eqn,x) %解析解
solveNum1 = vpasolve(eqn) %数值解

2.无约束优化
注意设置函数时，要将多元自变量写为x(1)、x(2)的形式。
优化工具箱提供的fminunc()优化速度更快。
求最大值只需要前面乘以一个负号。

[x,y] = meshgrid(-2:0.1:2,0:0.1:3);
z = 100*(y-x.^2).^2+(1-x).^2;
surf(x,y,z);
shading interp;
fun = @(x) 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
x0 = [-1.2,1];
options = optimset('Plotfcns',@optimplotfval);
x1 = fminsearch(fun,x0,options)
x2 = fminunc(fun,x0,options)

3.线性规划
$$\min {\mathbf{f}}^{\mathrm{T}}\mathbf{x}\text{ s. }\mathrm{t}.\{\begin{array}{l}\mathbf{A}\mathbf{x}⩽\mathbf{b}\\ {\mathbf{A}}_{\mathrm{e}\mathrm{q}}\mathbf{x}={\mathbf{b}}_{\mathrm{e}\mathrm{q}}\\ \mathbf{l}\mathbf{b}⩽\mathbf{x}⩽\mathbf{u}\mathbf{b}\end{array}$$
$$\begin{array}{l}f(x)=-2x_1-x_2\text{ 其中: }2x_2⩽7;x_1+4x_2⩽9;x_1⩽4;x_2⩽8\\ \text{ 求: }\min f(x)\end{array}$$

f = [-2 -1]';% 函数定义
A = [0 2;1 4];b = [7;9];
Aeq = []; beq = [];
lb = []; ub = [4;8];
[x,fval] = linprog(f,A,b,Aeq,beq,lb,ub)

4.非线性规划
$$\min f(\mathbf{x})\text{ s.t. }\{\begin{array}{l}\mathbf{A}\mathbf{x}⩽\mathbf{b}\\ {\mathbf{A}}_{\mathrm{e}\mathrm{q}}\mathbf{x}={\mathbf{b}}_{\mathrm{e}\mathrm{q}}\\ \mathbf{l}\mathbf{b}⩽\mathbf{x}⩽\mathbf{u}b\\ \mathbf{C}(\mathbf{x})⩽0\\ {\mathbf{C}}_{\mathrm{e}\mathrm{q}}(\mathbf{x})=0\end{array}$$
一般用边界约束优化函数fmincon来解决。
$$f(x)=-x_1{x}_2{x}_3\text{ 其中: }0⩽x_1+2x_2+4x_3⩽12\text{ 求 }\min f(x)$$

func = @(x) -x(1)*x(2)*x(3);
x0 = [1;1;1];% 起始值
A = [-1 -2 -4;1 2 4];
b = [0;12];
[x1,va11] = fmincon(func,x0,A,b)

5.最大值最小化问题
用[x,fval]=fminimax(fun,x0,A,b,Aeq,beq,lb,ub)格式的函数。

6.微分方程

1.常微分方程(ODE)解析解
无初始条件微分方程

syms a t y(t) z(t)
eqn = diff(y,t,2) == a*y %符号方程
ySol(t) = dsolve(eqn)

有初始条件的微分方程

syms a t y(t) z(t)
eqn = diff(y,t) == a*y
cond = y(0) == 5;
ySol(t) = dsolve(eqn,cond)

微分方程组

syms a t y(t) z(t)
eqns = [diff(y,t)==z,diff(z,t)==-y]
sol = dsolve(eqns);
solY = sol.y,solZ = sol.z

2.常微分方程数值解
用ode45()函数
3.微分方程Simulink求解
4.偏微分方程工具箱PDE
pdetool命令打开PDE工具箱。

7.概率统计

1.概率分布，连续随机变量的概率分布函数F(x)，又叫累计分布函数，即CDF，概率密度函数，简称PDF。
用法：

x = 0:20; lamda = 10;
pdf_P = pdf('Poisson',x,lamda);
cdf_P = cdf('Poisson',x,lamda);
ax(1) = subplot(2,2,1);stem(x,pdf_P),line(x,pdf_P);
%离散点用stem画针状图，连续点用plot画连续图
title('Poisson PDF');
ax(2) = subplot(2,2,2);plot(x,cdf_P);
title('Poisson CDF');

工具箱：命令行输入：disttool

2.伪随机数
产生均匀分布随机数的rand()
产生均匀分布随机整数的randi()
产生标准正太分布随机数的randn()
统计和机器学习工具箱SMLT中还提供了更为通用和强大的随机数产生函数random()，代码形式：R=random('name',A,B,C,D,[size])
3.统计量分析
mean()均值函数、var()方差函数、std()标准差函数、moment()中心矩函数
以上函数的输入参数不仅可以是向量，还可以是矩阵，相当于对矩阵中所有列进行计算，因而输出的是一个行向量，如果需要对矩阵中所有元素进行计算，则可以采用mean(x(:))

x = random('Normal',2,4,[1e4,1]);
meanX = mean(x) %均值
s2X = var(x) %方差
sX = std(x) %标准差
% 构造概率分布对象
pd = makedist('Normal','mu',2,'sigma',4);
meanX2 = pd.mean %均值
meanX2_mu = pd.mu %均值
s2X2 = pd.var %方差
s2X2 = pd.std %标准差
sX2_sigma = pd.sigma %标准差
A1 = sum(x.^1)/length(x) %一阶原点矩
A2 = sum(x.^2)/length(x) %二阶原点矩
B1 = moment(x,1) %一阶中心距
B2 = moment(x,2) %二阶中心矩

4.参数估计
通用参数估计函数pd=fitdist(x,distname)
distributionFitter——参数估计器APP
5.假设检验
对于疑似符合正态分布的数据，当已知标准差时，可以采用Z检验，函数为ztest()，若不已知标准差时，可以使用T检验。h为检验结论，为0时表示不拒绝假设，也就是接受假设；p为该检验的显著性水平；ci为置信区间。

x = random('Normal',2,4,[1e4,1]);
sX = std(x);
%Z检验—已知标准差
[hZ,pZ,ciZ,~] = ztest(x,2,sX)
%T检验—不知标准差
[hT,pT,ciT,~] = ttest(x,2)

.