高文華，吳培浩，崔超宇，張宇風(fēng)，廖嘉祺

(華南理工大學(xué)　a.數(shù)學(xué)學(xué)院;b.機(jī)械與汽車學(xué)院，廣州　510640)

?

基于隨機(jī)模擬法的數(shù)學(xué)實(shí)驗(yàn)設(shè)計(jì)

高文華a，吳培浩b，崔超宇b，張宇風(fēng)a，廖嘉祺a

(華南理工大學(xué)a.數(shù)學(xué)學(xué)院;b.機(jī)械與汽車學(xué)院，廣州510640)

摘要該文基于隨機(jī)模擬法和Matlab編程，設(shè)計(jì)了布朗運(yùn)動(dòng)的模擬和對(duì)伊藤型隨機(jī)微分方程解的數(shù)值仿真兩個(gè)數(shù)學(xué)實(shí)驗(yàn)。設(shè)計(jì)的實(shí)驗(yàn)在豐富數(shù)學(xué)實(shí)驗(yàn)課程的教學(xué)內(nèi)容的同時(shí)，使學(xué)生對(duì)隨機(jī)過(guò)程微分方程有更深入的理解，有利于以創(chuàng)造性思維能力和自學(xué)能力為重點(diǎn)的學(xué)生綜合能力的培養(yǎng)。

關(guān)鍵詞數(shù)學(xué)實(shí)驗(yàn)；Monte Carlo方法；Matlab 仿真；隨機(jī)微分方程

數(shù)學(xué)實(shí)驗(yàn)課是數(shù)學(xué)與計(jì)算機(jī)相結(jié)合的課程，更注重對(duì)學(xué)生應(yīng)用數(shù)學(xué)能力的培養(yǎng)，強(qiáng)調(diào)學(xué)生自己動(dòng)手做數(shù)學(xué)，也更能激發(fā)學(xué)生學(xué)習(xí)數(shù)學(xué)的積極性和興趣。在數(shù)學(xué)實(shí)驗(yàn)中，由于計(jì)算機(jī)的引入和數(shù)學(xué)軟件包的應(yīng)用，為數(shù)學(xué)的思想與方法注入了更多、更廣泛的內(nèi)容，促進(jìn)了數(shù)學(xué)同其他學(xué)科之間的結(jié)合，從而使學(xué)生有時(shí)間去做更多的創(chuàng)造性工作[1]。

隨機(jī)現(xiàn)象在實(shí)際生活中廣泛存在。隨時(shí)間變化的各類系統(tǒng)受到不確定因素的作用，這些不確定因素往往服從某種統(tǒng)計(jì)規(guī)律，通常把這種具有統(tǒng)計(jì)規(guī)律的不確定因素稱為隨機(jī)因素。隨機(jī)系統(tǒng)就是指用于描述這類受隨機(jī)因素作用的時(shí)間過(guò)程的一類數(shù)學(xué)模型。這類數(shù)學(xué)模型一般是某些含隨機(jī)過(guò)程的差分或微分方程[2]。隨機(jī)微分方程能更真實(shí)、更準(zhǔn)確地反映實(shí)際工程技術(shù)中的系統(tǒng)運(yùn)動(dòng)規(guī)律。隨機(jī)控制理論也廣泛地應(yīng)用于經(jīng)濟(jì)、人口系統(tǒng)等社會(huì)領(lǐng)域以及航空航天、導(dǎo)航與控制、制造工程等工程領(lǐng)域。隨機(jī)系統(tǒng)的研究已成為現(xiàn)代控制理論研究中的一個(gè)熱點(diǎn)問(wèn)題。

隨機(jī)模擬法也叫Monte Carlo方法，它是用計(jì)算機(jī)模擬隨機(jī)現(xiàn)象，通過(guò)大量的仿真試驗(yàn)，進(jìn)行分析推斷。概率論中的大數(shù)定理是數(shù)字特征隨機(jī)模擬的理論根據(jù)。本文利用Matlab編程[3]，設(shè)計(jì)了兩個(gè)數(shù)學(xué)實(shí)驗(yàn)：布朗運(yùn)動(dòng)[4-6]的模擬和伊藤隨機(jī)微分方程[7]的解的數(shù)值仿真。通過(guò)這些實(shí)驗(yàn)，可使學(xué)生對(duì)隨機(jī)過(guò)程和隨機(jī)微分方程有更深入的理解，也將激發(fā)學(xué)生學(xué)習(xí)數(shù)學(xué)的興趣。

1 布朗運(yùn)動(dòng)的模擬實(shí)驗(yàn)

