读书笔记《数学建模算法与应用》第4-6章

注：章节序号与书中保持一致，省略部分为基础概念

4. 图与网络模型及方法

4.2 最短路问题

4.2.1 两个指定顶点之间的最短路径

Dijkstra 算法 function[mydistance,mypath]=mydijkstra(a,sb,db);

• a 邻接矩阵
• a(i,j) i-j的距离，可以是有向的
• sb 起点的标号
• db 终点的标号
• mydistance 最短路距离
• mypath 最短路路径

4.2.3 每对顶点之间的最短路径

Dijkstra 算法是时间复杂度是 $O(n^3)$，第二种解决该问题的方法是Floyd 算法

Floyd 算法原理
最短路径-Floyd算法

Matlab 实现

% floyd.m
% 采用floyd算法计算图a中每对顶点最短路
% d是矩离矩阵
% r是路由矩阵
function [d,r]=floyd(a)
n=size(a,1);
% 初始化距离矩阵
d=a;
% 初始化路由矩阵
for i=1:n
    for j=1:n
        r(i,j)=j;
    end 
end 
r;

% Floyd算法开始
for k=1:n
    for i=1:n
        for j=1:n
            if d(i,k)+d(k,j)<d(i,j)
                d(i,j)=d(i,k)+d(k,j);
                r(i,j)=r(i,k);
            end 
        end 
    end
    k;
    d;
    r;
end
d
r

4.3 最小生成树问题

实际问题举例： 把无向图的每个顶点看作村庄，计划修建道路使得可以在所有村庄之间通行。把每个村庄之间修建道路的费用看作权值，那么我们就可以得到一个求解修建道路的最小费用的问题。

4.3.1 最小生成树

最小生成树特点：

• 无回路，且包含原图中的n-1条边。
• 包含原图中的全部顶点。
• 边的权重和在所有其他生成树中最小。
• 最小生成树存在，则该图一定连通。反过来一样，图连通，则最小生成树一定存在。

最小生成树：Prim算法
算法导论–最小生成树（Kruskal和Prim算法）

prim 算法构造最小生成树
操作方法： 从某一顶点出发，逐步构建，让一棵小树逐渐长大。

Kruskal 算法构造最小生成树
操作方法： 此算法可以称为“加边法”，初始最小生成树边数为0，每迭代一次就选择一条满足条件的最小代价边，加入到最小生成树的边集合里。

4.4 网络最大流问题

图割-最大流最小切割的最直白解读

实际问题举例： 将油从 s 运送到 t，中间的四个点为中转站，边的系数为管道容量。求 s 到 t 的最大流。

4.5 最小费用最大流问题

实际问题举例： 将油从 s 运送到 t，中间的四个点为中转站，边的系数为管道容量。求 s 到 t 的最大流。在该题干基础上，考虑不同管道上的运输费用也不相同，因为除了考虑输油管道的最大流外，还需要考虑输油管道输送最大流的最小费用。

4.6 Matlab 的图论工具箱

实例1： 求解 $v_1$ 到 $v_{1}1$ 的最短路径及长度

clc,clear
a(1,2)=2;a(1,3)=8;a(1,4)=1;
a(2,3)=6;a(2,5)=1;
a(3,4)=7;a(3,5)=5;a(3,6)=1;a(3,7)=2;
a(4,7)=9;
a(5,6)=3;a(5,8)=2;a(5,9)=9;
a(6,7)=4;a(6,9)=6;
a(7,9)=3;a(7,10)=1;
a(8,9)=7;a(8,11)=9;
a(9,11)=2;a(9,10)=1;
a(10,11)=4;
a=a'; % matlab工具箱要求数据为下三角矩阵
[i,j,v]=find(a);  % 找出a中非0元素的行和列，分别存储在i,j中，并将结果放在v中
b=sparse(i,j,v,11,11)  % 构造稀疏矩阵
[x,y,z]=graphshortestpath(b,1,11,'Directed',false);  % 声明是无向图'Directed',false（或写0）

[dist, path, pred]=graphshortestpath(G,S,T)

