fleury MATLAB,数学建模部分源码分享

写在前面

数学建模竞赛已经结束了。由于种种原因最终无缘国奖，还是很难受的。结合我的经历，有一些想法想告诉所有正在准备数学建模竞赛的同学。

数学建模竞赛它不会绝对相信实力，运气真的是必须考虑的一部分。

无论数学建模是否拿奖都受益终生，它会极大地提高工科生思维和建模能力。

如果没有必要千万不要熬夜，状态会很差很差很差，得不偿失。

心态爆炸是导致数学建模竞赛失败的大部分原因。

找个好队友真的extremely重要

正题

最近准备换系统，把我们在历次比赛中的代码拿出。反正删了也是删了，还不如传上来以后还能回味一下。如果中间有bug或error欢迎反馈。原谅我不写注释。。。

在此真诚感谢数学建模32种方法的整理者，这里我只提供一个下载链接

层次分析法AHP

clear, clc;

prin = xlsread('AHP.xlsx', 'A1:D4');

price = xlsread('AHP.xlsx', 'A6:C8');

oc = xlsread('AHP.xlsx', 'A10:C12');

conf = xlsread('AHP.xlsx', 'A14:C16');

surf = xlsread('AHP.xlsx', 'A18:C20');

[Vpric, Dpric] = eig(prin);

DpricMax = max(max(Dpric));

CIpric = (DpricMax - length(prin)) ./ (length(prin) - 1);

RI = 0.90;

CRpric = CIpric ./ RI;

if CRpric < 0.10 %#ok

disp(['准则层CR = ', num2str(CRpric), ',接受误差']);

end

wCR = sum(prin, 2)';

indexwCRmax = find(wCR == max(wCR));

switch indexwCRmax

case 1

disp('最大影响因素：price');

case 2

disp('最大影响因素：oil consumption');

case 3

disp('最大影响因素：conformance');

case 4

disp('最大影响因素：surface');

end

wPrice = sum(wCR(1) .* price, 2);

wOc = sum(wCR(2) .* oc, 2);

wConf = sum(wCR(3) .* conf, 2);

wSurf = sum(wCR(4) .* surf, 2);

wScheme = sum([wPrice, wOc, wConf, wSurf], 2);

disp(['最倾向第',num2str(find(wScheme == max(wScheme))) , '辆车']);

单源最短路Dijkstra算法

function [pathMat, pathWeightMat] = Dijkstra(weightmat, startpoint)

%function Dijkstra receives two arguments, the first one is the matrix of path weight,

%the second one is the startpoint.

%function Dijkstra returns two results, pathMat is the matrix of the

%shorest path, pathWeightMat is the matrix of weight

%an example:

% >> m. at = [

% 0 1 5 Inf Inf Inf Inf Inf Inf

% 1 0 3 7 5 Inf Inf Inf Inf

% 5 3 0 Inf 1 7 Inf Inf Inf

% Inf 7 Inf 0 2 Inf 3 Inf Inf

% Inf 5 1 2 0 3 6 9 Inf

% Inf Inf 7 Inf 3 0 Inf 5 Inf

% Inf Inf Inf 3 6 Inf 0 2 7

% Inf Inf Inf Inf 9 5 2 0 4

% Inf Inf Inf Inf Inf Inf 7 4 0]

% >> [path, weight] = Dijkstra(mat, 1)

% path =

% 1 0 0 0 0 0 0 0 0

% 1 2 0 0 0 0 0 0 0

% 1 2 3 0 0 0 0 0 0

% 1 2 3 5 4 0 0 0 0

% 1 2 3 5 0 0 0 0 0

% 1 2 3 5 6 0 0 0 0

% 1 2 3 5 4 7 0 0 0

% 1 2 3 5 4 7 8 0 0

% 1 2 3 5 4 7 8 9 0

weightMat = weightmat;