布朗運(yùn)動(dòng)是英國(guó)植物學(xué)家布朗在1827年首先發(fā)現(xiàn)的，他觀察到懸浮在液體中的小粒子的擴(kuò)散運(yùn)動(dòng)是極不規(guī)則的。1905年，著名物理學(xué)家愛(ài)因斯坦提出，布朗運(yùn)動(dòng)是由微小粒子與液體分子之間的碰撞引起的。1923年，維納對(duì)布朗運(yùn)動(dòng)給出了嚴(yán)格的數(shù)學(xué)分析。布朗運(yùn)動(dòng)的第一個(gè)嚴(yán)格研究是維納給出的，因此這種隨機(jī)過(guò)程也稱作維納過(guò)程。維納最重要的貢獻(xiàn)是他對(duì)布朗運(yùn)動(dòng)過(guò)程的樣本函數(shù)性質(zhì)的研究，他證明了布朗運(yùn)動(dòng)的樣本函數(shù)或軌道是連續(xù)的而且?guī)缀跆幪幉豢晌⒌暮瘮?shù)。

自布朗運(yùn)動(dòng)發(fā)現(xiàn)以來(lái)，這個(gè)隨機(jī)過(guò)程已在數(shù)理統(tǒng)計(jì)、經(jīng)濟(jì)學(xué)、量子力學(xué)、生物學(xué)、通信理論和管理科學(xué)等領(lǐng)域得到廣泛運(yùn)用。

標(biāo)量布朗運(yùn)動(dòng)又稱標(biāo)準(zhǔn)維納過(guò)程，是指定義在區(qū)間[0，T]上的一個(gè)隨機(jī)變量W(t)，且滿足下述3個(gè)條件：

(1)W(0)=0;

(3)對(duì)于0≤s

圖1　對(duì)標(biāo)量布朗運(yùn)動(dòng)的模擬

T=1;

N=500;

dt=T/N;

dW=zeros(1，N);

W=zeros(1，N);

dW(1)=sqrt(dt)*randn;

W(1)=dW(1);

for j=2:N

W(j)=W(j-1)+dW(j);

end

plot([0:dt:T]，[0，W]，′r-′);

xlabel(′t′，′FontSize′，16);

ylabel(′W(t)′，′FontSize′，16，′Rotation′，0);

可以將上述代碼保存為m文件，程序的第一行為randn(′state′，200)時(shí)，不論在哪臺(tái)計(jì)算機(jī)上運(yùn)行得到的結(jié)論都如圖1所示。如果把這個(gè)語(yǔ)句換為randn(′state′，sum(200*clock))，進(jìn)行盡可能多次的模擬，會(huì)發(fā)現(xiàn)每一次的模擬得到的圖形都不一樣，都與之前的模擬無(wú)關(guān)，即體現(xiàn)了布朗運(yùn)動(dòng)的隨機(jī)性。

2伊藤隨機(jī)微分方程的解的數(shù)值仿真實(shí)驗(yàn)

伊藤型隨機(jī)微分方程在工業(yè)、社會(huì)、經(jīng)濟(jì)領(lǐng)域都有廣泛的應(yīng)用，如在金融領(lǐng)域描述收益率波動(dòng)的期權(quán)定價(jià)數(shù)學(xué)模型就是采用伊藤型隨機(jī)微分方程，在工業(yè)領(lǐng)域?yàn)槊枋霭自肼晫?duì)系統(tǒng)的影響也是采用伊藤型隨機(jī)微分方程。

對(duì)于隨機(jī)微分方程而言，一般很難求出其具體解過(guò)程。因此，許多學(xué)者提出不同的數(shù)值策略，如Euler-Maruyama方法，向前、向后Euler-Maruyama方法，Taylor展式方法等。結(jié)合布朗運(yùn)動(dòng)的模擬，本文中主要運(yùn)用Euler-Maruyama方法對(duì)伊藤隨機(jī)微分方程的解進(jìn)行數(shù)值仿真。

2.1簡(jiǎn)單的伊藤型隨機(jī)微分方程解的數(shù)值仿真

例1簡(jiǎn)單的伊藤型隨機(jī)微分方程

實(shí)驗(yàn)要求是使用Matlab編程對(duì)上述隨機(jī)微分方程給出方程解的仿真曲線。

仿真程序如下，程序運(yùn)行結(jié)果如圖2所示。

randn(′state′，500)

A=[1.1，-1.3;1.8，-2];X1zero = 1;

T = 1;N = 2^8;dt = 1/N;

randn(′state′，400);

dWy=sqrt(dt)*randn(1，N);

X(:，1)=1;

for j = 1:L-1

Winc = [sum(dWx(R*(j-1)+1:R*j));sum(dWy(R*(j-1)+1:R*j))];

X(:，j+1) = X(:，j)+Dt*A*X(:，j)+ B*Winc;

end

plot(0:Dt:T，[X(1，1)，X(1，:)]，′r--*′)，hold on

plot(0:Dt:T，[1，X(2，:)]，′g--+′)，hold off

xlabel(′t′，′FontSize′，12)

ylabel(′X′，′FontSize′，16，′Rotation′，0，′HorizontalAlignment′，′right′)

圖2　簡(jiǎn)單的伊藤型隨機(jī)微分方程解的仿真曲線

