清华数学实验复习试题八蒙特卡罗方法龙格库塔方法
龙格-库塔方法（Runge-Kutta）
龙格-库塔⽅法（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值(即)，⽽k由前⾯(i-1)个已算出的表⽰，故公式是显式的.例i如当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个，所应满⾜的⽅程为这是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 。
82第二节 龙格—库塔法
k
(1)
h h k 若令 yn1 y xn hy xn y xn y xn (2) 2! k! 则 y xn1 yn1 O hk 1
y0 k1 2 k2 hf x0 h 2, y0 k1 2
y0 k3
k4 hf x0 h, y0 k3
y0 k2 2 k3 hf x0 h 2, y0 k2 2
k
x1 x0 h y1 y0 k
数学学院 信息与计算科学系
0.1832292
0.1584376
数学学院 信息与计算科学系
接上图
0.4 0.5 0.5 0.6 0.6 1.341667 1.416026 1.412676 1.482627 1.483281 0.0745394 0.0710094 0.0708400 0.0673253
0.1416245
数学学院 信息与计算科学系
数学学院 信息与计算科学系
由表8-4可见，虽然四阶龙格-库塔方法每步要 计算四次 f 的值，但以h=0.2为步长ቤተ መጻሕፍቲ ባይዱ计算结果就
有5 位有效数字，而欧拉法与预估计-校正方法以
h=0.1为步长的计算结果才具有2 位与3 位有效数字.
如果步长 h 也取0.2，则结果的精度会更低.
即公式(2)为k 阶方法.
数学学院 信息与计算科学系
二、龙格-库塔方法(R-K方法)
R-K方法不是通过求导数的方法构造近似公式, 而是通过计算不同点上的函数值, 并对这些函数值作 线性组合, 构造近似公式, 再把近似公式与解的泰勒 展开式进行比较, 使前面的若干项相同 , 从而使近似 公式达到一定的阶数.
龙格库塔
数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
则，对于该问题的RK4由如下方程给出：其中这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：∙k1是时间段开始时的斜率；∙k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点t n + h/2的值；∙k3也是中点的斜率，但是这次采用斜率k2决定y值；∙k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值：RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
显式龙格库塔法显示龙格－库塔法是上述RK4法的一个推广。
它由下式给出其中如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(h p+1)时的条件。
这些可以从舍入误差本身的定义中导出。
例如，一个2阶精度的2段方法要求b1 + b2 = 1, b2c2 = 1/2, 以及b2a21 = 1/2。
在Matlab下输入：edit，然后将下面两行百分号之间的内容，复制进去，保存%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%function dxdt=ode_Miss_ghost(t,x)%分别用x(1),x(2),x(3),x(4)代替N1,P1,N2,P2N1=x(1);P1=x(2);N2=x(3);P2=x(4);K=2;tau_c=3e-9;tan_p=6e-12;beta =5e-5;delta=0.692;eta =0.0001;fm =8e6;Ith =26e-3;Ib =1.5*Ith;Im =0.3*Ith;I1=Ib+Im*sin(2*pi*fm*t)+K*P2;I2=Ib+Im*sin(2*pi*fm*t)+K*P1;dxdt=[(I1/Ith-N1-(N1-delta)/(1-delta)*P1)/tau_e;((N1-delta)/(1-delta)*(1-eta*P1)*P1-P1+beta*N1)/tau_p;(I2/Ith-N2-(N2-delta)/(1-delta)*P2)/tau_e;((N2-delta)/(1-delta)*(1-eta*P2)*P2-P2+beta*N2)/tau_p;]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在Matlab下面输入：t_start=0;t_end=2e-9;y0=[1e-3;1e-4;0;0]; %初值[x,y]=ode15s('ode_Miss_ghost',[0,t_end],y0);plot(x,y);legend('N1','P1','N2','P2');xlabel('x');Mathlab定步长龙格－库塔法[345阶]、自适应步长rkf45经常看到很多朋友问定步长的龙格库塔法设置问题，下面吧定步长三阶、四阶、五阶龙格库塔程序贴出来，有需要的可以看看ODE3 三阶龙格－库塔法CODE:function Y = ode3(odefun,tspan,y0,varargin)%ODE3 Solve differential equations with a non-adaptive method of order 3.% Y = ODE3(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates% the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE3(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN% but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Bogacki-Shampine Runge-Kutta method of order 3.%% Example% tspan = 0:0.1:20;% y = ode3(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1, % and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,3);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:});Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3));endY = Y.';ODE4 四阶龙格－库塔法CODE:function Y = ode4(odefun,tspan,y0,varargin)%ODE4 Solve differential equations with a non-adaptive method of order 4.% Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the classical Runge-Kutta method of order 4.%% Example% tspan = 0:0.1:20;% y = ode4(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,4);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:}); F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:}); F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));endY = Y.';定步长RK4(自编):CODE:function Y=RungeKutta4(f,xn,y0)% xn=0:.1:1;% y0=1;% y_n=[];% f=@(X1,Y1) Y1-2*X1/Y1;y_n=[];h=diff(xn(1:2));for i=1:length(xn)-1K1=f(xn(i),y0);K2=f(xn(i)+h/2,y0+h*K1/2);K3=f(xn(i)+h/2,y0+h*K2/2);K4=f(xn(i)+h,y0+h*K3);y_n=[y_n;y0+h/6*(K1+2*K2+2*K3+K4)];y0=y_n(end);endY=y_n;ODE5 五阶龙格－库塔法CODE:function Y = ode5(odefun,tspan,y0,varargin)%ODE5 Solve differential equations with a non-adaptive method of order 5.% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Dormand-Prince method of order 5 in a general % framework of explicit Runge-Kutta methods.%% Example% tspan = 0:0.1:20;% y = ode5(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);% Method coefficients -- Butcher's tableau%% C | A% --+---% | BC = [1/5; 3/10; 4/5; 8/9; 1];A = [ 1/5, 0, 0, 0, 03/40, 9/40, 0, 0, 044/45 -56/15, 32/9, 0, 019372/6561, -25360/2187, 64448/6561, -212/729, 0 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]; % More convenient storageA = A.';B = B(:);nstages = length(B);F = zeros(neq,nstages);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);% General explicit Runge-Kutta frameworkF(:,1) = feval(odefun,ti,yi,varargin{:});for stage = 2:nstageststage = ti + C(stage-1)*hi;ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));F(:,stage) = feval(odefun,tstage,ystage,varargin{:});endY(:,i) = yi + F*(hi*B);endY = Y.';-------------------------------------------------------------------------------------------------------------------自适应步长RKF45(相当于ode45)ODE45 是4阶方法提供候选解，5阶方法控制误差。
龙格库塔法求解方程
例题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]搜索。
龙格-库塔方法-文档资料
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 阶精度，只需
龙格库塔法介绍
h 0.005稳定.
2) 改进欧拉法(预测 — 校正，即二阶R K法)：
yn
1
yn
h
k1 2
k2 2
,
k1
f (xn, yn ) yn,
k2 f (xn h, yn hk1) ( yn hyn )，
即yn1 故
yn
当初值准确即e0 0时，整体误差为en O(h p ).
证明
考察单步法的收敛性归结为验证增量函数(x, y,h)是否
满足Lipschitz条件.
欧拉法 : (x, y,h) f (x, y)，L L.
改进的欧拉法：
yn1
yn
h[ 2
f
(xn, yn )
f
( xn 1,
xn
f
( x,
y ( x))dx
r
h ci
i 1
f
( xn
ih,
y ( xn
ih)).
yn1 yn h(xn, yn, h),
(3.4)
其中
r
(xn, yn, h) ciki ,
(3.5)
i 1
k1 f (xn, yn )，
欧拉法r 1, p 1.改进 欧拉法r 2, p 2.
k1)
k3）
k3 f (xn h, yn hk1 2hk2 )
称为库塔三阶方法.
阶数p和段数r（计算函数值次数）的关系
r12 p1 2
3 4 5 6 7 r≥8 3 4 4 5 6 r-2
常用的经典四阶龙格 库塔方法:
龙格-库塔(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( k1 + k2 ) yn +1 = yn + 2 hk1 1 1 3 3 k1 = f [ xn + ( + )h, yn + )hk2 ] +( + 2 6 4 4 6 hk2 1 3 1 3 k2 = f [ xn + ( − )h, yn + ( − )hk1 + ] 2 6 4 6 4
1 库塔三阶方法 库塔三阶方法
h h k1 = f ( xn , yn ); k2 = f ( xn + , yn + k1 ) 2 2
h yn +1 = yn + ( k1 + 4k2 + k3 ) 6
k3 = f ( xn + h, yn − hk1 + 2hk2 )
四级方法： 四级方法：N = 4 常见的2种四阶方法： 常见的 种四阶方法： 种四阶方法 经典龙格 库塔方法 龙格1经典龙格-库塔方法
1 中点法（修正的Euler法） 取 c1 = 0, c2 = 1, a2 = 1 中点法（修正的 法 2 h h
y n +1 = y n + hf ( x n +
常见的3种二级方法： 常见的 种二级方法： 种二级方法
二阶龙格库塔方法 2 二阶龙格库塔方法
h yn +1 = yn + [ f ( xn , yn ) + f ( xn + h, yn + hf ( xn , yn ))] 2
1 取 c1 = c2 = , a2 = 1 2
2
, yn +
2
f ( x n , y n ))
类似于N 的推导方法， 的推导方法 三级方法： 三级方法：N = 3 类似于 = 2的推导方法，可得到 1 c1 + c 2 + c 3 = 1; c 2 a 2 + c 3 a 3 = ; 2 O ( h4 ) 1 1 2 2 c 2 a 2 + c 3 a 3 = ; c 3 a 2 b32 = 6 3 常见的2种三阶方法： 常见的 种三阶方法： 种三阶方法
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
考试课程数学实验下午班级姓名学号得分[说明]（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 命令是______ 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)]输出结果：三、 (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__。
线性代数方程组解的误差分析： 故其误差上限为：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) 输出结果:A=[5 -7 0 1 ;-3 22 6 2 ;5 -1 31 -1 ;2 1 0 23];D=diag(diag(A)); %从稀疏矩阵A 中提取D L=-tril(A,-1); %从稀疏矩阵A 中提取LU=-triu(A,1); %从稀疏矩阵A 中提取U b=[6 3 4 7]'; %设定方程组右端项向量x= zeros(4,1); %设定方程组初始向量 m= inv(D-L)*U;n= inv(D-L)*b; %高斯-赛德尔迭代法 for j2=1:5 y=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.当瞄准目标的中心发射炮弹时，在众多随机因素的影响下，弹着点服从以目标中心为均值的二维正态分布，设弹着均值假设检验: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=3)用蒙特卡罗方法求炮弹落在椭圆形区域内的概率（取10000个数据点；请附程序）。
%炮弹命中椭圆形区域概率源程序：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 =考试课程数学实验下午班级学号姓名得分[说明]（1）第一、二、三题的答案直接填在试题纸上；（2）第四题将数学模型、简要解题过程和结果写在试题纸上；卷面空间不够时，请写在背面；（3）除非特别说明，所有计算结果小数点后保留4位数字。
（4）考试时间为120分钟。
一、（10分）某厂生产A、B两种产品，1千克原料在甲类设备上用12小时可生产3件A，可获净利润64元；在乙类设备上用8小时可生产4件B，可获净利润56元。
该厂每天可获得55千克原料，每天总的劳动时间为480小时，且甲类设备每天至多能生产80件A。
试为该厂制订生产计划使每天的净利润最大。
1）以生产A、B产品所用原料的数量x1、x2（千克）作为决策变量，建立的数学规划模型是：2）每月的最大净利润是_____________元。
若要求工人加班以增加劳动时间，则加班费最多为每小时__________元。
若A获利增加到27元/件，应否改变生产计划？_____________二、(10分) 已知常微分方程组初值问题试用数值方法求()8y π=_______________________(保留小数点后5位数字)。
你用的MATLAB 命令是__________________________________________，其精度为_____________。
三、(10分) 已知线性代数方程组Ax =b , 其中5701322625131121023A -⎡⎤⎢⎥-⎢⎥=⎢⎥--⎢⎥⎣⎦, ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=4321x x x x x , 8325b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,若方程组右端项有小扰动]1.0,0,0,0['=b δ，试根据误差估计式估计11x x δ≤_____________ （x x δ,分别表示原问题的解和右端项小扰动后对应的解的变化量）；若取初值]0,0,0,0[)0('=x ，则用高斯-赛德尔迭代法求解Ax =b 时，=)5(x ____________________________________________；对本题而言，此迭代方法是否收敛_________________，原因是__________________________________。