python算积分蒙特卡罗_菜鸟之路——数学建模之蒙特卡罗积分（投点法，平均值法）+牛顿法解方程组MATLAB实现...

蒙特卡罗方法这里不再赘述

1，例题

Matlab代码：

%蒙特卡罗法求积分

N=1000; %随机选取1000个点

x=rand(1,N);

y=rand(1,N);

S=sum(y<=sin(x)/x)/N %比较每一个元素，y<=sin(x)/x则为1

运行结果：

S = 0.9150

2，例题

Matlab代码：

N=100000; %随机选取100000个点

x=rand(1,N);

y=rand(1,N);

S=sum(exp((x+y).^2))/length(x) %(1-0)*(1-0)*函数的平均值(也就是长乘宽乘高)

运行结果：

S = 4.9167

3，练习题

1，投点法

N=1000; %随机选取1000个点

x=rand(1,N);

y=rand(1,N);

S=sum(y<=sqrt(1-x.^2))/N

运行结果

S = 0.7940 (半径为1 的1/4圆的面积)

2，平均值法

N=10000; %随机选取10000个点

x=rand(1,N);

y=rand(1,N);

S=sum(sqrt(1-x.^2))/N

运行结果

S = 0.7881 (半径为1 的1/4圆的面积)

Newton法解二元一次方程

Matlab代码

function [ x,k,index ] = Newtons( x , eq,it_max )

index=0;

k=-1;

while(k<=it_max)

x1=x;

[J,F]=funs(x);

x=x-(J\F')';

norm=sqrt((x-x1).^2);

if(norm

index=-1;break

end

k=k+1;

end

end

function [ J,F] = funs(x)

F=[x(1)^2+x(2)^2-5,(x(1)+1)*x(2)-(3*x(1)+1)];

J=[2*x(1),2*x(2);x(2)-3,x(1)+1];

end

[x,k,index]=Newtons([0,1],1e-5,100)

计算结果

x =

1.0000 2.0000

k =

4

index =

-1

知识点

1，函数的建立与嵌套

函数嵌套时直接使用即可，不需要考虑不同的m文件里，需不需要声明一下。

遇到问题：输入参数的数目不足。

就是字面意思，仔细检查一下输入的参数。

2，二元一次方程的求解，也就是求解雅克比方程时。直接Ax=B，令x=A\B,或者x=B*A^(-1)

当时傻了，还想着找个函数解方程。

例题

function [ x,k,index ] = Newtons( x , eq,it_max )

index=0;

k=-1;

while(k<=it_max)

x1=x;

[J,F]=funs(x);

x=x-J\F';

norm=sqrt((x-x1).^2);

if(norm

index=-1;break

end

k=k+1;

end

end

function [J,F] = funs(x)

F=[x(1)^2+x(2)^2-1,x(1)^3-x(2)];

J=[2*x(1),2*x(2);3*x(1)^2,1];

end

[ x,k,index ] = Newtons( [-0.8;0.6] , 10^(-3),10000000 )

运行结果

x =

0.8257

0.5642

k =1563483

index = -1