清华数学实验复习试题八(蒙特卡罗方法-龙格-库塔方法)
清华数学实验复习试题八(蒙特卡罗方法-龙格-库塔方法)考试课程数学实验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位数字。
龙格-库塔实习题
欧拉法的程序：
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正式推向市场。