蒙特卡罗法算法的应用

一、蒙特卡罗法的概念及应用

二、蒙特卡罗法的实质

三、求非线性方程一个实根的蒙特卡罗法

四、求实函数或复函数方程一个复根的蒙特卡罗法

蒙特卡罗法的概念及应用


蒙特卡洛法 (又称统计试验法)是描述装备运用过程中各种随机现象的基本方法，而且它特别适用于一些解析法难以求解甚至不可能求解的问题，因而在装备效能评估中具有重要地位。 用蒙特卡洛法来描述装备运用过程是1950年美国人约翰逊首先提出的。

这种方法能充分体现随机因素对装备运用过程的影响和作用。 更确切地反映运用活动的动态过程。 在装备 效能评估 中，常用蒙特卡洛法来确定含有随机因素的效率指标，如发现概率、 命中概率 、平均毁伤目标数等；模拟 随机服务系统 中的随机现象并计算其数字特征；对一些复杂的装备运用行动，通过合理的分解，将其简化成一系列前后相连的事件，再对每一事件用随机抽样方法进行模拟，最后达到模拟装备运用活动或运用过程的目的。

  **蒙特卡罗法的实质**
  蒙特卡罗法的实质：当无法求得精确解时，进行随机抽样，根据统计试验求近似解。

通俗理解：

1、假如有一万个苹果，需要挑选出其中最大的，但只能闭着眼睛挑，手里最多保留一个苹果。

2、初始时闭着眼睛随机挑选，拿到一个苹果。

3、然后可继续闭着眼睛挑一个，并与手里现有的比较，留下较大的、扔掉较小的。

4、循环重复上一步，则挑的次数越多，挑出最大苹果的可能性也就越大

5、但人的时间精力有限，除非把十万苹果都挑一遍，否则无法确定挑出来的就是最大的。

6、所以挑了3000后，就把此时手里的苹果视为十万个苹果中最大的近似解。

求非线性方程一个实根的蒙特卡罗法

算法描述：
① 令i = 0 i=0i=0，任选一个初值x ( 0 ) x^{(0)}x 
(0)，计算F ( 0 ) = f ( x ( 0 ) ) F^{(0)}=f(x^{(0)})F 
(0) =f(x (0) )
②产生区间上均匀分布的随机数a aa，以x ( 0 ) + a x^{(0)}+ax 
(0)+a作为自变量并计算F ( 1 ) = f ( x ( 1 ) + a ) F^{(1)}=f(x^{(1)}+a)F (1) =f(x (1) +a)
③若F ( 0 ) ⩽ F ( 1 ) F^{(0)} \leqslant F^{(1)}F 
(0)⩽F (1) 则返回上一步并重新产生随机数，直到找到某个a aa使得F ( 1 ) ⩽ F ( 0 ) F^{(1)} \leqslant F^{(0)}F 
(1) ⩽F (0)
 ④ 若F ( 1 ) {F^{(1)}}F (1)的绝对值⩽ σ \leqslant \sigma⩽σ，则x ( 1 ) + a x^{(1)}+ax (1) +a作为方程的一个根，否则转至②

function root = Mont(f,B,x0,eps)
% f：方程表达式
% B：随机数区间
% x0：初始值
% eps：根的精度
format long
if(nargin == 3)
	eps = 1.0e-4;
end
Fx = subs(sym(f),findsym(sym(f)),x0);
while abs(Fx)>eps
	Fx = subs(sym(f),findsym(sym(f)),x0);
	Fx1 = abs(Fx) + 10;
	while(abs(Fx1) >= abs(Fx))
		a = 2*B*rand() - B;
		x1 = x0 + a;
		Fx1 = subs(sym(f),findsym(sym(f)),x1);
	end
	x0 = x1;
end
root = x1;

function x=qiushigeng(fname,a,b,e)
%fname 为函数名，a,b 为区间端点，e 为精度
fa=feval(fname,a); %把 a 端点代入函数，求 fa
fb=feval(fname,b); %把 b 端点代入函数，求 fb
if fa*fb<0 error('两端函数值为同号');
end
%如果 fa*fb<0，则输出两端函数值为同号
k=0
x=(a+b)/2
while(b-a)>(2*e) %循环条件的限制
fx=feval(fname,x);%把 x 代入代入函数，求 fx
if fa*fx>0%如果 fa 与 fx 同号，则把 x 赋给 b，把 fx 赋给 fb
b=x;
fb=fx;
else
%如果 fa 与 fx 异号，则把 x 赋给 a，把 fx 赋给 fa
a=x;
fa=fx;
end
k=k+1
x=(a+b)/2 
end

在命令行窗口输入以下代码，即可求出一个非线性方程的一个实根。

fun=inline('exp(-x*x*x)-sin(x)/cos(x)+800');x=qiushigen(fun,-0.5,0.5,0.5.*10.^-3)

运行结果如下：




求实函数或复函数方程一个复根的蒙特卡罗法


定义两个函数

function root = Mont(f,B,x0,eps)
% f：方程表达式
% B：随机数区间
% x0：初始值
% eps：根的精度
format long
if(nargin == 3)
	eps = 1.0e-4;
end
Fx = subs(sym(f),findsym(sym(f)),x0);
while abs(Fx)>eps
	Fx = subs(sym(f),findsym(sym(f)),x0);
	Fx1 = abs(Fx) + 10;
	while(abs(Fx1) >= abs(Fx))
		a = 2*B*rand() - B;
		x1 = x0 + a;
		Fx1 = subs(sym(f),findsym(sym(f)),x1);
	end
	x0 = x1;
end
root = x1;

%求复数的正整数次方根，并且返回所有复数根。
% X可以是标量，也可以是向量。
% n必须是正整数。
% Y是L*n的矩阵，其中L是X的长度。
function Y=complex_sqrt(X,n)
[~,b]=size(X);
if b>1
    X=X.';%进行转置但复数元素的值不会改变，使用"'"转置会造成复数元素的值改变
end
pre_angle=angle(X)+2*pi*(angle(X)<0);%angle的输出范围是[-pi,pi]，我们希望的范围是[0,2pi)
complex_angle=(pre_angle+(0:2*pi:2*(n-1)*pi))/n;
Y=repmat(abs(X).^(1/n),[1,n]).*(cos(complex_angle)+1i*sin(complex_angle));

在命令行窗口输入以下代码：

 X=[4*sqrt(2)*(-1+1i),-i]

>>complex_sqrt(X,3)

在次声明：本次实验是为了完成测试，不具有很大的参考价值。

相关链接：



(https://blog.csdn.net/weixin_55105724/article/details/122767399)



相关链接：



(https://blog.csdn.net/qq_67421354/article/details/127683589)