下载文档原格式
计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。
它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。
龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。
在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。
龙格库塔法的核心是使用差分逼近计算函数的斜率。
假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。
龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。
接下来以龙格库塔法的经典四阶形式为例进行说明。
该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。
使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。
具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。
2.初始化：设置yy=y₀，yy=y₀。
3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。
根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。
更新当前点的y值，即yy+1=yy+ℎ。
4.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。
欢迎阅读考试课程 数学实验下午班级 姓名 学号 得分[说明]（1）第一、二、三题的答案直接填在试题纸上；（2）第四题将数学模型、简要解题过程和结果写在试题纸上；卷面空间不够时，请写在背面；（3）除非特别说明，所有计算结果小数点后保留4位数字。
（4）考试时间为120分钟。
一、（10分）某厂生产A 、B 两种产品，1千克原料在甲类设备上用12小时可生产3件A ，可获净利润64元；在乙类设备上用8小时可生产4件B ，可获净利润54元。
该厂每天可获得55千克原料，每天总的劳动时间为480小时，且甲类设备每天至多能生产80件A 。
试为该厂制订生产计划1 生产s.t. c=[64 A1=[ [x,z,ef,out,lag]=linprog(-c,A1,b1,[],[],v1)lag.ineqlin输出结果:x =10.000000004005848ans =0.0000000002784052）每天的最大净利润是___3070__元。
若要求工人加班以增加劳动时间，则加班费最多 为每小时__2.5__元。
若A 获利增加到26元/件，应否改变生产计划？____不变___c=[78 54];A1=[ 1 1 ;12 8 ;3 0];b1=[55;480;80];v1=[0 0];[x,z,ef,out,lag]=linprog(-c,A1,b1,[],[],v1)45.000000000000625二、(10分) 已知常微分方程组初值问题试用数值方法求()6yπ=__ 1.73205____(保留小数点后5位数字)。
你用的MATLAB命令是______ode45(@ff, ts,y0)______，其精度为____四阶__。
%待解常微分方程组函数M文件源程序：function dy=ff (x,y)dy=[y(2);-y(2)./x-y(1)*(x.^2-0.25)/(x.^2)];%应用欧拉方法和龙格-库塔方法求解该常微分方程：三、A分别-赛德敛___d=cond(A,1)*norm(db,1)/norm(b,1)输出结果:A=[5 -7 0 1 ;-3 22 6 2 ;5 -1 31 -1 ;2 1 0 23];D=diag(diag(A)); %从稀疏矩阵A中提取DL=-tril(A,-1); %从稀疏矩阵A中提取LU=-triu(A,1); %从稀疏矩阵A中提取Ub=[6 3 4 7]'; %设定方程组右端项向量x= zeros(4,1); %设定方程组初始向量m= inv(D-L)*U;n= inv(D-L)*b; %高斯-赛德尔迭代法for j2=1:5y=m*(x(:,j2));for i=1:4x(i,j2+1)=y(i,:)+n(i,:);endendt2=x(:,end) %输出迭代法最终结果j2输出结果:t2 =判敛:lamda=eig(inv(D-L)*U)pubanjing=max(abs(lamda))输出结果:四、（20分）炮弹射击的目标为一椭圆形区域，在X方向半轴长110m，Y方向半轴长90m.当瞄准k=(n-1)*var(x)/(s0^2) %χ2分布检验方差if tail==0k1=chi2inv(alpha/2,n-1)k2=chi2inv(1-alpha/2,n-1)if k>=k1&k<=k2h=0;elseh=1;endendif tail==1k0=chi2inv(1-alpha,n-1)if k<=k0h=0;elseh=1;endendif tail==-1k0=chi2inv(alpha,n-1)if k>=k0h=0;elseh=1;endh1=ktest(x,70,0.05,0)h2=ktest(y,50,0.05,0)输出结果：h1 =0h2 =02)15.2];r=3)%b=90;z=0;u=exp(-0.5/(1-r^2)*(x(i)^2/sx^2-2*r*x(i)*y(i)/(sx*sy)+y(i)^2/sy^2));z=z+u;endendP=4*a*b*z/(2*pi*sx*sy*sqrt(1-r^2)*n)输出结果：P =考试课程数学实验下午班级学号姓名得分[说明]（1）第一、二、三题的答案直接填在试题纸上；（2）第四题将数学模型、简要解题过程和结果写在试题纸上；卷面空间不够时，请写在背面；（3）除非特别说明，所有计算结果小数点后保留4位数字。
四阶龙格库塔法计算例题哎呀，宝子们！今天咱们就来唠唠四阶龙格 - 库塔法计算例题这事儿。
这四阶龙格 - 库塔法啊，就像是我们数学世界里的一个小魔法。
咱先来说说它是咋来的呢？这可就得追溯到数学发展的长河里了。
这方法呀，是经过好多数学家的智慧结晶才形成的。
它就像是一个精心打造的工具，专门用来解决一些特定的数学计算问题。
那它到底怎么计算呢？比如说有个例题，给了一个一阶常微分方程。
咱就把这个方程想象成是一个小怪兽，而四阶龙格 - 库塔法就是咱们降伏小怪兽的魔法棒。
我们先把方程写成标准形式，就好像是给小怪兽定个型。
然后呢，根据这个方法的公式，一步一步地来计算。
这个过程就像是在给小怪兽套上各种枷锁，让它乖乖听话。
我给大家举个特别简单的例子哈。
假设我们有个方程是dy/dx = x + y，初始条件是y(0)=1。
这时候我们就可以运用四阶龙格 - 库塔法啦。
我们要先确定步长，这个步长就像是我们前进的小步伐，可不能太大也不能太小哦。
比如说我们设步长h = 0.1。
然后按照公式计算k1，k1就像是我们对小怪兽发起的第一轮攻击。
k1 = h f(xn,yn)，这里的f(xn,yn)就是我们方程的右边部分，也就是x + y。
所以k1 = 0.1 (0 + 1)=0.1。
接着算k2，k2就像是我们调整了攻击策略后的又一轮进攻。
k2 = h f(xn + h/2,yn + k1/2)。
代入数值算出来。
再算k3，k3的计算方式和k2有点像，但又有一点点小区别。
最后算k4。
最后呢，yn+1 = yn + (k1 + 2k2 + 2k3 + k4)/6。
就这样一步一步地，我们就能算出y在不同点的值啦。
这四阶龙格 - 库塔法啊，虽然计算的时候可能有点小复杂，但是一旦掌握了，就会觉得超级有趣。
就像是我们学会了一个超级厉害的游戏秘籍一样。
每次算出正确的结果，都有一种成就感油然而生呢。
而且这个方法在很多实际的科学计算中都超级有用，比如说物理里的一些运动方程的求解，工程里的一些模型计算等等。
例题1.11 用龙格-库塔法求解下列微分方程组：⎪⎪⎩⎪⎪⎨⎧⨯--==1222212.31.0dx x x dtdx x dt 当t=0时，1x =1，2x =0.求解区间为（0，1），步长为0.2。
C 语言程序如下：#include<stdio.h>#include<math.h>#include<stdlib.h>main(){int i,a,b;double h,n,y1,y2,x[100],z[100],y[100],k1,k2,k3,k4,g1,g2,g3,g4;printf(" 请输入区间（a,b ）:\n\t");scanf("%d %d",&a,&b);printf(" 请输入步长h:\n\t");scanf("%lf",&h);printf(" 请输入初始条件x1,x2:\n\t");scanf("%lf %lf",&y1,&y2);n=(b-a)/h;y[0]=y1;z[0]=y2;x[0]=a;for(i=0;i<=n;i++){x[i+1]=x[i]+h;k1=z[i];k2=y[i]+k1*h*0.5;k3=y[i]+k2*h*0.5;k4=y[i]+k3*h;g1=-0.1*z[i]-3.2*3.2*y[i];g2=z[i]+g1*h*0.5;g3=z[i]+g2*h*0.5;g4=z[i]+g3*h;z[i+1]=z[i]+h*(g1+2*g2+2*g3+g4)/6;y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}printf("\tT\tX0\tX1\n");for(i=0;i<=n;i++){printf("\t%1.3f\t%lf\t%lf\n",x[i],y[i],z[i]);}system("pause");}其运行结果如下：例题3.2 用黄金分割法求解下列最优化问题：目标函数：2016)(2+-=x x x f约束条件：⎩⎨⎧≤-≥+010010x x 根据约束条件，最优x 只能在区间 [-10,10]搜索。
实验 龙格—库塔方法实验目的： 深入理解龙格—库塔方法的原理，同时学会用迭代方法看待某些特定问题。
比较本方法与其它方法（尤其是欧拉方法）解题的不同处，分析两者的精度。
实验要求: 编写程序，用龙格—库塔方法求解第三章3题（修改）实验内容：题目： 取h =0.1，用四阶龙格-库塔方法求解初值问题： 41,211)(22≤≤-+='x y xx f 0)0(=y 并与精确解21x x y +=比较结果. 原理： 本题很明显为迭代方法的一种。
已知初值和函数的导数。
再根据四阶龙格—库塔公式：k 1=dy(x0,y0);k 2=dy((2 * x0 + h)/2 , y0 + h/2 * k 1);k 3=dy((2 * x0 + h)/2 , y0 + h/2 * k 2);k 4=dy(x0 + h , y0 + h * k 3);y0=y0 + h/6 * (k 1 + 2*k 2 + 2*k 3 + k 4);x0=x0+h;设计思想：由四阶龙格—库塔公式很容易看出，只要作一个循环即可完与此项操作。
每次更新k 1,k 2,k 3,k 4,y0以及x0的值。
这样即可得出y 。
源代码：disp('第三章3题：取h=0.1，用四阶龙格-库塔方法求解初值问题：')disp(' diff(y,x,1)=1/(1+x.^2)-2*y.^2 1<=x<=4') disp(' y(0)=0 ')disp(' 并与精确解y=x/(1+x.^2)比较结果.')x0=1; %初值为1y0=0.5; %初值为0.5h=input('请输入步长>>')for i=1:1:5k1=dy(x0,y0);k2=dy((2 * x0 + h)/2 , y0 + h/2 * k1);k3=dy((2 * x0 + h)/2 , y0 + h/2 * k2);k4=dy(x0 + h , y0 + h * k3);y0=y0 + h/6 * (k1 + 2*k2 + 2*k3 + k4);x0=x0+h;fprintf('\nk%d = %f ,k%d = %f ,k%d = %f ,k%d = %f ',1,k1,2,k2,3,k3,4,k4) fprintf('\n =>>y[%d]=%f ',i,y0);endfprintf('\n 以下为用公式：y=x/(1+x.^2)求得的精确值\n')x0=1;for i=0:1:5y0=f(x0);x0=x0+h;fprintf('y[%d]=%f ',i,y0);endfprintf('\n/****************************************************/')(下图为yy=dy(x)以及y=f(x)两函数的源代码：)文件名： f.mfunction y=f(x)y=x./(1+x.^2);文件名：dy.mfunction yy=dy(x,y)yy=1/(1+x.^2)-2*y.^2;实验结果：图形实验体会：本实验中，关于龙格—库塔方法的使用，使我个人又对迭代方法有了新的一种理念。
考试课程数学实验下午班级姓名学号得分[说明]（1）第一、二、三题的答案直接填在试题纸上；（2）第四题将数学模型、简要解题过程和结果写在试题纸上；卷面空间不够时，请写在背面；（3）除非特别说明，所有计算结果小数点后保留4位数字。
（4）考试时间为120分钟。
一、（10分）某厂生产A、B两种产品，1千克原料在甲类设备上用12小时可生产3件A，可获净利润64元；在乙类设备上用8小时可生产4件B，可获净利润54元。
该厂每天可获得55千克原料，每天总的劳动时间为480小时，且甲类设备每天至多能生产80件A。
试为该厂制订生产计划使每天的净利润最大。
1）以生产A、B产品所用原料的数量x1、x2（千克）作为决策变量，建立的数学规划模型是：决策变量：生产A原料x1；生产B原料x2目标函数：y=64*x1+54*x2约束条件：x1+x2 ≤5512*x1+8*x2≤4803*x1≤80x1,x2≥0基本模型：max(y)= 64*x1+54*x2s.t. x1+x2 ≤5512*x1+8*x2≤4803*x1≤80x1,x2≥0c=[64 54];A1=[ 1 1 ;12 8 ;3 0];b1=[55;480;80];v1=[0 0];[x,z,ef,out,lag]=linprog(-c,A1,b1,[],[],v1)lag.ineqlin输出结果:x =10.000000004005848ans =0.0000000002784052）每天的最大净利润是___3070__元。
若要求工人加班以增加劳动时间，则加班费最多为每小时__2.5__元。
若A获利增加到26元/件，应否改变生产计划？____不变___c=[78 54];A1=[ 1 1 ;12 8 ;3 0];b1=[55;480;80];v1=[0 0];[x,z,ef,out,lag]=linprog(-c,A1,b1,[],[],v1) 45.000000000000625二、(10分) 已知常微分方程组初值问题试用数值方法求()6y π=__ 1.73205____(保留小数点后5位数字)。
资料范本本资料为word版本，可以直接编辑和打印，感谢您的下载龙格库塔方法及其matlab实现地点：__________________时间：__________________说明：本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。
MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MatLab 是英文MATrix LABoratory（矩阵实验室）的缩写。
在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
关键词：龙格-库塔 matlab 微分方程前言1.1：知识背景龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。
该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。
Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。
经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。
清华数学实验复习试题八(蒙特卡罗方法-龙格-库塔方法)考试课程数学实验2005.6.15下午班级姓名学号得分[说明]（1）第一、二、三题的答案直接填在试题纸上；（2）第四题将数学模型、简要解题过程和结果写在试题纸上；卷面空间不够时，请写在背面；（3）除非特别说明，所有计算结果小数点后保留4位数字。
（4）考试时间为120分钟。
一、（10分）某厂生产A、B两种产品，1千克原料在甲类设备上用12小时可生产3件A，可获净利润64元；在乙类设备上用8小时可生产4件B，可获净利润54元。
该厂每天可获得55千克原料，每天总的劳动时间为480小时，且甲类设备每天至多能生产80件A。
试为该厂制订生产计划使每天的净利润最大。
1）以生产A、B产品所用原料的数量x1、x2（千克）作为决策变量，建立的数学规划模型是：决策变量：生产A原料x1；生产B原料x2目标函数：y=64*x1+54*x2约束条件：x1+x2 ≤5512*x1+8*x2≤4803*x1≤80x1,x2≥0基本模型：max(y)= 64*x1+54*x2s.t. x1+x2 ≤5512*x1+8*x2≤4803*x1≤80x1,x2≥0c=[64 54];A1=[ 1 1 ;12 8 ;3 0];b1=[55;480;80];v1=[0 0];[x,z,ef,out,lag]=linprog(-c,A1,b1,[],[],v1)lag.ineqlin输出结果:x =10.00000000400584844.999999993870908z =-3.069999999925403e+003ans =33.9999999989193572.5000000001404410.0000000002784052）每天的最大净利润是___3070__元。
若要求工人加班以增加劳动时间，则加班费最多为每小时__2.5__元。
若A获利增加到26元/件，应否改变生产计划？____不变___c=[78 54];A1=[ 1 1 ;12 8 ;3 0];b1=[55;480;80];v1=[0 0];[x,z,ef,out,lag]=linprog(-c,A1,b1,[],[],v 1)x = 9.999999999999400 45.000000000000625z =-3.209999999999987e+003二、(10分) 已知常微分方程组初值问题221'''()0,4x y xy x y ++-=()2,2y π=2'().2y ππ=-试用数值方法求()6y π=__ 1.73205____(保留小数点后5位数字)。
你用的MATLAB 命令是______ ode45(@ff, ts,y0)______，其精度为____四阶__。
%待解常微分方程组函数M 文件源程序：function dy=ff (x,y) dy=[y(2);-y(2)./x-y(1)*(x.^2-0.25)/(x.^2)];%应用欧拉方法和龙格-库塔方法求解该常微分方程：ts=pi/2:- pi/12:pi/6; ！！！！步长必须是可以整除步长区间长度的数 y0=[2,-2/pi];[x,y]=ode45(@ff, ts,y0); %龙格-库塔方法求数值解 [x, y(:,1)] 输出结果：0.523598775598299 1.732050795523993 三、 (10分) 已知线性代数方程组Ax =b , 其中5701322625131121023A -⎡⎤⎢⎥-⎢⎥=⎢⎥--⎢⎥⎣⎦, ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=4321x x x x x ,6347b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,若方程组右端项有小扰动]1.0,0,0,0['=b δ，试根据误差估计式估计11x x δ≤___0.0743___（x x δ,分别表示原问题的解和右端项小扰动后对应的解的变化量）；若取初值]0,0,0,0[)0('=x ，则用高斯-赛德尔迭代法求解Ax =b 时，=)5(x _(1.7160, 0.3926, -0.1306, 0.1381)_；对本题而言，此迭代方法是否收敛___是__，原因是__谱半径ρ(B)=0.397<1__。
线性代数方程组解的误差分析：1()x b bA A cond A x b bδδδ-≤⋅⋅=⋅故其误差上限为：A=[5 -7 0 1 ;-3 22 6 2 ;5 -1 31 -1 ;2 1 0 23]; b=[6 3 4 7]; db=[0 0 0 0.1];d=cond(A,1)*norm(db,1)/norm(b,1) 输出结果:d =0.074339065208930A=[5 -7 0 1 ;-3 22 6 2 ;5 -1 31 -1 ;2 1 0 23];D=diag(diag(A)); %从稀疏矩阵A 中提取DL=-tril(A,-1); %从稀疏矩阵A 中提取LU=-triu(A,1); %从稀疏矩阵A 中提取U b=[6 3 47]';%设定方程组右端项向量x=zeros(4,1);%设定方程组初始向量m= inv(D-L)*U;n=inv(D-L)*b;%高斯-赛德尔迭代法for j2=1:5y=m*(x(:,j2));for i=1:4x(i,j2+1)=y(i,:)+n(i,:);endendt2=x(:,end) %输出迭代法最终结果j2输出结果:t2 =1.7159723472264450.392646824062879-0.130571*********0.138061238325401判敛:lamda=eig(inv(D-L)*U)pubanjing=max(abs(lamda))输出结果:pubanjing =0.396832340862002四、（20分）炮弹射击的目标为一椭圆形区域，在X方向半轴长110m，Y方向半轴长90m.当瞄准目标的中心发射炮弹时，在众多随机因素的影响下，弹着点服从以目标中心为均值的二维正态分布，设弹着点偏差的均方差在X方向和Y方向分别为70m和50m。
今测得一组弹着点的横纵坐标如下：X -6.3 -71.6 65.6-79.2-49.7-81.974.6-47.6-120.856.9Y28. 9 1.6 61.7-68 -41.3-30.587 17.3-17.81.2X 100 .9 47 9.7 -60.1-52.786 80.6-42.656.4 15.2Y -12. 6 39.185 32.728.1-9.3-4.55.1 -32 -9.51)根据这组数据对X方向和Y方向的均值和均方差进行假设检验（设显著性水平为0.05）。
均值假设检验:H0:μ=0;H1:μ≠0；x=[-6.3 -71.6 65.6 -79.2 -49.7 -81.9 74.9 -47.6 -120.8 56.9 100.9 47 9.7 -60.1 -52.7 86 80.6 -42.6 56.4 15.2];y=[28.9 1.6 61.7 -68 -41.3 -30.5 87 17.3 -17.8 1.2 -12.6 39.1 85 32.7 28.1 -9.3 -4.5 5.1 -32 -9.5];h1=ztest(x,0,70)h2=ztest(y,0,50)输出结果：h1 =0h2 =0方差假设检验:H0：σ2=σ02；H1：σ2≠σ02；x=[-6.3 -71.6 65.6 -79.2 -49.7 -81.9 74.9 -47.6 -120.8 56.9 100.9 47 9.7 -60.1 -52.7 86 80.6 -42.6 56.4 15.2];y=[28.9 1.6 61.7 -68 -41.3 -30.5 87 17.3 -17.8 1.2 -12.6 39.1 85 32.7 28.1 -9.3 -4.5 5.1 -32 -9.5]; function [h]=ktest(x,s0,alpha,tail)n=length(x);k=(n-1)*var(x)/(s0^2)%χ2分布检验方差if tail==0k1=chi2inv(alpha/2,n-1)k2=chi2inv(1-alpha/2,n-1)if k>=k1&k<=k2h=0;elseh=1;endendif tail==1k0=chi2inv(1-alpha,n-1)if k<=k0h=0;elseh=1;endendif tail==-1k0=chi2inv(alpha,n-1)if k>=k0h=0;elseh=1;endh1=ktest(x,70,0.05,0)h2=ktest(y,50,0.05,0)输出结果：h1 =0h2 =02)根据这组数据给出随机变量X和Y相关系数的一个点估计。
相关系数点估计：x=[-6.3 -71.6 65.6 -79.2 -49.7 -81.9 74.9 -47.6 -120.8 56.9 100.9 47 9.7 -60.1 -52.7 86 80.6 -42.6 56.4 15.2];y=[28.9 1.6 61.7 -68 -41.3 -30.5 87 17.3 -17.8 1.2-12.6 39.1 85 32.7 28.1 -9.3 -4.5 5.1 -32 -9.5]; r=corrcoef(x,y) 输出结果：r= 0.313412305102197 3) 用蒙特卡罗方法求炮弹落在椭圆形区域内的概率（取10000个数据点；请附程序）。
2222221(,)22(1)21x x y x x yx xy x p x y r r r σσσσπσσ⎧⎫⎡⎤⎪⎪=--+⎢⎥⎨⎬--⎢⎥⎪⎪⎣⎦⎩⎭%炮弹命中椭圆形区域概率源程序：a=110; b=90; sx=70; sy=50; r=0.3134123; z=0; n=10000;x=unifrnd(-a,a,1,n); y=unifrnd(-b,b,1,n); for i=1:nif (x(i)^2)/(a^2)+y(i)^2/(b^2)<=1u=exp(-0.5/(1-r^2)*(x(i)^2/sx^2-2*r*x(i)*y(i)/(sx*sy)+y(i)^2/sy^2)); z=z+u;endendP=4*a*b*z/(2*pi*sx*sy*sqrt(1-r^2)*n)输出结果：P =0.761272218724379考试课程数学实验2005.6.15下午班级学号姓名得分[说明]（1）第一、二、三题的答案直接填在试题纸上；（2）第四题将数学模型、简要解题过程和结果写在试题纸上；卷面空间不够时，请写在背面；（3）除非特别说明，所有计算结果小数点后保留4位数字。