2.2帶有時(shí)滯的伊藤型隨機(jī)微分方程解的數(shù)值仿真

例2帶有時(shí)滯的伊藤型隨機(jī)微分方程

實(shí)驗(yàn)要求使用Matlab編程對(duì)上述隨機(jī)微分方程給出方程解的仿真曲線。

仿真程序如下，程序運(yùn)行結(jié)果如圖3所示。

圖3　帶有時(shí)滯的伊藤型隨機(jī)微分方程解的仿真曲線

randn(′state′，200)

A=[-0.7，-1.2;-0.5，-0.2];

B=[-1.6，-1.8;2.8，-1.4];

A1=[2，-0.9;-0.9，0.8];B1=[1.4，-2.5;0.07，-1.5];

T = 1;N = 2^8;dt = 1/N;

Xem = zeros(1，L+5);

Yem=zeros(1，L+5);

X=[Xem;Yem];

for j = 1:L+5

if j<6;

X(1，j)=sin(j-5)+1;

X(2，j)=sin(2*(j-5))+1;

else

Winc = sum(dW(R*(j-6)+1:R*(j-5)));

[X(:，j+1)]=X(:，j)+C*X(:，j+1-4)-C*X(:，j-4)+(A*X(:，j)+A1*X(:，j-4))*Dt+(B*X(:，j)+B1*X(:，j-4))*Winc;

end

end

plot(0:Dt:T，X(1，(5:L+5))，′r--*′)，hold on;

plot(0:Dt:T，X(2，(5:L+5))，′b--+′)，hold off;

legend(′ x1圖像′，′ x2圖像′，1);

xlabel(′t′，′FontSize′，12);

ylabel(′x′，′FontSize′，16，′Rotation′，0，′HorizontalAlignment′，′right′)

3結(jié)束語(yǔ)

數(shù)學(xué)實(shí)驗(yàn)課程能夠提高學(xué)生學(xué)習(xí)數(shù)學(xué)的積極性，培養(yǎng)學(xué)生應(yīng)用數(shù)學(xué)的能力?；诳蒲信c教學(xué)資源共享的原則，本文把科研中用到的仿真算法和例子編成了數(shù)學(xué)實(shí)驗(yàn)的案例，可以拓寬學(xué)生的視野，豐富數(shù)學(xué)實(shí)驗(yàn)課程的教學(xué)內(nèi)容。

參 考 文 獻(xiàn)

[1]焦光虹.數(shù)學(xué)實(shí)驗(yàn)[M].北京：科學(xué)出版社，2010.

[2]韓崇昭，王月娟，萬(wàn)百五.隨機(jī)系統(tǒng)理論[M].西安：西安交通大學(xué)出版社，1987.

[3]趙景波.Matlab控制系統(tǒng)仿真與設(shè)計(jì)[M].北京:機(jī)械工業(yè)出版社，2010.

[4]BRZEZNIAK Z，ZASTAWNIAK T.Basic stochastic processes[M].影印版.北京：清華大學(xué)出版社，2009.

[5]陳六新.概率論與隨機(jī)過(guò)程[M].北京：清華大學(xué)出版社，2013.

[6]伊藤清.隨機(jī)過(guò)程[M].上海：上?？茖W(xué)技術(shù)出版社，1961.

[7]HIGHAM D J.An algorithmic introduction to numerical simulation of stochastic differential equations[J].Siam Review.2001，43(3):525-546.

收稿日期:2015-03-05；修改日期: 2015-03-25

基金項(xiàng)目：廣東省高等學(xué)校教學(xué)質(zhì)量與教學(xué)改革項(xiàng)目(x2lxN9120650);華南理工大學(xué)教研重點(diǎn)項(xiàng)目(x2lxY1130010)。

作者簡(jiǎn)介：高文華(1974-) ，女，博士，副教授，主要從事隨機(jī)時(shí)滯系統(tǒng)的穩(wěn)定性與控制方面的研究。

中圖分類號(hào)O175;TP391.9

文獻(xiàn)標(biāo)志碼A

doi:10.3969/j.issn.1672-4550.2016.03.003

Mathematical Experiment Design Based on Random Simulation Method

GAO Wenhuaa，WU Peihaob，CUI Chaoyub，ZHANG Yufenga，LIAO Jiaqia

(a.School of Mathematics；b.School of Mechanical and Automotive Engineering，South China University of Technology，Guangzhou 510640，China)

AbstractBased on Monte Carlo method and Matlab programming，two mathematical experiment cases are designed in this paper.One is the simulation of the Brownian movement，the other is the numerical simulation of the It? stochastic differential equation.These cases will enrich the content of mathematical experiment course as well as deepen the understanding of differential equation of random process.It would help the development of comprehensive ability of the students centering on creative thought ability and self-learning ability.

Key wordsmathematical experiment;Monte Carlo method;Matlab simulation;stochastic differential equation