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

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 模型三

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