清华数学实验复习试题八(蒙特卡罗方法-龙格-库塔方法)
龙格-库塔实习题
欧拉法的程序：
function x=Euler(f,x0,y0,xn,N) x=zeros(1,N+1); y=zeros(1,N+1); y(1)=x0; y(1)=y0; h=(xn-x0)/N; for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(f,x(n),(yn)); end T=[x’,y ’] 欧拉法的思路： （1）给出自变量 x 的定义域[a,b],初始值 y0 及步长 h。 （2）对 k=0,1，...(b-a)/h，计算 y(k+1)=y(k)+hf(x(k),y(k)).
数学 1403 徐彤 学号：2014016544 二阶龙格-库塔方法及其实习题
二阶龙格-库塔方法的程序：
N=(b-a)/h; y=zeros(N+1,1); y(1)=y0; x=a:h:b; var=findsym(f); for i=2:N+1 v1= Funval(f,varvec,[x(i-1) y(i-1)]); t=y(i-1)+h*v1/2; V2=Funval(f,varvec,[x(i)+h/2 t]); Y(i)=y(i-1)+h*v2; end format short; function fv=Funval(f,varvec,varval) var=findsym(f); if length(var)<4 if var(1)==varvec(1) fv=subs(f,varvec(1),varvel(1)); else fv=subs(f,varvec(2),varvel(2)); end else fv=subs(f,var1403 徐彤 学号：2014016544 实习题：用欧拉法求一阶常微分方程
计算方法上机作业——龙格库塔法
计算方法上机作业——龙格库塔法龙格库塔法（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.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。
龙格库塔方法
• 可以通过检查步长折半前后两次计
算结果的偏差
=
(h)
y2 n1
y(h) n1
变步长方
来判断所选取的步长是否合法适
（1）对于给定的精度，如果 ,则反复将步长折半 进行计算直至 为止，这时取折半以后的“ 新值” 作为结果；
（2）如果 ，则反复将步长加倍，直至 为止， 这时取步长加倍前的“ 老值” 作为结果
四阶RungeKutta方法
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔 (h)和一个估算的斜率的乘积 决定。该斜率是以下斜率的加权平均：
k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率 k1 来 决定 y在点 tn + h/2的值； k3也是中点的斜率，但是这次采用斜率 k2决定 y值； k4是时间段终点的斜率，其 y值用 k3 决定。 当四个斜率取平均时，中点的斜率有更大的权值：
y=y0;z=zeros(c,6); %z生成c行，6列的零矩阵存放结 果；
不断迭代运算： for x=a:h:b
if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4);
五、变步长Runge-Kutta方
法• 从每一步看，步长越小，截断误差
越小；但随着步长的缩小，在一定 求解范围内所要完成的步数就会增 加，步数的增加不但引起计算量的 增大，而且可能导致舍入误差的严 重积累，因此需要选择步长
清华数学实验复习试题八(蒙特卡罗方法_龙格-库塔方法)
欢迎阅读考试课程 数学实验下午班级 姓名 学号 得分[说明]（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]搜索。
清华数学实验实验五蒙特卡罗方法
03 蒙特卡罗方法在清华数学 实验实验五中的应用
模拟随机过程
随机过程模拟
蒙特卡罗方法可以模拟各种随机 过程，如股票价格波动、气象变 化等，通过模拟这些过程，可以 更好地理解和预测实际现象。
概率分布模拟
蒙特卡罗方法可以生成符合特定 概率分布的随机数，用于模拟和 研究各种概率分布的性质和行为 。
求解数学问题
蒙特卡罗方法的优缺点
误差和不确定性
蒙特卡罗方法的精度取决于抽样次数，抽样次数越多，精 度越高，但计算成本也越高。同时，由于是随机模拟，结 果存在一定的不确定性。
对离散问题处理不佳
对于一些离散或非连续的问题，蒙特卡罗方法的精度可能 会受到影响。
对参数敏感
蒙特卡罗方法的参数选择对结果影响较大，需要谨慎选择。
02 清华数学实验实验五内容
实验目的
掌握蒙特卡罗方法的原理和应用。 学会使用蒙特卡罗方法解决实际问题。 培养数学建模和计算能力。
实验原理
蒙特卡罗方法是一种基于概率统 计的数值计算方法，通过随机抽
样和统计模拟来求解问题。
该方法适用于具有随机性和不确 定性的问题，通过大量模拟实验
来获得近似解。
蒙特卡罗方法的精度取决于模拟 实验的次数和随机抽样的质量。
金融工程
蒙特卡罗方法在金融工程中广泛应用于 风险评估、资产定价和衍生品定价等问
题。
工程设计
蒙特卡罗方法在工程设计中用于优化 设计参数、模拟系统性能和可靠性分
析等。
物理科学
在物理科学中，蒙特卡罗方法被用于 模拟分子运动、材料性质和量子力学 等领域。
社会科学
在社会科学中，蒙特卡罗方法被用于 模拟社会现象、预测人口变化和评估 政策效果等。
蒙特卡罗方法的优缺点
实验龙格—库塔方法
实验 龙格—库塔方法实验目的： 深入理解龙格—库塔方法的原理，同时学会用迭代方法看待某些特定问题。
比较本方法与其它方法（尤其是欧拉方法）解题的不同处，分析两者的精度。
实验要求: 编写程序，用龙格—库塔方法求解第三章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位数字)。
龙格库塔方法及其matlab实现
资料范本本资料为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正式推向市场。
龙格-库塔方法-文档资料
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
龙格-库塔方法(Runge-Kutta)
龙格-库塔方法〔Runge-Kutta〕3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即由前面(i-1)个已算出的表示，故公式是显式的.例)，而ki如当r=2时，公式可表示为(3.2.5)其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5)，要求选择参数，使公式的精度阶p 尽量高，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，由〔3.2.5〕式可知，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶精度R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)与中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数共有8个，所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法：(3.2.11)附：R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法与步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
龙格-库塔(Runge-Kutta)方法
——常微分方程及方程组的数值解 常微分方程及方程组的数值解
基本概念 只有一个自变量的微分方程为常微分方程 ODE, 否则称为偏微分方程 PDE。 。 方程中未知函数导数的最高阶数称为方程的阶
[1] [2] C. Runge, "Ueber die numerische Auflsung von Differentialgleichungen" Math. Ann. , 46 (1895) pp. 167–178 W. Kutta, "Beitrag zur naherungsweisen Integration von Differentialgleichungen" Z. Math. und Phys. , 46 (1901) pp. 435–453
2. 递推化 在具有唯一解的条件下， 在具有唯一解的条件下，通过步进法逐步计算出解 在一系列离散点上的值。从而得到原 的数值近似解。 在一系列离散点上的值。从而得到原ODE的数值近似解。 的数值近似解
初值问题解法 讨论一阶ODE（高阶可化为一阶ODEs） （高阶可化为一阶 讨论一阶 ） 的初值问题。初值问题可以一般地写成： 的初值问题。初值问题可以一般地写成：
欧拉（ 欧拉（Euler）方法 ）
Euler方法是求解初值问题的最简单方法， 方法是求解初值问题的最简单方法， 方法是求解初值问题的最简单方法 精度差。然而对理论分析很有用。 精度差。然而对理论分析很有用。Runge-Kutta 法是对Euler法的改进 法是对 法的改进
Euler方法： 方法： ym+1 = ym + hK1 + O(h ) K1 = f ( xm , ym )
K 1 = f ( xm , ym ) K i = f ( xm + hai , ym + h∑ bil K l )
龙格-库塔方法
h 按照数学学 习一般方法， yn 1 yn （k1 2k 2 2k3 +k 4） 6 继续上面的 做法，经过 k1 f ( xn , yn ) 较复杂的数 h 学推理和计 k （16） 2 f ( xn 1 , yn k1 ) 算，可以导 2 2 出四阶龙格 h －库塔格式， k f ( x 1 , yn k 2 ) 这里有一个 3 n 2 2 较有用的四 k f ( x , y hk ) 阶格式。 n 1 n 3 4
？
3、欧拉格式
yn 1 yn hf ( xn , yn ),
取点 xn 的斜率 k
*
(3)
n 0,1, 2,
f ( xn , yn )作为区间 [ xn , xn1 ] 上的平均斜率，精度低。
4、改进的欧拉格式
于是：
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn ))] 2 h yn 1 yn [k1 k2 ] 2
龙格－库塔方法
考虑一维经典初值问题
dy f ( x , y ) , y( x0 ) y0 , x [a, b] dx
(1)
龙格－库塔方法
一、设计思想
加权平均斜率。
1、微分中值定理
如果 f (x) 在 [a, b]上连续，在 ( a, b) 内可导，则在( a, b) 内至 少存在一点 ，有 f (b) f (a) f ' ( ) ba
xnq xn qh, 0 p q 1
令
yn1 yn h[(1 )k1 k2 k3 ] (14)
k1 , k 2
龙格库塔法
实验四 龙格库塔法【实验内容】1. 用标准四阶龙格－库塔方法设计算法2. 编程解微分方程初值问题；【实验方法与步骤】1． 龙格－库塔方法的基本思路设法计算(),f x y 在某些点上的函数值，然后对这些函数值做线性组合，构造近似计算公式；再把近似公式和解的泰勒展开式相比较，使前面的若干项吻合，从而达到较高的精度。
2． 四阶龙格－库塔方法的计算步骤求解'0()(,)()()y x f x y a x b y a y ⎧=≤≤⎨=⎩ 对上述给定的(,)f x y ，用四阶龙格-库塔法求解常微分方程初值问题112341213243(22)6(,)11(,)2211(,)22(,)n n n n n n n n n n h y y k k k k k hf x y k hf x h y k k hf x h y k k hf x h y k +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ 程序如下#include<stdio.h>void Runge_Kutta(float( * f)(float x,float y),float a,float b,float y0,int N){ float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;printf("x[0]=%f\ty[0]=%f\n",x,y);for(i=1;i<=N;i++){ K1=( * f)(x,y);K2=( * f)(x+h/2,y+h*K1/2);K3=( * f)(x+h/2,y+h*K2/2);K4=( * f)(x+h,y+h*K3);y=y+h*(K1+2*K2+2*K3+K4)/6;x=a+i*h;printf("x[%d]=%f\ty[%d]=%f\n",i,x,i,y); }}float f(float x,float y){ return y-2*x/y;}void main(){float a=0,b=1,y0=1;Runge_Kutta(f,a,b,y0,5);}[实验内容]课本Ｐ１８２页例７．３．２【实验结果】。
第3章03龙格-库塔方法
这样有计算公式：
yn1 yn h[(1 )k1 k2 )
k1 f (xn , yn ), k2 f (xn ph, yn phk1)
y(h) n1
1 16
y( xn1)
(h)
y2 n1
1 15
(h)
y2 n1
y(h) n1
误差事后估计
可以通过检查步长折半前后两次计算结果的偏差
y y ( h ) 2
(h)
n1
n1
来判断所选取的步长是否合适.
1) 对于给定的精度ε,如果δ> ε,则反复将步长折半 进行计算,直到δ< ε为止,这时取步长折半后的 新值作为结果;
k4
k3 f ( xn 2h, yn 21hk1 22hk2 ) f ( xn 3h, yn 31hk1 32hk2 33hk3)
四阶龙格－库塔公式每步要四次计算函数值，具有四阶精 度，局部截断误差是O(h5) .
四阶龙格－库塔格式
下列经典格式是其中常用的一种:
yn1 k2
yk
1.000000 1.005000 1.019025 1.041218 1.070800 1.107076 1.149404 1.197210 1.249975 1.307228 1.368514
y( xk ) yk
0.0 1.6×10－4 2.9×10－4 4.0×10－4 4.8×10－4 5.5×10－4 5.9×10－4 6.2×10－4 6.5×10－4 6.6×10－4 6.6×10－4
调用 2 次函数。有 2 阶精度。
龙格库塔法-原理及程序实现
龙格库塔法一、基本原理：“龙格-库塔(Runge-Kutta)方法”是一种在工程上应用广泛的“高精度单步算法”。
由于此算法精度高，采取措施对误差进行抑制，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即：yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)通常所说的龙格-库塔法就是指四阶——龙格库塔法，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。
（1）计算公式（1）的局部截断误差是。
龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。
二、小程序#include<stdio.h>#include<math.h>#define f(x,y) (-1*(x)*(y)*(y))void main(void){double a,b,x0,y0,k1,k2,k3,k4,h;int n,i;printf("input a,b,x0,y0,n:");scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n); printf("x0\ty0\tk1\tk2\tk3\tk4\n");for(h=(b-a)/n,i=0;i!=n;i++){k1=f(x0,y0);k2=f(x0+h/2,y0+k1*h/2);k3=f(x0+h/2,y0+k2*h/2);k4=f(x0+h,y0+h*k3);printf("%lf\t%lf\t",x0,y0);printf("%lf\t%lf\t",k1,k2);printf("%lf\t%lf\n",k3,k4);y0+=h*(k1+2*k2+2*k3+k4)/6;x0+=h;}printf("xn=%lf\tyn=%lf\n",x0,y0);}运行结果：input a,b,x0,y0,n:0 5 0 2 20x0y0k1k2k3 k40.000000 2.000000-0.000000-0.500000-0.469238-0.8861310.250000 1.882308-0.885771-1.176945-1.129082-1.2800600.500000 1.599896-1.279834-1.295851-1.292250-1.2227280.750000 1.279948-1.228700-1.110102-1.139515-0.9901621.000000 1.000027-1.000054-0.861368-0.895837-0.7528521.2500000.780556-0.761584-0.645858-0.673410-0.5621891.5000000.615459-0.568185-0.481668-0.500993-0.4205371.7500000.492374-0.424257-0.361915-0.374868-0.3178552.0000000.400054-0.320087-0.275466-0.284067-0.2435982.2500000.329940-0.244935-0.212786-0.218538-0.1894822.5000000.275895-0.190295-0.166841-0.170744-0.1495632.7500000.233602-0.150068-0.132704-0.135399-0.1197033.0000000.200020-0.120024-0.106973-0.108868-0.0970483.2500000.172989-0.097256-0.087300-0.088657-0.0796183.5000000.150956-0.079757-0.072054-0.073042-0.0660303.7500000.132790-0.066124-0.060087-0.060818-0.0553054.0000000.117655-0.055371-0.050580-0.051129-0.0467434.2500000.104924-0.046789-0.042945-0.043363-0.0398334.5000000.094123-0.039866-0.036750-0.037072-0.0342024.7500000.084885-0.034226-0.031675-0.031926-0.029571xn=5.000000yn=0.076927。
Matlab中龙格-库塔（Runge-Kutta）方法原理及实现
函数功能ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。
ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)³。
解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.使用方法[T,Y] = ode45(odefun,tspan,y0)odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf]y0 是初始值向量T 返回列向量的时间点Y 返回对应T的求解列向量[T,Y] = ode45(odefun,tspan,y0,options)options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)在设置了事件参数后的对应输出TE 事件发生时间YE 事件解决时间IE 事件消失时间sol =ode45(odefun,[t0 tf],y0...)sol 结构体输出结果应用举例1 求解一阶常微分方程程序:一阶常微分方程odefun=@(t,y) (y+3*t)/t^2; %定义函数tspan=[1 4]; %求解区间y0=-2; %初值[t,y]=ode45(odefun,tspan,y0);plot(t,y) %作图title('t^2y''=y+3t,y(1)=-2,1<t<4')legend('t^2y''=y+3t')xlabel('t')ylabel('y')% 精确解% dsolve('t^2*Dy=y+3*t','y(1)=-2')% ans =一阶求解结果图% (3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)2 求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]程序:function Testode45tspan=[3.9 4.0]; %求解区间y0=[2 8]; %初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y'' ''=-t*y + e^t*y'' +3sin2t')xlabel('t')ylabel('y')function y=odefun(t,x)y=zeros(2,1); % 列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend高阶求解结果图相关函数ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tbMatlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
龙格库塔法
第九章 常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时nab h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法1.欧拉法 1．对常微分方程初始问题(9.2))((9.1)),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
杨环宇小组试卷
一、填空题(共4小题12分，每小题3分)1、初值问题在数值求解过程中的误差有两个来源：和1、MATLAB都是当作处理。
如果是矩阵，MATLAB会根据输入的数值自动扩充以满足需要。
3、将变量范围为无限的连续的周期函数变换成无限的离散的傅里叶频谱序列。
4、一维矢量的梯度就是用计算导数，在端点直接用作为近似值。
二、单选题(共4小题12分，每小题3分)1、在MTLAB中，泰勒展开的计算指令是（）A: diff B：taylor C: int D: limit2、表示插值函数可以视为两个插值基函数的线性组合，其组合系数是对应点上的函数值。
这种形式的插值称为（）。
A.牛顿插值B.埃尔米特插值C.拉格朗日插值D.不确定3、下列选项中哪一个选项代表当前工作目录（）mand WindowB.Current DirectorymandHistoryD.Workspace4、二阶常微分方程y″=f(x,y,y′),a≤x≤b的边值问题中，第三类边值问题的条件是（）。
A、y′(a)-α0y(a)=α1y′(b)+β0y(b)=β1,B、y(a)=α,y(b)=βC、y′(a)+α0y(a)=α1y′(b)-β0y(b)=β1,D、y′(a)=α,y′(b)=β三、判断题，对的打√，错的打×（共4小题12分，每小题3分）1、变量名用字母打头，后面可以跟字母、数字、下划线，数目不限。
但只有变量名的前面64个符号有限。
（）2、龙格-库塔法的基本思想是为了提高精度，可多取几点的斜线值作加权平均当做平均斜率k ave.（）3、由于端点的导数值并不是已知的数据，所以就可以对它进行选择（）4、洛巴托法在数值计算中也称为高斯积分法（）简答题(共5小题32分)1、数据的预处理有哪些方法？2、对于变步长的龙格-库塔法，可以通过检查步长折半后前后两次计算结果的偏差⊿=|y i+1(h/2)- y i+1(j)|来判断所选取的步长是否合适，具体可分为哪两种情况？3、数值计算难免产生误差，请说出产生误差的原因。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
清华数学实验复习试题八(蒙特卡罗方法-龙格-库塔方法)考试课程数学实验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位数字。