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

作业一、 蒙特卡洛法求定积分：𝜃 =∫40(cos𝑥 + 2.0)𝑑𝑥 整体思路，蒙特卡洛求定积分的主要思想就是通过在取值范围内大量随机数的随机选取对 函数进行求值，进而除以所取次数并乘以其区间，即为所积值 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 countNew—>Function,并点击Function，弹 出Editor窗口，将其中内容修改为 将二阶函数化为一阶过程中，定义如下：将二阶函数化为一阶过程中，定义如下： function dy=func(~,y)dy=zeros(2,1); %初始化，2行一列，即（）=（ ）𝑑𝑦(1) 𝑑𝑦(2)0 0dy(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坐标024681012-10-5051015值 值 值xy方法二： 2.2.最常用的最常用的 R-KR-K 公式公式 ———— 标准标准 4 4 阶阶 R-KR-K 公式为：公式为：{𝑦𝑖 + 1= 𝑦𝑖+ℎ 6(𝐾1+ 2𝐾2+ 2𝐾3+ 𝐾4)𝐾1= 𝑓(𝑥𝑖+ 𝑦𝑖)𝐾2= 𝑓(𝑥𝑖+ℎ 2,𝑦𝑖+ℎ 2𝐾1)𝐾3= 𝑓(𝑥𝑖+ℎ 2,𝑦𝑖+ℎ 2𝐾2) 𝐾4= 𝑓(𝑥𝑖+ ℎ,𝑦𝑖+ ℎ𝐾3)?在在editoreditor窗口中打开一个新函数，窗口中打开一个新函数，对龙格库塔函数定义：对龙格库塔函数定义： 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)且步长为hk=1;%初始化k=1y(:,1)=y0(:); %初始化y(:,1)for i=1:length(x)-1 %从第一次循环开始，共循环(length(x)-1)次K1=feval(odefun,x(k),y(:,k));%feval(函数名，x1，x2，…)K2=feval(odefun,x(k)+h/2,y(:,k)+h/2*K1);K3=feval(odefun,x(k)+h/2,y(:,k)+h/2*K2);K4=feval(odefun,x(k)+h,y(:,k)+h*K3); %以上K1，K2，K3，K4均为龙哥库塔定义中的式子y(:,k+1)=y(:,k)+h/6*(K1+2*K2+2*K3+K4); %利用龙哥库塔定义通过y(:,k)的迭代求y(:,k+1)k=k+1;endend回到Matlab程序中，将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。

在Command Window里输入以下源程序 >> f=@(x,y)[y(2);-4*y(1)-29*y(1)];%利用刚刚化简的一次方程 >> [x,y]=lgkt4j(f,[0,12],[0,15],1)%利用定义的龙哥库塔，[0,12]为x区间，[0,0]为（）=（ ），1为𝑑𝑦(1) 𝑑𝑦(2)0 0步长二、求微分方程的特征解：关于求微分方程（组）的解析解命令: dsolve(‘方程 1’, ‘方程 2’,…‘方程 n’, ‘初始条件’, ‘自变量’) 关于此公式的注解：在表达微分方程时，用字母 D 表示求微分,D2、D3 等可以表示求高阶 微分，任何 D 后所跟字母为因变量，自变量可以指定或由系统规则选定为缺省 步骤：在 matlab 中直接输入源程序： 源程序： y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') %直接利用解析函数求特解 y =(3*sin(5*x))/exp(2*x) %结果。