一.目录

二.引例1：布丰投针实验

1.投针步骤

步骤：
1：取一张白纸，在上面画出许多条间距为a的平行线
2：去一根长度为l的针，随机的向纸张上投掷n次，观察针与直线相交的次数记为m
3：计算针与直线相交的概率

2.思路

3. 历史实验者

二.蒙特卡罗模拟概述

1.定义

蒙特卡罗方法又称统计模拟法，是一种随机的模拟方法，以概率论和数理统计为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。

2.提出

于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”中首先提出。1977年Buffon提出用投针实验求圆周率。被认为是起源。

3.原理

有大数定理而知，当样本的容量足够大时，事件的发生频率就是时间的概率。

4.讨论

讨论
（1）：蒙特卡洛只能是一种方法，并不能算是一种算法。如果所求解的问题与概率模型有一定的关联，那么我们就可以使用计算机多次模拟事件的发生，以获得问题的近似解。
（2）：计算机仿真是一门用随机数解决随即问题的方法，其主要是用在复杂问题的数值模拟上。在数学建模上就是使用计算机对问题进行模拟。
（3）：枚举法是把所有的有可能发生的情况都考虑进去，最终算出的一个确定的结果。与蒙特卡洛显然易见不同。蒙特卡洛是模拟的次数越多，理论上计算的也就越精准。可以看成蒙特卡洛方法是枚举法的一种变异。

三.应用实例

1.三门问题

【1】问题

你参加⼀档电视节⽬，节⽬组提供了ABC三扇⻔，主持⼈告诉你，其中⼀扇⻔后边有辆汽⻋，其它两扇⻔后是空的。
假如你选择了B⻔，这时，主持⼈打开了C⻔，让你看到C⻔后什么都没有，然后问你要不要改选A⻔？

【2】代码

预备知识

%% （1）预备知识
%  rand(m,n)函数产生由在[0,1]之间均匀分布的随机数组成的m行n列的矩阵（或称为数组）。
rand(5,4)
%     0.8300    0.1048    0.2396    0.4398
%     0.5663    0.1196    0.8559    0.5817
%     0.9281    0.2574    0.3013    0.9355
%     0.3910    0.3173    0.2108    0.1676
%     0.3645    0.4372    0.8819    0.9232
rand(3) % 若只给一个输入，则会生成一个方阵
%     0.1709    0.4951    0.0541
%     0.9374    0.8500    0.6155
%     0.2400    0.3156    0.5741
% a + rand(m,n)*(b-a) 可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。
-2 + rand(3,2) * (2 - (-2))  % 输出在[-2,2]之间均匀分布的随机数组成的3行2列的矩阵。
%    -1.2743    0.6013
%    -1.3084    0.0766
%     1.5075    0.7563
% a + rand(m,n)*(b-a)等价于unifrnd(a,b,m,n)
unifrnd(-2,2,3,2)

示例代码

%% （2）代码部分
l =  0.520;     % 针的长度（任意给的）
a = 1.314;    % 平行线的宽度(大于针的长度l即可)
n = 1000000;    % 做n次投针试验，n越大求出来的pi越准确
m = 0;    % 记录针与平行线相交的次数
x = rand(1, n) * a / 2 ;   % 在[0, a/2]内服从均匀分布随机产生n个数， x中每一个元素表示针的中点和最近的一条平行线的距离
phi = rand(1, n) * pi;    % 在[0, pi]内服从均匀分布随机产生n个数，phi中的每一个元素表示针和最近的一条平行线的夹角
% axis([0,pi, 0,a/2]);   box on;  % 画一个坐标轴的框架，x轴位于0-pi，y轴位于0-a/2， 并打开图形的边框
for i=1:n  % 开始循环，依次看每根针是否和直线相交if x(i) <= l / 2 * sin(phi (i))     % 如果针和平行线相交m = m + 1;    % 那么m就要加1
%         plot(phi(i), x(i), 'r.')   % 模仿书上的那个图，横坐标为phi，纵坐标为x , 用红色的小点进行标记
%         hold on  % 在原来的图形上继续绘制end
end
p = m / n;    % 针和平行线相交出现的频率
mypi = (2 * l) / (a * p);  % 我们根据公式计算得到的pi
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])%% （3） 由于一次模拟的结果具有偶然性，因此我们可以重复100次后再来求一个平均的pi
result = zeros(100,1);  % 初始化保存100次结果的矩阵
l =  0.520;     a = 1.314;
n = 1000000;
for num = 1:100m = 0;  x = rand(1, n) * a / 2 ;phi = rand(1, n) * pi;for i=1:nif x(i) <= l / 2 * sin(phi (i))m = m + 1;endendp = m / n;mypi = (2 * l) / (a * p);result(num) = mypi;  % 把求出来的myphi保存到结果矩阵中
end
mymeanpi = mean(result);  % 计算result矩阵中保存的100次结果的均值
disp(['蒙特卡罗方法得到pi为：', num2str(mymeanpi)])