if ismember(1, weightMat ~= weightMat')

[row, col] = find((weightMat ~= weightMat') == 1);

error(['Matrix input illegal, ', '(', num2str(row(1)), ',', num2str(col(1)),...

')and(', num2str(row(2)) , ',', num2str(col(2)), ')']);

end

unfinishMat = ones(1, length(weightMat(1, :)));

pathWeightMat = inf .* ones(length(weightMat(1, :)), 1);

pathMat = zeros(size(weightMat));

startPoint = startpoint;

presPoint = startPoint;

presPathWeight = 0;

path = presPoint;

pathMat(presPoint, 1) = presPoint;

while ismember(1, unfinishMat)

i = 1;

unfinishMat(presPoint) = 0;

tempPathMat = weightMat(presPoint, :) .* unfinishMat;

tempPathMat(tempPathMat == 0) = inf;

nextPoint = find(tempPathMat == min(tempPathMat), 1, 'last');

if presPoint == nextPoint

break;

end

while i <= length(unfinishMat)

if tempPathMat(i) ~= inf && ~isnan(tempPathMat(i)) && i ~= presPoint

tempPath = [path, i];

tempWeight = presPathWeight + tempPathMat(i);

if tempWeight < pathWeightMat(i)

pathWeightMat(i) = tempWeight;

pathMat(i, :) = [tempPath, zeros(1, length(pathMat(i, :)) - length(tempPath))];

end

i = i + 1;

end

path = pathMat(nextPoint, pathMat(nextPoint, :) ~= 0);

presPathWeight = pathWeightMat(nextPoint);

presPoint = nextPoint;

end

pathWeightMat(startPoint) = 0;

欧拉回路Fleury算法

参考后封装的，原文太久忘了，抱歉

function EulerCircuit = Fleury(G, st)

G(G == inf) = 0;

if nargin == 2

presPoint = st;

circuit = st;

else

presPoint = 1;

circuit = 1;

end

while sum(any(G))

altDeg = (sum(G) - G(presPoint, :)) .* (G(presPoint, :) > 0)

if ~any(altDeg)

[~, nextPoint] = find(G ~= 0)

circuit = [circuit, nextPoint(1)]

break

end

nextPoint = find(altDeg == max(altDeg), 1)

circuit = [circuit, nextPoint]

G(presPoint, nextPoint) = G(presPoint, nextPoint) - 1

G(nextPoint, presPoint) = G(presPoint, nextPoint)

presPoint = nextPoint

altDeg = (sum(G, 2) - G(:, presPoint)) .* (G(:, presPoint) > 0)

if ~any(altDeg)

[nextPoint, ~] = find(G ~= 0)

circuit = [circuit, nextPoint(1)]

break

end

nextPoint = find(altDeg == max(altDeg), 1)

circuit = [circuit, nextPoint]

G(presPoint, nextPoint) = G(presPoint, nextPoint) - 1

G(nextPoint, presPoint) = G(presPoint, nextPoint)

presPoint = nextPoint

end

EulerCircuit = circuit;

完全最短路Flory算法

function [pathMat, pathWeightMat] = Floyd(weightMat, startpoint, endpoint)

% Floyd find out the shorest path between two given vertex

% [a, b] = Floyd(G, s, e)

% G: the graph

% s: start point

% e: end point

% a: the shorest path between s and e

% b: the path length of a

% Example:

% >> G = [ 0 12 19 28 40 59

% 12 0 13 20 29 41

% 19 13 0 14 21 30

% 28 20 14 0 15 22

% 40 29 21 15 0 16

% 59 41 30 22 16 0];

% s = 1; e = 6;

% >> [a, b] = Floyd(G, s, e)

% a =

% 1 3 6

% b =

% 49

if ismember(1, weightMat ~= weightMat')

[row, col] = find((weightMat ~= weightMat') == 1);

error(['Matrix input illegal, ', '(', num2str(row(1)), ',', num2str(col(1)),...

')and(', num2str(row(2)) , ',', num2str(col(2)), ')']);

end

n=length(weightMat);

U=weightMat;

m=1;

%Floyd algorithm

while m<=n

for i=1:n

for j=1:n

if U(i,j)>U(i,m)+U(m,j)

U(i,j)=U(i,m)+U(m,j);

end

m=m+1;

end

pathWeightMat = U(startpoint,endpoint);

P1=zeros(1,n);

k=1;

P1(k)=endpoint;

V=ones(1,n)*inf;

kk=endpoint;

while kk~=startpoint

for i=1:n

V(1,i)=U(startpoint,kk)-weightMat(i,kk);

if V(1,i)==U(startpoint,i)

P1(k+1)=i;

kk=i;

k=k+1;

end

k=1;

wrow=find(P1~=0);

for j=length(wrow):(-1):1

pathMat(k)=P1(wrow(j));

k=k+1;

end

遗传算法(求解哈密顿图的一种优化算法)

注意，此算法可能主要部分不是遗传算法，而是其中的改良圈算法。学校教材代码，无法验证但结果不错，注释部分可删

clc,clear;

close all;

load data;

data = [data; [data(1, 1), data(1, 2)]]; % 数据部分，其中第一行为起点和终点(环)

x = data(:, 1);

y = data(:, 2);

fullData = data;

codeLen = length(x);

libSize = 50;

dist = zeros(codeLen);

for i = 1: codeLen - 1

for j = i + 1: codeLen

pathLen = sqrt((x(i) - x(j)) ^ 2 + (y(i) - y(j)) ^ 2);

dist(i,j) = pathLen;

end

dist = dist + dist';

%通过改良圈算法选取优良父代A

% rec = [];

for k = 1: libSize

c = randperm(codeLen - 2);

c1 = [1, c + 1, codeLen];

flag = 1;

while flag > 0

flag = 0;

for m = 1: codeLen - 3

for n = m + 2: codeLen - 1

if dist(c1(m), c1(n)) + dist(c1(m + 1), c1(n + 1)) < dist(c1(m), c1(m + 1)) + dist(c1(n), c1(n + 1))

flag = 1;

c1(m + 1:n) = c1(n: - 1:m + 1);

end

% temp = 0;

% for i = 1: length(c1) - 1

% temp = temp + dist(c1(i), c1(i + 1));

% end

% rec = [rec, temp];

end

lib(k,c1) = 1: codeLen;

end

[~,index] = sort(lib,2);

pathLen = zeros(length(lib(:, 1)), 1);

for j = 1:length(lib(:, 1))

for i = 1: codeLen - 1

pathLen(j) = pathLen(j) + dist(index(j, i),index(j, i + 1));

end

path = index(pathLen == min(pathLen),:);

len = min(pathLen)

X = fullData(path,1);

Y = fullData(path,2);

figure;

plot(X, Y, ' - o')

% % figure;

% % plot(rec)

数据表data.mat

1304 2312

3639 1315

4177 2244

3712 1399

3488 1535

3326 1556

3238 1229

4196 1004

4312 790

4386 570

3007 1970

2562 1756

2788 1491

2381 1676

1332 695

3715 1678

3918 2179

4061 2370

3780 2212

3676 2578

4029 2838

4263 2931

3429 1908

3507 2367

3394 2643

3439 3201

2935 3240

3140 3550

2545 2357

2778 2826

2370 2975

最小生成树Kruskal算法

function MST = Kruskal(G, method)

%Kruskal Get the MST of graph.

% MST = Kruskal(G, method)

% Kruskal has two arguments. G is necessary while method can be omitted.

% G is the matrix of the graph, in which the elements stands for the path

% weight between the vertexes, and G(i, i) = 0.

% if vertex(i) and vertex(j) are non-adjacent, G(i, j) and G(j, i) could

% be both 0 or inf.

% method points out whether the function will find out the Minimum Spanning Tree

% or the Maximum Spanning Tree. if method = 'min' or omitted, function

% will returns Minimum Spanning Tree, if method = 'max', Maximum Spanning

% Tree will be returned.

% Kruskal returns a Ax3 matrix ,A is the number of edges, column 1 and 2

% are the vertexes of edges, column 3 is the weight of the relevant edge.

% Example:

% if G = [0 5 Inf Inf

% 5 0 2 3

% Inf 2 0 4

% Inf 3 4 0]

% then Kruskal(G) is

% 3 2 2

% 4 2 3

% 2 1 5

% ## [4 2 3] means from vertex4 to vertex2 and the distance is 3.

% Kruskal(G, 'max') is

% 2 1 5

% 4 3 4

% 4 2 3

if nargin == 1 || strcmp(method, 'min')

weightMat = tril(G);

weightMat(weightMat == 0) = inf;

checkMat = zeros(size(weightMat));

path = [];

edgeNum = 0;

while ~all(all(weightMat == inf))

[row, col] = find(weightMat == min(min(weightMat)), 1);

checkMat(row, col) = 1;

checkMat(col, row) = 1;

tempMat = checkMat;

while 1

index = find(sum(tempMat, 2) == 1 | sum(tempMat, 2) == 0, 1);

if isempty(index)

checkMat(row, col) = 0;

break;

end

tempMat(index, :) = [];

tempMat(:, index) = [];

if isempty(tempMat)

path = [path; row, col, weightMat(row, col)];

edgeNum = edgeNum + 1;

break;

end

weightMat(row, col) = inf;

end

elseif nargin == 2 || strcmp(method, 'max')

weightMat = tril(G);

weightMat(weightMat == inf) = -inf;

weightMat(weightMat == 0) = -inf;

checkMat = zeros(size(weightMat));

path = [];

edgeNum = 0;

while ~all(all(weightMat == -inf))

[row, col] = find(weightMat == max(max(weightMat)), 1);

checkMat(row, col) = 1;

checkMat(col, row) = 1;

tempMat = checkMat;

while 1

index = find(sum(tempMat, 2) == 1 | sum(tempMat, 2) == 0, 1);

if isempty(index)

checkMat(row, col) = 0;

break;

end

tempMat(index, :) = [];

tempMat(:, index) = [];

if isempty(tempMat)

path = [path; row, col, weightMat(row, col)];

edgeNum = edgeNum + 1;

break;

end

weightMat(row, col) = -inf;

end

MST = path;

最大流问题

这个代码没有验证过，因为这是培训的时候讲的，回来就顺便仿了

function f=maxFlow(u)

f = zeros(size(u));

n = length(u);

maxf(n) = 1;

while maxf(n)>0

maxf = zeros(1,n);

pred = zeros(1,n);

list = 1;

record = list;

maxf(1) = inf;

%list是未检查邻接点的标号点，record是已标号点

while (~isempty(list))&&(maxf(n) == 0)

flag = list(1);

list(1) = [];

label1 = find(u(flag, :) - f(flag, :));

label1 = setdiff(label1, record);

list = union(list, label1);

pred(label1) = flag;

maxf(label1) = min(maxf(flag), u(flag,label1) - f(flag, label1));

record = union(record, label1);

label2 = find(f(:, flag));

label2 = label2';

label2 = setdiff(label2, record);

list = union(list, label2);

pred(label2) = - flag;

maxf(label2) = min(maxf(flag), f(label2, flag));

record = union(record, label2);

end

if maxf(n) > 0

v2 = n;

v1 = pred(v2);

while v2 ~= 1

if v1 > 0

f(v1,v2) = f(v1, v2) + maxf(n);

else

v1 = abs(v1);

f(v2, v1) = f(v2, v1) - maxf(n);

end

v2 = v1;

v1 = pred(v2);

end

主成分分析PCA

这个算法好像matlab有自带的函数，自行help pca。这个代码只作为原理步骤展示。原文参考地址见代码

%example refer to https://blog.csdn.net/zhongkelee/article/details/44064401

clear, clc;

score = xlsread('PCA.xlsx', 'B2:G10');

%第一步，分别求每列平均值，然后对于所有的样例，都减去对应的均值。

for i = 1:length(score(1,:))

score(:, i) = score(:, i) - mean(score(:, i));

end

%第二步，求特征协方差矩阵

covMatrix = zeros(length(score(1,:)), length(score(1,:)));

for i = 1:length(score(1,:))

for j = 1: length(score(1,:))

temp = cov(score(:, i), score(:, j));

covMatrix(i ,j) = temp(2, 1);

end

%第三步，求协方差的特征值和特征向量

[VcovMatrix, DcovMatrix] = eig(covMatrix);

%第四步，将特征值按照从大到小的顺序排序，选择其中最大的k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵

projectDimension = 3; %提取的主成分个数

projectMatrix = zeros(6, projectDimension);

for k = 1:projectDimension

index = find(sum(DcovMatrix == max(max(DcovMatrix))) == 1);

projectMatrix(:, k) = VcovMatrix(:, index);

DcovMatrix(DcovMatrix == max(max(DcovMatrix))) = 0;

end

%第五步，做乘积

dataProjected = score * projectMatrix

排队论/生灭过程

这个代码好像只适合生灭过程，其他类型自行查询

t_arr = zeros(1,21); %每位顾客到达时间

t_leave = zeros(1,21); %每位顾客离去时间

time_wait_c = zeros(1,21); %顾客等待时间累加

time_free_s = zeros(1,21); %收款台空闲时间累加

t_exp = exprnd(2,1,21); %服从指数分布的随机数

t_norm = normrnd(1,1/3,1,21); %服从正态分布的随机数

for i = 1:1:20

t_arr(i+1) = t_arr(i)+t_exp(i);

if t_arr(i+1) >= t_leave(i)

time_wait_c(i+1) = time_wait_c(i);

t_leave(i+1) = t_arr(i+1)+t_norm(i+1);

time_free_s(i+1) = t_arr(i+1)-t_leave(i)+time_free_s(i);

else

time_wait_c(i+1) = t_leave(i)-t_arr(i+1)+time_wait_c(i);

t_leave(i+1) = t_leave(i)+t_norm(i+1);

time_free_s(i+1) = time_free_s(i);

end

b = [t_arr',t_leave',time_wait_c',time_free_s']; %队伍情况

[brow,bcol] = size(b);

b = [b,zeros(brow,1)];

for j = 2:brow

l = 0;

if j-l-1 > 0

while b(j,1) <= b(j-l-1,2) l = l+1;

end;

b(j,bcol+1) = l;

end;

ave_wait_time = time_wait_c(end)/length(time_wait_c) %平均等待时间

ave_stay_time = sum(t_leave-t_arr)/length(t_arr) %平均逗留时间

served_per_min = t_leave(end)/length(t_leave) %平均每分钟服务的顾客人数

画六面体的函数

这个函数完全是闲得蛋疼写的，当初做了个堆箱子的题的可视化，没别的办法了就只好上matlab，比较吃内存。所以能找到别的方法尽量还是别用这个吧

function simplot(arg)

x0 = arg(1);

y0 = arg(2);

z0 = arg(3);

x1 = arg(4);

y1 = arg(5);

z1 = arg(6);

d = 1;

color = {'blue', 'red', 'green', 'black', 'yellow'};

color = color{randi(5)};

x = x0: d: x1; % x起始坐标: 间隔: 终止坐标

y = y0: d: y1; % y起始坐标: 间隔: 终止坐标

z = z0: d: z1; % z起始坐标: 间隔: 终止坐标

hold on;

% xoy

[X, Y] = meshgrid(x, y);

Z = z0 * ones(size(X));