matlab逆变换法产生随机数_matlab数值积分方法(一)

一元函数定积分的数学表示为

在被积函数f(x)理论上不可积时，即无法求出该积分的解析解，所以往往要采用数值方法来求解。求解定积分的数值方法是多种多样的，如简单的梯形法、Simpson 法、Romberg 法等算法都是数值分析课程中经常介绍的方法。它们的基本思想都是将整个积分空间[a, b]分割成若干个子空间。

,其中

,这样整个积分问题就分解为下面的求和形式。

粗区间划分

细区间划分

% 语法:
% 版本9.6.0.1072779 (R2019a) on PCWIN64
% shi z.y. (155****@qq.com), 2019-09-14 12:36
%-------------------------------------------------------------------------clear variable;close all
clc;format compact
%% -------------------------------------------------------------------------
a= 0;
b = pi;
n =50;
x = linspace(a,b,n);
y = x.^2;
plot(x,y,'r-')
hold on
for i = 1:length(x)-1plot([x(i) x(i) x(i+1)],[0 y(i) y(i)],'b-')hold on
end
title('y=x^2,n=50')

梯形求积

而在每一个小的子空间上都可以近似地求解出来，当然最简单的求每-一个小的子空间的积分方法是采用梯形近似的方法。梯形方法还可以应用于已知数据样本点的数值积分问题求解。假设在实验中测得一组数据(x1,y1), (x2,y2), (x3,y3), ... (xN,yN),且x为严格单调递增的数值，直接求取这些点对应曲线的数值积分最直观的方法就是用梯形方法.

matlab积分实现

%例: 试用梯形法求出x∈(0,π)区间内，函数sin x的定积分值。x = linspace(0,pi,1000);
y= sin(x);
S= sum((2*y(1:end-1)+diff(y)).*(diff(x)))/2

辛普森算法求积

单变量函数的数值积分还可以采用一般数值分析中介绍的其他算法进行求解。例如，可以采用下面给出的Simpson方法求解出

上的积分的近似值为

y = quad(fun,a,b)

其中，Fun为描述被积函数的字符串变量，可以是一个Fun.m函数文件名，该函数的一般格式为y=Fun(x)，还可以用inline()函数和匿名函数直接定义。a, b分别为定积分的上限和下限.

多重积分

使用MATLAB提供的dblquad()函数就可以直接求出上述双重定积分的数值解。该函数的调用格式为

y=db1quad (Fun ,Xm,IM , Ym ,YM)%矩形区域的双重积分
y=dblquad(Fun,xm,TM,Ym,YM,E) % 限定精度的双重积分

注意，本函数不能返回被积函数调用次数，故用户可以自已在被积函数中设置一个计数器，从而测出调用次数。

蒙特卡洛积分

Monte Carlo法是通过大量实验来求取随机变量近似值的一种常用的方法，在现代科学研究中经常用来求解一些建模困难的问题。

考虑图3中给出的示意图。假设正方形的边长为1，可见，四分之一圆的面积是π/4,其面积和正方形面积的比是

，换句话说，如果产生一个均匀分布的随机数，它落到四分之一圆的概率为

。生成N组随机数x和y，使其均为区间[0,1]内均匀分布的随机数。这样记满足

概率为Ni，则对大量的实验数据，有Ni/N≈π/4,.如果N足够大，则可以通过下面的式子近似求出π的值。

% 语法:
% 版本9.6.0.1072779 (R2019a) on PCWIN64
% shi z.y. (155****@qq.com), 2019-09-14 12:51
%-------------------------------------------------------------------------clear variable;close all
clc;format compact
%% -------------------------------------------------------------------------
N=1e5;
x= rand(N,1);
y = rand(N,1);
i = (x.^2+y.^2)<1;
p = sum(i)/N*4
p =3.1387
>>

echo on
% 语法:
% 版本9.6.0.1072779 (R2019a) on PCWIN64
% shi z.y. (155****@qq.com), 2019-09-14 12:51
%-------------------------------------------------------------------------clear variable;close all
clc;format compact
%% -------------------------------------------------------------------------
N=1e3;
x= rand(N,1);
y = rand(N,1);
i = (x.^2+y.^2)<1;
j = (x.^2+y.^2)>=1;
p = sum(i)/N*4
p =3.0440
r =1;
theta=0:0.01:pi/2;
x1 = r*cos(theta);
y1 =r*sin(theta);
plot(x1,y1,'r-','linewidth',2)
axis equal
axis([0 1 0 1])
box on
grid on
hold on
scatter(x(i),y(i),'b')scatter(x(j),y(j),'g')
>>