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

作业一、 蒙特卡洛法求定积分： θ=04(cosx+2.0)dx 整体思路，蒙特卡洛求定积分的主要思想就是通过在取值范围内大量随机数的随机选取对函数进行求值，进而除以所取次数并乘以其区间，即为所积值。 Step 1： 在执行程序前，打开matlab，执行以下操作打开File—New—Function,并点击Function，弹出Editor窗口，将其中内容修改为 function [ y ] = f( x ) y=cos(x)+2.0; end （如图所示）。 Step 2： 在Editor窗口中点击File—Save As，保存至所建立的文件夹内，保存名称必须为f.m。 Step 3： 回到Matlab程序中，将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。 Step 4： 在Command Window里输入以下源程序。 源程序： n=10^6; %用来表示精度，表示需要执行次数。 y=0; %初始化y=0，因为f(x)=cos(x)+2.0，下面出现y=y+f(x),故尔y用来统计计算出的所有f(x)值的和。 count=1; %从1开始计数，显然要执行n次。 while count=n %当执行次数少于n次时，继续执行 x=unifrnd(0,4); % 在matlab中，unifrnd是在某个区间内均匀选取实数（可为小数或整数），表示x取0到4的任意数。 y=y+f(x); %求出到目前为止执行出的所有f(x)值的和，并用y表示 count=count+1; %count+1表示已执行次数再加一，将来通过与n的比对，判断是否执行下一次 end %如果执行次数已达n次，那么结束循环 z=y/n*4 %用y除以执行次数n,那么取得平均值,并用它乘以区间长度4,即可得到z的近似值 z=7.2430 由于matlab中不能使用θ字符，故用z代替，即可得到θ=7.2395。 在执行过程中，发现每一次执行结果都会不一样，这是因为是随机选取，所以不同次数的计算结果会有偏差，但由于执行次数很多，从概率的角度来讲都是与真实值相近似的值。 作业二、 用四阶龙格库塔法解微分方程： d2ydx2+4dydx+29y=0y0=0,y0=15 解： 求微分方程的数值解： 方法一： 1.将方程化为1次，即化为： y=y1; y2=y1y2+4y2+29y1=0y10=0; y20=15 在执行程序前，打开matlab，执行以下操作打开File—New—Function,并点击Function，弹出Editor窗口，将其中内容修改为 将二阶函数化为一阶过程中，定义如下： function dy=func(~,y) dy=zeros(2,1); %初始化，2行一列，即（dy(1)dy(2)）=（00） dy(1)=y(2); %将上述方程组化简得来 dy(2)=-4*y(2)-29*y(1); %将上述方程组化简得来 end %定义结束 回到Matlab程序中，将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。 在Command Window里输入以下源程序。 [x,y]=ode45(func,[0 12],[0 15]) 绘图 y=size(y) y = 229 2% 数组y包含的元素数 reshape(y,458,1)%将含有229*2=458格元素的数组y化为458行1列的数组 plot(x,y) %绘图 title(数值图)%图名 hold on%保持当前图形 xlabel(x)%加x坐标 ylabel(y)%加y坐标 方法二： 2.最常用的R-K公式 —— 标准4阶R-K公式为： yi+1=yi+h6(K1+2K2+2K3+K4)K1=f(xi+yi)K2=f(xi+h2,yi+h2K1)K3=f(xi+h2,yi+h2K2)K4=f(xi+h,yi+hK3) 在editor窗口中打开一个新函数，对龙格库塔函数定义： function [x,y]=lgkt4j(odefun,xspan,y0,h) %odefun为微分方程的函数描述，xspan为解区间[x0,xn],y0为初始条件，h为迭代步长 x=xspan(1):h:xspan(2); %定义x为从xspan(1)开始到xspan(2)且步长为h k=1;%初始化k=1 y(:,1)=y0(:); %