• G是稀疏矩阵，S是起点，T是终点。
• dist表示最短距离
• path表示最短距离经过的路径节点
• pred表示从S到每个节点的最短路径中，目标节点的先驱，即目标节点的前面一个节点。

实例2： （最短路算法）

• 用四维向量表示状态，(人，狼，羊，菜)，在此岸时取1，在对岸时取0。
• 穷举所有可行状态，(1,1,1,1)；(1,1,1,0) (1,1,0,1) (1,0,1,1) ；(1,0,1,0) (0,1,0,1)；(0,0,1,0) (0,0,0,1) (0,0,0,0)
• 构造赋权图，顶点为以上10种可行状态，状态间可以转换时边对应权重为1，反正则为 $\infty$
• 问题转化为，初始点(1,1,1,1)，结束点(0,0,0,0)的最短路径

本题难点在于邻接矩阵的表示，因为摆渡一次 就会改变现有的状态，为此再引入一个状态转移向量。

• 1表示过河，0表示未过河，状态转移有4种情况，人自己过河(1,0,0,0)；人带狼(1,1,0,0)；(1,0,1,0)；(1,0,0,1)
• 状态向量 与 转移向量运算：0+0=0；1+0=1；0+1=1；1+1=0
• 如果可行状态+转移向量=可行状态，那么这两种可行状态对应的顶点就存在边。

% 4.10
clc,clear
% 可行状态行
a=[1 1 1 1;1 1 1 0;1 1 0 1;1 0 1 1; 1 0 1 0;0 1 0 1;0 1 0 0;0 0 1 0;0 0 0 1;0 0 0 0];
% 转移向量行
b=[1 0 0 0;1 1 0 0;1 0 1 0;1 0 0 1];
w=zeros(10); % 邻接矩阵初始化
for i=1:9 % 循环某顶点
    for j=i+1:10 % 循环某顶点的下一个顶点
        for k=1:4
            % xor() 异或运算，参数同为0/不为0 返回0，参数有0有1时 返回1
            if findstr(xor(a(i,:),b(k,:)),a(j,:)) % 判断初始顶点+状态转移 是否等于 下一个顶点
                w(i,j)=1; % 是的话，初始顶点-下一个顶点的边权重=1
            end
        end
    end
end
w=w';  % 下三角矩阵
c=sparse(w);  % 构造稀疏矩阵
[x,y,z]=graphshortestpath(c,1,10,'Directed',0);
h=view(biograph(c,[],'ShowArrows','off','ShowWeights','off'));  % 画出无向图
Edges=getedgesbynodeid(h);  % 提取句柄h种的边集
set(Edged,'LineColor',[0 0 0]);  % 边为黑色
set(Edged,'LineWidth',1.5);  % 边宽度1.5

实例3： （注意是有向图）

clc,clear
a=zeros(7);
a(1,2)=4;a(1,3)=2;
a(2,3)=3;a(2,4)=2;a(2,5)=6;
a(3,4)=5;a(3,6)=4;
a(4,5)=2;a(4,6)=7;
a(5,6)=4;a(5,7)=8;
a(6,7)=3;
b=sparse(a); % 稀疏矩阵
[x,y,z]=graphshortestpath(b,1,7,'Directed',1,'Method','Bellman-Ford'); % 有向图；
view(biograph(b,[]))

实例4：

求总电缆长度最短的问题实际上就是求图 $G$ 的最小生成树

matlab的mandist函数

[Tree, pred] = graphminspantree(G,R)

• Tree：给出最后的树的边集答案
• pred：答案中每个边对应的权值
• G：稀疏矩阵
• R：根节点（可以不传）

clc,clear
x=[0 5 16 20 33 23 35 25 10];
y=[15 20 24 20 25 11 7 0 3];
xy=[x;y]; % 拼接成2*9的矩阵
% mandist（A,B)函数是用来求A中的每个行向量与B中的每个列向量的绝对距离,mandist(A) = mandist(A’,A)
d=mandist(xy);  
d=tril(d); % 截取下三角矩阵
b=sparse(d); % 转化为稀疏矩阵
[ST,pred]=graphminspantree(b,'Method','Kruskal'); % 调用最小生成树
st=full(ST); % 将最小生成树的稀疏矩阵转化为普通矩阵
% sum(st)为每列之和，sum(sum(st))返回元素之和
TreeLength=sum(sum(st)); % 求最小生成树的长度
view(biograph(ST,[],'ShowArrows','off'))