2.模拟排队问题

【1】题目

假设某银行工作时间只有⼀个服务窗口，工作人员只能逐个的接待顾客。当来的顾客较多时，⼀部分顾客就需要排队等待。
假设：
1) 顾客到来的间隔时间服从参数为0.1的指数分布
2) 每个顾客的服务时间服从均值为10，方差为4的正态分布(单位为分钟，若服务时间小于1分钟，则按1分钟计算)
3) 排队按先到先服务的规则，且不限制队伍的长度，每天工作时长为8小时。

试回答下⾯的问题：1) 模拟⼀个工作日，在这个工作日共接待了多少客户，客户平均等待的时间为多少? 2) 模拟100个工作日，计算出平均每日接待客户的个数以及每日客户的平均等待时长

【2】引入参数

CCC_i：第iii个客户到达时间
bbb_i：第iii个客户开始服务时间
eee_i：第iii个客户服务结束时间

【3】总结题目信息

1. 第i−1i-1i−1个客户和第iii个客户到达的间隔时间xxx_i~E(0,1)E(0,1)E(0,1),参数lamdba=0.1lamdba=0.1lamdba=0.1的指数分布（均值为10）
   
2. 第iii个客户的服务持续时间yyy_{i$N(10,4)$,若$y$}i_{$<1$=》令$y$}i~=1=1=1

【4】过程可视化

【5】一些推导公式

【6】代码

1.预备代码

%% （1）预备知识
% normrnd(MU,SIGMA):生成一个服从正态分布(MU参数代表均值,SIGMA参数代表标准差,方差开根号是标准差)的随机数
normrnd(10,2)  % 均值为10 标准差为2（方差为4）的正态分布随机数
% exprnd(M)表示生成一个均值为M的指数分布随机数(其对应的参数为1/M)
exprnd(5)  % 均值为5的指数分布随机数（对应的参数为0.2）
% mean函数是用来求解均值的函数（第一期视频第五讲）
mean([1,2,3])
% tic函数和toc函数可以用来返回代码运行的时间，例如我们要计算一段代码的运行时间，就可以在这段代码前加上tic，在这段代码后加上toc (我的微信公众号"数学建模学习交流"中有一篇推送《为什么要对代码初始化》中使用过这对函数)
tic
a = 2^100
toc

2.示例代码

