蒙特卡洛求定积分及龙哥库塔解微分方程

蒙特卡洛求定积分及龙哥库塔解微分方程

ID:42045978

大小:135.45 KB

页数:4页

时间:2019-09-06

蒙特卡洛求定积分及龙哥库塔解微分方程_第1页
蒙特卡洛求定积分及龙哥库塔解微分方程_第2页
蒙特卡洛求定积分及龙哥库塔解微分方程_第3页
蒙特卡洛求定积分及龙哥库塔解微分方程_第4页
资源描述:

《蒙特卡洛求定积分及龙哥库塔解微分方程》由会员上传分享,免费在线阅读,更多相关内容在工程资料-天天文库

1、作业一、蒙特卡洛法求定积分:0=(cosx+2.0)dx整体思路,蒙特卡洛求定积分的主要思想就是通过在取值范围内大量随机数的随机选取对函数进行求值,进而除以所取次数并乘以其区间,即为所积值。Step1:FileEditTextGoCellToolsDebugDesktc♦o回1.0+■1.1X必癣O1Efunction[y3=f(X)2尸cos(x)+2.0:1•3end:'Editor-Untitled*在执行程序前,打开matlab,执行以下操作打开File—>New—>Functioix并点击Function,弹出Editor窗口,将其中内容修改为function[y]=f(x

2、)y=cos(x)+2•0;end(如图所示)。Step2:在Editor窗口中点击File->SaveAs,保存至所建立的文件夹内,保存名称必须为f.m。Step3:冋到Matlab程序中,将CurrentFolder改为与刚刚Function函数定义的保存路径一致的路径°Step4:在CommandWindow里输入以卜•源程序。源程序:»n=10A6;»y=0;%用来表示精度,表示需要执行次数。%初始化y=0,因为f(x)二cos(x)+2.0,下面出现y=y+f(x),故尔y用来统计计算出的所有f(x)值的和。»count=l;»whilecountv二nx=unifrnd(

3、0,4);%从1幵始计数,显然要执行n次。%当执行次数少于n次时,继续执行%在matlab中,unifrnd是在某个区间内均匀选取实数(可为小数或整数),表示x取0到4的任意数。y=y+f(x);count=count+1;%求出到目前为止执行出的所有f(x)值的和,并用y表示%count+l表示已执行次数再加一,将来通过与n的比对,判断是否执行下一次endz=y/n*4%如果执行次数已达n次,那么结束循环%用y除以执行次数①那么取得平均值,并用它乘以区间长度4,即可得到z的近似值z=7.2430由于matlab中不能使用&字符,故用z代替,即可得到o=7.2395o在执行过程中,发

4、现每一次执行结果都会不一样,这是因为是随机选取,所以不同次数的计算结果会有偏差,但由于执行次数很多,从概率的角度来讲都是与真实值相近似的值。作业二、用四阶龙格库塔法解微分方程:fd2yNew—'Function,并点击Function,弹出Editor窗口,将其中内容修改为将二阶函数化为一阶过程中,定义如下:functiondy=f

5、unc(~,y)dy=zeros(2,1);密初始化,2行一列,即(;器)=C)dy(l)=y(2);%将上述方程组化简得来dy(2)=-4*y(2)-29*y(1);拓将上述方程组化简得来end%定义结束回到Matlab程序中,将CurrentFolder改为与刚刚Functional数定义的保存路径一致的路径。在CommandWindow里输入以下源程序。>>[xzy]=ode45(*func1,[012]z[015])绘图>>y=size(y)y=2292%数组y包含的元素数>>reshape(y,458,1)%将含有229*2=458格元素的数组y化为458行1列的数组>>p

6、lot(x,y)%绘图>>title(1数值图T%图名>>holdon%保持当前图形>>xlabel%力Ux坐标>>ylabel(fy*)%加丫坐标数值图15-1%281012方法二2.最常用的R-K公式标准4阶R-K公式为:Yi+i=并+石(K」+2K2+2K3+KJKi=f(Xi+yi)"fg雋,刃将Kt)K4=f(Xi+h,yi+hK3)在editor窗口中打开一个新函数,对龙格库塔函数定义:function[x,y]=lgkt4j(odefun,xspan,yO,h)%odefun为微分方程的函数描述,xspan为解区间[xO,xn],y0为初始条件,h为迭代步长x=xspa

7、n(1):h:xspan(2);%定义区为从xspan(1)开女台至Ijxspan(2)且步长为hk=l;%初始化k=ly(:zl)=yO(:);%初始化y(:,l)fori=l:length(x)-1为从第一次循坏开始,共循坏(length(x)-1)次Kl=feval(ode£un,x(k),y(:,k));%feva丄(函数名,xl,x2,...)K2=feval(odefunzx(k)+h/2,y(:rk)+h/2*Kl);K3=feval(ode

当前文档最多预览五页,下载文档查看全文

此文档下载收益归作者所有

当前文档最多预览五页,下载文档查看全文
温馨提示:
1. 部分包含数学公式或PPT动画的文件,查看预览时可能会显示错乱或异常,文件下载后无此问题,请放心下载。
2. 本文档由用户上传,版权归属用户,天天文库负责整理代发布。如果您对本文档版权有争议请及时联系客服。
3. 下载前请仔细阅读文档内容,确认文档内容符合您的需求后进行下载,若出现内容与标题不符可向本站投诉处理。
4. 下载文档时可能由于网络波动等原因无法下载或下载错误,付费完成后未能成功下载的用户请联系客服处理。