% 解出的树的边集
ST =
   (2,1)       10
   (4,2)       15
   (4,3)        8
   (5,4)       18
   (6,4)       12
   (7,6)       16
   (8,6)       13
   (9,8)       18

实例5：

Matlab 图论工具箱求解最大流命令，只能解决权重都为正值，且两个顶点之间不能有两条弧。本题中3-4点有两条弧，处理：

• 删掉4→3权重为2的弧
• 引入虚拟顶点9
• 设置为4→9权重为2，9→3权重为2

[MaxFlow, FlowMatrix, Cut] = graphmaxflow(G, SNode, TNode)

clc,clear,a=zeros(9);
a(1,2)=6;a(1,3)=4;a(1,4)=5;
a(2,3)=3;a(2,5)=9;a(2,6)=9;
a(3,4)=5;a(9,3)=2;a(3,5)=6;a(3,6)=7;a(3,7)=3;
a(4,7)=5;a(4,9)=2;
a(5,8)=12;
a(6,5)=8;a(6,8)=10;
a(7,6)=4;a(7,8)=15;
b=sparse(a);
[x,y,z]=graphmaxflow(b,1,8);

5. 插值与拟合

线性插值(Linear Interpolation)原理及使用

定义： 简单说，利用插值 可以通过函数在有限个点处的取值状况，估算出函数在其他未知点的近似值。
作用：

• 用于补齐缺失值