%% （2）模型中出现的变量的说明
% x(i)表示第i-1个客户和第i个客户到达的间隔时间，服从参数为0.1的指数分布
% y(i)表示第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布 (若小于1则按1计算)
% c(i)表示第i个客户的到达时间，那么c(i) = c(i-1) + x(i)，初始值c0=0
% b(i)表示第i个客户开始服务的时间
% e(i)表示第i个客户结束服务的时间，初始值e0=0
% 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
% 即：e(i) = b(i) + y(i）
% 第i个客户开始服务的时间取决于该客户的到达时间和上一个客户结束服务的时间
% 即：b(i) = max(c(i),e(i-1))，初始值b1=c1;
% 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
% 即：wait(i) = b(i) - c(i)
% w表示所有客户等待时间的总和
% 假设一天内银行最终服务了n个顾客，那么客户的平均等待时间t = w/n%% （3）问题1的代码
clear
tic  %计算tic和toc中间部分的代码的运行时间
i = 1;  % i表示第i个客户，最开始取i=1
w = 0;  % w用来表示所有客户等待的总时间，初始化为0
e0 = 0;  c0 = 0;   % 初始化e0和c0为0
x(1) = exprnd(10);  % 第0个客户(假想的)和第1个客户到达的时间间隔
c(1) = c0 + x(1);  % 第1个客户到达的时间
b(1) = c(1); % 第1个客户的开始服务的时间
while b(i) <= 480  % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布if y(i) < 1  % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算y(i) = 1;ende(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间w = w + wait(i); % 更新所有客户等待的总时间i = i + 1; % 增加一名新的客户x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间
end
n = i-1; % n表示银行一天8小时一共服务的客户人数
t = w/n; % 客户的平均等待时间
disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])
disp(['客户的平均等待时间为: ',num2str(t)])
toc  %计算tic和toc中间部分的代码的运行时间%% （4）问题2的代码
clear
tic  %计算tic和toc中间部分的代码的运行时间
day = 100;  % 假设模拟100天
n = zeros(day,1); % 初始化用来保存每日接待客户数结果的矩阵
t = zeros(day,1); % 初始化用来保存每日客户平均等待时长的矩阵
for k = 1:dayi = 1;  % i表示第i个客户，最开始取i=1w = 0;  % w用来表示所有客户等待的总时间，初始化为0e0 = 0;  c0 = 0;   % 初始化e0和c0为0x(1) = exprnd(10);  % 第0个客户(假想的)和第1个客户到达的时间间隔c(1) = c0 + x(1);  % 第1个客户到达的时间b(1) = c(1); % 第1个客户的开始服务的时间while b(i) <= 480  % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布if y(i) < 1  % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算y(i) = 1;ende(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间w = w + wait(i); % 更新所有客户等待的总时间i = i + 1; % 增加一名新的客户x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间endn(k) = i-1; % n(k)表示银行第k天服务的客户人数t(k) = w/n(k); % t(k)表示该银行第k天客户的平均等待时间
end
disp([num2str(day),'个工作日中，银行每日平均服务的客户人数为: ',num2str(mean(n))])
disp([num2str(day),'个工作日中，银行每日客户的平均等待时间为: ',num2str(mean(t))])
toc  %计算tic和toc中间部分的代码的运行时间

3.有约束的非线性规划问题

【1】问题的分类

【2】例子

【3】推导过程

【4】代码

1.预备知识

%% （1）预备知识
%  (1) format long g  可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
5/7
5895*514100
format long g
5/7
5895*514100
%  (2)unifrnd(a,b,m,n)可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。(等价于 a + rand(m,n)*(b-a))
unifrnd(0,5,4,3)
%           4.07361843196589          3.16179623112705          4.78753417717149
%            4.5289596853781         0.487702024997048          4.82444267599638
%           0.63493408146753          1.39249109433524         0.788065408387741
%            4.5668792806951          2.73440759602492          4.85296390880308

2.示例代码

%%  蒙特卡罗求解有约束的非线性规划问题
% max f(x) = x1*x2*x3
% s.t.
% (1) -x1+2*x2+2*x3>=0
% (2) x1+2*x2+2*x3<=72
% (3) x2<=20 & x2>=10
% (4) x1-x2 == 10

%% （2）代码部分
clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(20,30,n,1);  % 生成在[20,30]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(-10,16,n,1);  % 生成在[-10,16]之间均匀分布的随机数组成的n行1列的向量构成x3
fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）
for i=1:nx = [x1(i), x2(i), x3(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2, x3]if (-x(1)+2*x(2)+2*x(3)>=0)  &  (x(1)+2*x(2)+2*x(3)<=72)     % 判断是否满足条件result = x(1)*x(2)*x(3);  % 如果满足条件就计算函数值if  result  > fmax  % 如果这个函数值大于我们之前计算出来的最大值fmax = result;  % 那么就更新这个函数值为新的最大值X = x;  % 并且将此时的x1 x2 x3保存到一个变量中endend
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间%% （3）缩小范围重新模拟得到更加精确的取值
clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(22,23,n,1);  % 生成在[22,23]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(11,13,n,1);  % 生成在[11,13]之间均匀分布的随机数组成的n行1列的向量构成x3
fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）
for i=1:nx = [x1(i), x2(i), x3(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2, x3]if (-x(1)+2*x(2)+2*x(3)>=0)  &  (x(1)+2*x(2)+2*x(3)<=72)     % 判断是否满足条件result = x(1)*x(2)*x(3);  % 如果满足条件就计算函数值if  result  > fmax  % 如果这个函数值大于我们之前计算出来的最大值fmax = result;  % 那么就更新这个函数值为新的最大值X = x;  % 并且将此时的x1 x2 x3保存到一个变量中endend
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

4.书店买书问题（0-1规划）

【1】题目

【2】代码

预备知识

%% （1）预备知识
% (1)unique函数: 剔除一个矩阵或者向量的重复值，并将结果按照从小到大的顺序排列
% adj.  唯一的; 独一无二的   [ju'ni:k]
unique([1 2 5; 6 8 9;2 4 6])
unique([5 6 8 8 4 1 6 2 2 4 8 4 5 6])% (2)randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
randi([-5,5],2,6)

%% （2）代码求解
min_money = +Inf;  % 初始化最小的花费为无穷大，后续只要找到比它小的就更新
min_result = randi([1, 6],1,5);  % 初始化五本书都在哪一家书店购买，后续我们不断对其更新
%若min_result = [5 3 6 2 3]，则解释为：第1本书在第5家店买，第2本书在第3家店买，第3本书在第6家店买，第4本书在第2家店买，第5本书在第3家店买
n = 100000;  % 蒙特卡罗模拟的次数
M = [18     39 29  48  5924    45  23  54  4422    45  23  53  5328    47  17  57  4724    42  24  47  5927    48  20  55  53];  % m_ij  第j本书在第i家店的售价
freight = [10 15 15 10 10 15];  % 第i家店的运费
for k = 1:n  % 开始循环result = randi([1, 6],1,5); % 在1-6这些整数中随机抽取一个1*5的向量，表示这五本书分别在哪家书店购买index = unique(result);  % 在哪些商店购买了商品，因为我们等下要计算运费money = sum(freight(index)); % 计算买书花费的运费% 计算总花费：刚刚计算出来的运费 + 五本书的售价for i = 1:5   money = money + M(result(i),i);  endif money < min_money  % 判断刚刚随机生成的这组数据的花费是否小于最小花费，如果小于的话min_money = money  % 我们更新最小的花费min_result = result % 用这组数据更新最小花费的结果end
end
min_money   % 18+39+48+17+47+20
min_result

5.导弹追踪问题

【1】题目

【2】解题

【3】代码

预备代码

%% （1）预备知识
% mod(m,n)表示求m/n的余数
mod(8,3)
mod(1000,50)% 设置横纵坐标的范围并标上字符
x = 1:0.01:3;
y = x .^ 2;
plot(x,y)  % 画出x和y的图形
axis([0 3 0 10])  % 设置横坐标范围为[0, 3] 纵坐标范围为[0, 10]
pause(3)  % 暂停3秒后再继续接下来的命令
text(2,4,'清风')  % 在坐标为(2,4)的点上标上字符串：清风
close % 关闭图形窗口

%% (2) 代码求解
% 1. 不画追击的示意图
clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
while(dd>=0.001)  % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）t=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离tan_alpha=(y(2)-y(1))/(x(2)-x(1));   % 计算斜率，即tan(α)cos_alpha=sqrt(1/(1+tan_alpha^2));   % sec(α)^2 = (1+tan(α)^2)sin_alpha=sqrt(1-cos_alpha^2);  % sin(α)^2 +cos(α)^2 = 1x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置if d>50  % 导弹的有效射程为50个单位disp('导弹没有击中B船');break;  % 退出循环endif d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）disp(['导弹飞行',num2str(d),'单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])end
end% 2. 画追击的示意图
clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
for i=1:2plot(x(i),y(i),'.k','MarkerSize',1);  % 画出导弹和B船所在的坐标，点的大小为1，颜色为黑色(k)，用小点表示grid on;  % 打开网格线hold on;  % 不关闭图形，继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量  为了控制画图的速度（因为Matlab中画图的速度超级慢）
while(dd>=0.001)  % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）t=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离tan_alpha=(y(2)-y(1))/(x(2)-x(1));   % 计算斜率，即tan(α)cos_alpha=sqrt(1/(1+tan_alpha^2));   % 利用公式：sec(α)^2 = (1+tan(α)^2)  计算出cos(α)sin_alpha=sqrt(1-cos_alpha^2);  % 利用公式： sin(α)^2 +cos(α)^2 = 1  计算出sin(α)x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha;   % 计算新的导弹的位置k = k +1 ;  if mod(k,500) == 0   % 每刷新500次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数for i=1:2plot(x(i),y(i),'.k','MarkerSize',1);hold on; % 不关闭图形，继续画图endpause(0.001);  % 暂停0.001s后再继续下面的操作endif d>50  % 导弹的有效射程为50个单位disp('导弹没有击中B船');break;  % 退出循环endif d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）disp(['导弹飞行',num2str(d),'个单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])end
end

6.旅行商问题

【1】题目

【2】代码

预备代码

%% （1）预备知识
plot([1,2],[5,10],'-o') % 画出一条线段，x范围是[1, 2] ，y范围是[5,10]
text(1.5,7.5,'清风') % 在坐标(1.5,7.5)处标上文本：清风
close% randperm函数的用法
randperm(5)  % 生成1-5组成的一个随机序列(类似于洗牌的操作)
%      3     5     1     2     4
%      1     4     5     3     2

%% （2）代码求解
clear;clc
% 只有10个城市的简单情况coord =[0.6683 0.6195 0.4    0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488 ;0.2536 0.2634 0.4439 0.1463 0.2293 0.761  0.9414 0.6536 0.5219 0.3609]' ;  % 城市坐标矩阵，n行2列
% 38个城市，TSP数据集网站(http://www.tsp.gatech.edu/world/djtour.html) 上公测的最优结果6656。% coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;11600.000000,43150.000000;11690.555600,42686.666700;11715.833300,41836.111100;11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;12421.666700,42895.555600;12645.000000,42973.333300];n = size(coord,1);  % 城市的数目figure(1)  % 新建一个编号为1的图形窗口
plot(coord(:,1),coord(:,2),'o');   % 画出城市的分布散点图
for i = 1:ntext(coord(i,1)+0.01,coord(i,2)+0.01,num2str(i))   % 在图上标上城市的编号（加上0.01表示把文字的标记往右上方偏移一点）
end
hold on % 等一下要接着在这个图形上画图的d = zeros(n);   % 初始化两个城市的距离矩阵全为0
for i = 2:n  for j = 1:i  coord_i = coord(i,:);   x_i = coord_i(1);     y_i = coord_i(2);  % 城市i的横坐标为x_i，纵坐标为y_icoord_j = coord(j,:);   x_j = coord_j(1);     y_j = coord_j(2);  % 城市j的横坐标为x_j，纵坐标为y_jd(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2);   % 计算城市i和j的距离end
end
d = d+d';   % 生成距离矩阵的对称的一面min_result = +inf;  % 假设最短的距离为min_result，初始化为无穷大，后面只要找到比它小的就对其更新
min_path = [1:n];   % 初始化最短的路径就是1-2-3-...-n
N = 10000000;  % 蒙特卡罗模拟的次数
for i = 1:N  % 开始循环result = 0;  % 初始化走过的路程为0path = randperm(n);  % 生成一个1-n的随机打乱的序列for i = 1:n-1  result = d(path(i),path(i+1)) + result;  % 按照这个序列不断的更新走过的路程这个值endresult = d(path(1),path(n)) + result;  % 别忘了加上从最后一个城市返回到最开始那个城市的距离if result < min_result  % 判断这次模拟走过的距离是否小于最短的距离，如果小于就更新最短距离和最短的路径min_path = path;min_result = resultend
end
min_path
min_path = [min_path,min_path(1)];   % 在最短路径的最后面加上一个元素，即第一个点（我们要生成一个封闭的图形）
n = n+1;  % 城市的个数加一个（紧随着上一步）
for i = 1:n-1 j = i+1;coord_i = coord(min_path(i),:);   x_i = coord_i(1);     y_i = coord_i(2); coord_j = coord(min_path(j),:);   x_j = coord_j(1);     y_j = coord_j(2);plot([x_i,x_j],[y_i,y_j],'-')    % 每两个点就作出一条线段，直到所有的城市都走完pause(0.5)  % 暂停0.5s再画下一条线段hold on
end

数学建模更新10（蒙特卡罗模拟）相关推荐

1. Python小白的数学建模课-10.微分方程边值问题

   小白往往听到微分方程就觉得害怕,其实数学建模中的微分方程模型不仅没那么复杂,而且很容易写出高水平的数模论文. 本文介绍微分方程模型边值问题的建模与求解,不涉及算法推导和编程,只探讨如何使用 Pytho ...

2. 数学建模优化和仿真模拟的区别001

   如果我们过于爽快地承认失败,就可能使自己发觉不了我们非常接近于正确!---卡尔·波普尔 数学建模优化和仿真模拟的区别 几十年来,工程师们一直使用模型(Model)来帮助他们理解系统过程并确定最佳解决方 ...

3. 关于投篮的数学建模模型_投篮问题的数学建模[共10页]

   投篮问题的数学建模[共10页] 摘 要如今全民大爱篮球运动,投球的命中率是一场比赛输赢的关键所在,能否投入篮筐与投球时运动员所处的位置.投球时的角度和投球时的出手速度有很大关系,该论文主要以罚球为出发 ...

4. 数学建模更新12（非线性规划）

   非线性规划 一.非线性规划的标准型 1.概述 2.例题 3.MATLAB的非线性规划的函数 4.基础例题 [1]默认方法求例1 [2]用其他方法求例2 [3]改变初始值的影响 [4]使用蒙特卡罗的方法 ...

5. 2021年暑假数学建模第三次模拟赛：全国各省人均消费分析与预测（很好的SPSS统计训练）

   本系列赛题.数据获取: 2021年暑假数学建模模拟赛(赛题+数据+分析) 不直接提供论文等资料,分析已经很详细了 整理不易,欢迎点赞+关注+收藏 赛题 分析 这次的赛题完全可以全用SPSS来实现,以前 ...

6. 【数学建模】基于matlab模拟疫情SEIRS模型【含Matlab源码 2214期】

   一.⛄SEIR模型简介 1 SEIR模型简介 如果所研究的传染病有一定的潜伏期,与病人接触过的健康人并不马上患病,而是成为病原体的携带者,归入 E 类.此时有: 仍有守恒关系 S(t) + E(t) ...

7. Matlab基础编程知识处理(2)(数学建模中模型的模拟与数据提取,本篇全干货)

   今天我终于没有再当鸽子了. 如何提取运算中的数值信息 这里主要是介绍矩阵存储数据的思想. 先假设一个情形: 我在分析模拟电网波动时,想了解究竟有多少次电频率超过我设置的舒适值,那么我该如何统计超过的次 ...

8. 数学建模更新13（MATLAB绘制三维图【上】）

   MATLAB绘制三维图 一.mesh函数以及拓展函数 1.mesh(X,Y,Z)的用法 [1]X是n维向量,Y是m维向量,Z是m*n维的矩阵 [2]X.Y和Z都是m*n维的矩阵 2.mesh(Z)的用 ...

9. 数学建模学习笔记——蒙特卡洛模拟

   n=3  abca acba n=4   3!

10. 数学建模——更新1——excel直方图

    目录 1.插入直方图 2.绘制频率分布直方图 步骤 1. 2. 3.生成表 4.生成频率 5.生成区间字符串 6.画柱状图 1.插入直方图 2.绘制频率分布直方图 步骤 1. 2. 3.生成表 4.生 ...

最新文章

1. 90 后利用平台漏洞薅羊毛，获利 45 万被抓捕！网友们却争论不休……
2. linux 删除乱码文件
3. lua工具库penlight--06数据（一）
4. 【Java】异常处理的目的
5. 华为鸿蒙5g售价,华为首款5G手机售价公布，余承东透露鸿蒙将用于连接家庭设备...
6. 八皇后问题（非递归版）
7. Memcached(四)Memcached的CAS协议
8. 牛客网暑期ACM多校训练营(第三场): E. Sort String（KMP）
9. 大白话告诉你什么是java
10. iOS录音、播放、WAV以及caf转成MP3上传后台
11. 快速搜索Wox工具之Everything Service没有运行报错，解决办法！
12. PSO算法及其对函数优化问题的处理+PSO算法改进
13. python tokenize()_tokenize- 用于Python源代码的Python – Python语言服务（Python教程）（参考资料）...
14. Excel的文件打开特别慢,xls文件特别大解决一例
15. python 提取出图片特定区域的平均rgb值
16. [Python3]pandas.merge用法详解
17. python学习 之 pyqt5前后端分离试验（进度条）
18. 每日一题-Acwing2058笨拙的手指-位运算-枚举
19. Swin Transformer v2实战：使用Swin Transformer v2实现图像分类（一）
20. 小行星大小计算机,一颗有着吉萨大金字塔尺寸的小行星今天掠过地球

热门文章

1. 【论文阅读】VulCNN: An Image-inspired Scalable Vulnerability Detection System
2. 淘宝网或阿里巴巴模拟登陆获取数据
3. 大点再大点个性化超大MSN头像(转)
4. 毕设-基于SSM仓库管理系统
5. 毕向东java ppt下载_毕向东Java基础ppt课件
6. 用python代码制作视频
7. 配置alexa skill（二）
8. PostgreSQL客户端验证
9. Premiere Pro Guru: Speed Changes Premiere Pro 大师教程之改变速度 Lynda课程中文字幕
10. 修改IP、DNS、MAC工具VC源码实现