• 对数据进行放大或缩小
  
  最简单的线性插值问题，就是两点式，已知点(x0,y0)、(x1,y1)，试问在x处插值，y的值是多少？
  $$L_1(x)=\frac{x-x_1}{x_0-x_1}{y}_0+\frac{x-x_0}{x_1-x_0}{y}_1$$
  显然，$L_1(x_0)=y_0$, $L_1(x_1)=y_1$ 满足插值条件，所以$L_1(x)$就是线性插值函数

  如果，我们令 $l_0(x)=\frac{x-x_1}{x_0-x_1}{l}_1(x)=\frac{x-x_0}{x_1-x_0}$，那么可以得到$L_1(x)=y_0{l}_0(x)+y_1{l}_1(x)$，我们可以看到：
  $$\{\begin{array}{l}{\mathbf{l}}_{\mathbf{o}}\left({\mathbf{x}}_{\mathbf{o}}\right)=1,{\mathbf{l}}_{\mathbf{o}}\left({\mathbf{x}}_1\right)=0\\ {\mathbf{l}}_1\left({\mathbf{x}}_0\right)=0,{\mathbf{l}}_1\left({\mathbf{x}}_1\right)=1\end{array}$$
  于是，我们将 $l_0(x),l_1(x)$ 称为关于 $x_0,x_1$ 的线性插值基函数。总结一下：线性插值函数$L_1(X)$ 可以通过插值基函数 $l_0(x),l_1(x)$ 的线性组合得到，且组合系数 $y_0,y_1$ 恰为所给数据

5.1 插值方法

5.1.1 分段线性插值

分段越多，插值误差越小。因此将其推广到一般形式：

5.1.2 拉格朗日插值多项式

5.1.3 样条插值

三次样条插值法
为了使通过不同方法得到的插值函数的衔接点平滑，使每个分段函数都采用高次函数形式构造（三次样条差值 就是用 $x^3$ 形式构造）

5.1.4 Matlab 插值工具箱

• 一维插值函数 y=interpl(x0,y0,x,'method') 其中（1）x0 是单调的；（2）method 指定插值方法
• 三次样条插值

y=interpl(x0,y0,x,'spline'); % x0,y0是已知数据点；x是插值点；y是插值点的函数值
y=spline(x0,y0,x);
pp=csape(x0,y0.conda);
pp=csape(x0,y0.conda,valconda);y=fnval(pp,x);
% conda 指定插值的边界条件

x0=[0   3   5   7   9   11   12   13   14  15];
y0=[0  1.2  1.7  2.0  2.1  2.0  1.8  1.2   1.0  1.6];
x=0:0.1:15; % 需要得到x每改变0.1时的y坐标

% 分段线性插值
y1=interp1(x0,y0,x);

% 样条插值1
y2=interp1(x0,y0,x,'spline'); % 'linear'线性插值；'spline'逐段3次样条插值

% 样条插值2
pp1=csape(x0,y0); % pp返回样条函数表达式，默认边界条件是lagrange边界
y3=fnval(pp1,x);  % 必须通过调用fnval返回插值点的函数值

pp2=csape(x0,y0,'second'); % conds='second'指定插值的边界条件为二阶导数
y4=fnval(pp2,x);

% 绘图
[x',y1',y2',y3',y4'] % 将行数据转置
subplot(1,3,1) % 子图1
plot(x0,y0,'+',x,y1)
title('Piecewise linear')

subplot(1,3,2) % 子图2
plot(x0,y0,'+',x,y2)
title('Spline1')

subplot(1,3,3) % 子图3
plot(x0,y0,'+',x,y3)
title('Spline2')

% 求x=0处的曲线斜率
dx=diff(x); % diff求前后相邻元素之差
dy=diff(y3);
dy_dx=dy./dx;  % a./b就是a、baib中对应的每个元素相除,如果a、b是两个数，那么a./b就是普通的除法
dy_dx0=dy_dx(1)

% [13,15]范围内的y最小值
ytemp=y3(131:151);
ymin=min(ytemp);
index=find(y3==ymin);
xmin=x(index);
[xmin,ymin]

clc, clear
x0=0.15:0.01:0.18;
y0=[3.5	1.5	2.5	2.8];
pp=csape(x0,y0)   %默认的边界条件，Lagrange边界条件
format long g % 将数据显示为长整型科学计数
xishu=pp.coefs   %显示每个区间上三次多项式的系数
s=quadl(@(t)fnval(pp,t),0.15,0.18)  %求积分
format  %恢复短小数的显示格式

• 二位插值：若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的高程（节点值），为了画出较准确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）。

插值节点为网格节点

clear,clc
x=100:100:500; % x-m维
y=100:100:400; % y-n维
z=[636    697    624    478   450  
   698    712    630    478   420
   680    674    598    412   400
   662    626    552    334   310]; % z-n*m维
pp=csape({x,y},z') % 注意这里的z0需要是m*n维
xi=100:10:500;yi=100:10:400; % 插值点
cz=fnval(pp,{xi,yi});
[i,j]=find(cz==max(max(cz)))  %找最高点的地址，max(cz)按列求最大，max(max(z))整个矩阵最大
% a=find(z==max(max(z)))则返回该最大值出现的序列位数
% [i,j]=find(z==max(max(z)))则返回该最大值出现的行数，列数
x=xi(i),y=yi(j),zmax=cz(i,j)  %求最高点的坐标 

插值节点为散乱节点

clc, clear
x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];
y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xmm=minmax(x)  %求x的最小值和最大值
ymm=minmax(y)  %求y的最小值和最大值
xi=xmm(1):xmm(2);
yi=ymm(1):ymm(2);

% 立方插值
zi1=griddata(x,y,z,xi,yi','cubic'); % xi,yi给定的网格点的横坐标和纵坐标，zi1为网格(xi,yi)处的函数值
% 最近点插值
zi2=griddata(x,y,z,xi,yi','nearest'); 

zi=zi1;  %立方插值和最近点插值的混合插值的初始值
% isnan()返回一个与A相同维数的数组,若A的元素为NaN（非数值），在对应位置上返回逻辑1（真），否则返回逻辑0（假）
zi(isnan(zi1))=zi2(isnan(zi1))  %把立方插值中的不确定值换成最近点插值的结果
subplot(1,2,1), plot(x,y,'*')
subplot(1,2,2), mesh(xi,yi,zi) %mesh命令绘制双变量的三维图

5.2 曲线拟合的线性最小二乘法

5.2.1 线性最小二乘法

• 系数 $a_k$ 的确定，记$$\begin{array}{l}\mathbf{R}={\left[\begin{array}{ccc}{r}_1\left(x_1\right)& \cdots & r_m\left(x_1\right)\\ ⋮& ⋮& ⋮\\ r_1\left(x_n\right)& \cdots & r_m\left(x_n\right)\end{array}\right]}_{n\times m}\\ \mathbf{A}={\left[a_1,\cdots ,a_m\right]}^{\mathrm{T}},\mathbf{Y}={\left[y_1,\cdots ,y_n\right]}^{\mathrm{T}}\end{array}$$
• 函数 $r_k(x)$ 的选取，常用的曲线有：
  

5.2.2 最小二乘法的Matlab 实现

解方程组法

x=[19     25    31     38    44]';
y=[19.0   32.3   49.0   73.3   97.8]';
r=[ones(5,1),x.^2]; % r(x)
ab=r\y % 两个系数
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')

多项式拟合方法

% 根据已知点做散点图，观察到年生产利润几乎直线上升，因此可以用y=a1*x+a0作为拟合函数
x0=[1990  1991  1992  1993  1994  1995  1996];
y0=[70   122   144   152   174   196   202];
plot(x0,y0,'*')

x0=[1990  1991  1992  1993  1994  1995  1996];
y0=[70   122   144   152   174   196   202];
a=polyfit(x0,y0,1) % 1代表1次多项式
% 注意，a返回两个值，第一个值是最高项x的系数，最后的值是常数项
y97=polyval(a,1997)
y98=polyval(a,1998)

5.3 最小二乘优化

5.3.1 lsqlin 函数

x=[19     25    31     38    44]';
y=[19.0   32.3   49.0   73.3   97.8]';
r=[ones(5,1),x.^2];
ab=lsqlin(r,y)
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')

5.3.2 lsqcurvefit 函数

% 定义函数F(x,xdata)
function f=fun1(canshu,xdata);
f= exp(-canshu(1)*xdata(:,1)).*sin(canshu(2)*xdata(:,2))+xdata(:,3).^2;  
%其中canshu(1)=k1，canshu(2)=k2，注意函数中自变量的形式

clc, clear
a=textread('data1.txt');
y0=a(:,[2,7]); %提出因变量y的数据，第2和第7行
y0=nonzeros(y0); %去掉最后的零元素，且变成列向量
x0=[a(:,[3:5]);a([1:end-1],[8:10])]; %由分块矩阵构造因变量数据的2列矩阵
canshu0=rand(2,1); %拟合参数的初始值是任意取的
%非线性拟合的答案是不唯一的，下面给出拟合参数的上下界，
lb=zeros(2,1); %这里是随意给的拟合参数的下界,无下界时，默认值是空矩阵[]
ub=[20;2]; %这里是随意给的上界，无上界时，默认值是空矩阵[]
canshu=lsqcurvefit(@fun1,canshu0,x0,y0,lb,ub)  

5.4 曲线拟合与函数逼近

syms x
base=[1,x^2,x^4];
y1=base.'*base % 3*3矩阵
y2=cos(x)*base.'
r1=int(y1,-pi/2,pi/2) % 求定积分，上下限为-pi/2,pi/2
r2=int(y2,-pi/2,pi/2)
a=r1\r2
xishu1=double(a)  %符号数据转化成数值型数据
xishu2=vpa(a,6)   %把符号数据转化成小数型的符号数据

6. 微分方程建模

将实际问题转化为微分方程的定解问题，大体上可分为以下几步：

• 根据实际要求确定要研究的量（自变量、未知函数、必要的参数等）并确定坐标系。
• 找出这些量满足的基本规律。
• 运用规律列出方程和定解条件。

微分方程模型：

• 人口预测模型
• 捕食者猎物模型
• 种群相互竞争模型
• 种群相互依存模型
• 传染病模型
• …

6.1 传染病模型

第一讲 微分方程模型a

6.1.1 模型一

先以一个简单的模型示例：

存在问题：当 $t\to \infty$ 时，$i\to \infty$ 但很显然地球人口有限不可能是 $\infty$，问题的根源在于该模型未区分人，应区分已感染者（病人）和未感染者（健康人）。

6.1.2 模型二

6.1.3 模型三

传染病无免疫性——病人治愈为健康人，健康人可能被再次感染