清华数学实验复习试题八蒙特卡罗方法龙格库塔方法.docx
《清华数学实验复习试题八蒙特卡罗方法龙格库塔方法.docx》由会员分享，可在线阅读，更多相关《清华数学实验复习试题八蒙特卡罗方法龙格库塔方法.docx（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≤55
12*x1+8*x2≤480
3*x1≤80
x1,x2
0
基本模型：
max（y）=64*x1+54*x2
s.t.x1+x2≤55
12*x1+8*x2≤480
3*x1≤80
x1,x2
0
c=[6454];
A1=[11;
128;
30];
b1=[55;480;80];
v1=[00];
[x,z,ef,out,lag]=linprog（-c,A1,b1,[],[],v1）
lag.ineqlin
输出结果:
x=10.000000004005848
ans=
0.000000000278405
2）每天的最大净利润是___3070__元。
若要求工人加班以增加劳动时间，则加班费最多
为每小时__2.5__元。
若A获利增加到26元/件，应否改变生产计划？
____不变___
c=[7854];
A1=[11;
128;
30];
b1=[55;480;80];
v1=[00];
[x,z,ef,out,lag]=linprog（-c,A1,b1,[],[],v1）
45.000000000000625
二、（10分）已知常微分方程组初值问题
试用数值方法求
__1.73205____（保留小数点后5位数字）。
你用的MATLAB命令是______ode45（@ff,ts,y0）______，其精度为____四阶__。
%待解常微分方程组函数M文件源程序：
functiondy=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,其中
若方程组右端项有小扰动
，试根据误差估计式估计
___0.0743___（
分别表示原问题的解和右端项小扰动后对应的解的变化量）；若取初值
，则用高斯-赛德尔迭代法求解Ax=b时，
_（1.7160,0.3926,-0.1306,0.1381）_；对本题而言，此迭代方法是否收敛___是__，原因是__谱半径ρ（B）=0.397<1__。
线性代数方程组解的误差分析：
故其误差上限为：
A=[5-701;-32262;5-131-1;21023];
b=[6347];
db=[0000.1];
d=cond（A,1）*norm（db,1）/norm（b,1）
输出结果:
A=[5-701;-32262;5-131-1;21023];
D=diag（diag（A））;%从稀疏矩阵A中提取D
L=-tril（A,-1）;%从稀疏矩阵A中提取L
U=-triu（A,1）;%从稀疏矩阵A中提取U
b=[6347]';%设定方程组右端项向量
x=zeros（4,1）;%设定方程组初始向量
m=inv（D-L）*U;
n=inv（D-L）*b;%高斯-赛德尔迭代法
forj2=1:
5
y=m*（x（:
j2））;
fori=1:
4
x（i,j2+1）=y（i,:
）+n（i,:
）;
end
end
t2=x（:
end）%输出迭代法最终结果
j2
输出结果:
t2=
判敛:
lamda=eig（inv（D-L）*U）
pubanjing=max（abs（lamda））
输出结果:
四、（20分）炮弹射击的目标为一椭圆形区域，在X方向半轴长110m，Y方向半轴长90m.当瞄准目标的中心发射炮弹时，在众多随机因素的影响下，弹着点服从以目标中心为均值的二维正态分布，设弹着点偏差的均方差在X方向和Y方向分别为70m和50m。
今测得一组弹着点的横纵坐标如下：
X
-6.3
-71.6
65.6
-79.2
-49.7
-81.9
74.6
-47.6
-120.8
56.9
Y
28.9
1.6
61.7
-68
-41.3
-30.5
87
17.3
-17.8
1.2
X
100.9
47
9.7
-60.1
-52.7
86
80.6
-42.6
56.4
15.2
Y
-12.6
39.1
85
32.7
28.1
-9.3
-4.5
5.1
-32
-9.5
1）根据这组数据对X方向和Y方向的均值和均方差进行假设检验（设显著性水平为0.05）。
均值假设检验:
H0:
μ=0;H1:
μ≠0；
x=[-6.3-71.665.6-79.2-49.7-81.974.9-47.6-120.856.9100.9479.7-60.1-52.78680.6-42.656.415.2];
y=[28.91.661.7-68-41.3-30.58717.3-17.81.2-12.639.18532.728.1-9.3-4.55.1-32-9.5];
h1=ztest（x,0,70）
h2=ztest（y,0,50）
输出结果：
h1=0
h2=0
方差假设检验:
H0：
σ2=σ02；H1：
σ2≠σ02；
x=[-6.3-71.665.6-79.2-49.7-81.974.9-47.6-120.856.9100.9479.7-60.1-52.78680.6-42.656.415.2];
y=[28.91.661.7-68-41.3-30.58717.3-17.81.2-12.639.18532.728.1-9.3-4.55.1-32-9.5];
function[h]=ktest（x,s0,alpha,tail）
n=length（x）;
k=（n-1）*var（x）/（s0^2）%χ2分布检验方差
iftail==0
k1=chi2inv（alpha/2,n-1）
k2=chi2inv（1-alpha/2,n-1）
ifk>=k1&k<=k2
h=0;
else
h=1;
end
end
iftail==1
k0=chi2inv（1-alpha,n-1）
ifk<=k0
h=0;
else
h=1;
end
end
iftail==-1
k0=chi2inv（alpha,n-1）
ifk>=k0
h=0;
else
h=1;
end
h1=ktest（x,70,0.05,0）
h2=ktest（y,50,0.05,0）
输出结果：
h1=0
h2=0
2）根据这组数据给出随机变量X和Y相关系数的一个点估计。
相关系数点估计：
x=[-6.3-71.665.6-79.2-49.7-81.974.9-47.6-120.856.9100.9479.7-60.1-52.78680.6-42.656.415.2];
y=[28.91.661.7-68-41.3-30.58717.3-17.81.2-12.639.18532.728.1-9.3-4.55.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）;
fori=1:
n
if（x（i）^2）/（a^2）+y（i）^2/（b^2）<=1
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;
end
end
P=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分）已知常微分方程组初值问题
试用数值方法求
_______________________（保留小数点后5位数字）。
你用的MATLAB命令是__________________________________________，其精度为_____________。
三、（10分）已知线性代数方程组Ax=b,其中
若方程组右端项有小扰动
，试根据误差估计式估计
_____________（
分别表示原问题的解和右端项小扰动后对应的解的变化量）；若取初值
，则用高斯-赛德尔迭代法求解Ax=b时，
____________________________________________；对本题而言，此迭代方法是否收敛_________________，原因是__________________________________。
四、（20分）炮弹射击的目标为一椭圆形区域，在X方向半轴长100m，Y方向半轴长80m.当瞄准目标的中心发射炮弹时，在众多随机因素的影响下，弹着点服从以目标中心为均值的二维正态分布，设弹着点偏差的均方差在X方向和Y方向分别为70m和50m。
今测得一组弹着点的横纵坐标如下：
X
-6.3
-71.6
65.6
-79.2
-49.7
-81.9
74.6
-47.6
-100.8
56.9
Y
28.9
1.6
61.7
-68
-41.3
-30.5
87
17.3
-17.8
1.2
X
100.9
47
9.7
-60.1
-52.7
86
80.6
-42.6
56.4
15.2
Y
-12.6
39.1
85
32.7
28.1
-9.3
-4.5
5.1
-32
-9.5
1）根据这组数据对X方向和Y方向的均值和均方差进行假设检验（设显著性水平为0.05）。
2）根据这组数据给出随机变量X和Y相关系数的一个点估计。
3）用蒙特卡罗方法求炮弹落在椭圆形区域内的概率（取10000个数据点；请附程序）。
考试课程数学实验参考答案与评分标准
A卷（班级-姓名-学号-得分）
一、1）
（如果进一步要求3x1和4x2为非负整数，不扣分）
2）3070元，2.5元；不变
二、1.73203（或1.73205）,ode23（或ode45），3级2阶（或5级4阶）
（不写3级（或5级）不扣分；个别同学可能用其他命令，则结果相应略有变化）
三、0.0743,[1.7160,0.3926,-0.1306,0.1381]’,收敛,谱半径为0.3968<1
（不写出谱半径的具体数值不扣分，但写错要扣分）
四、1）对均值做的假设为
（X,Y方向相同）X,Y方向均接受H0
对X方向的方差做的假设为
H0:
=4900,H0:
4900（如果做单侧检验，可不扣分）接受H0
对Y方向的方差做的假设为
H0:
=2500,H0:
2500（如果做单侧检验，可不扣分）接受H0
2）相关系数的点估计为0.313（用r=corrcoef（x,y）命令）
3）大约0.76，结果具有随机性
[附]主要程序示例：
%1）~2）
x=[-6.3-71.665.6-79.2-49.7-81.974.6-47.6-120.8（%B-100.8）56.9100.9479.7-60.1-52.78680.6-42.656.415.2];
h1=ztest（x,0,70）,%x方向均值检验
y=[28.91.661.7-68-41.3-30.58717.3-17.81.2-12.639.18532.728.1-9.3-4.55.1-32-9.5];
h2=ztest（y,0,50）,%Y方向均值检验
r=corrcoef（x,y）%相关系数的点估计
pause
n=20;
alpha=0.05;
sx2=var（x）,sx0=70;
chi2=（n-1）*sx2/（sx0^2）
chi2alpha=chi2inv（1-alpha,n-1）
ifchi2<=chi2alphaH0=0
elseH0=1
end
sy2=var（y）,sy0=50;
chi2=（n-1）*sy2/（sy0^2）
chi2alpha=chi2inv（1-alpha,n-1）
ifchi2<=chi2alphaH0=0
elseH0=1
end
%3）
a=0.7;b=0.5;m=0;z=0;
p=0.313;c=1.1;d=0.9;%A
%p=0.311;c=1;d=0.8;%B
n=10000;
fori=1:
n
x=2*rand（1,2）-1;
y=0;
ifx
（1）^2+x
（2）^2<=1
y=exp（-0.5/（1-p*p）*（c^2*x
（1）^2/a^2+d^2*x
（2）^2/b^2-2*p*c*d*x
（1）*x
（2）/a/b））;
z=z+y;
m=m+1;
end
end
P=4*c*d*z/2/pi/a/b/sqrt（1-p*p）/n,m
B卷（班级-学号-姓名）
一、1）
2）3160元；2元；不变
二、1.53077（或1.53073）,龙格-库塔方法ode23（或ode45），3级2阶（或5级4阶）
三、0.0826,[2.3416,0.5359,-0.2961,-0.0095]’,收敛,谱半径为0.3968<1
四、1）同A卷。
2）0.311;3）大约0.69，结果具有随机性
评分标准：
一、1）5分（每个式子一分），2）前两空每空2分，最后一空1分。
二、第一空6分，后两空每空2分。
三、第一空3分，中间一空5分，最后两空每空1分。
四、1）8分（4个检验每个2分）；2）2分；
3）10分（看程序及结果分析：
积分表达式2分，程序5分，结果3分）
考试课程数学实验参考答案与评分标准
A卷（班级-姓名-学号-得分）
一、1）
（如果进一步要求3x1和4x2为非负整数，不扣分）
2）3070元，2.5元；不变
二、1.73203（或1.73205）,ode23（或ode45），3级2阶（或5级4阶）
（不写3级（或5级）不扣分；个别同学可能用其他命令，则结果相应略有变化）
三、0.0743,[1.7160,0.3926,-0.1306,0.1381]’,收敛,谱半径为0.3968<1
（不写出谱半径的具体数值不扣分，但写错要扣分）
四、1）对均值做的假设为
（X,Y方向相同）X,Y方向均接受H0
对X方向的方差做的假设为
H0:
=4900,H0:
4900（如果做单侧检验，可不扣分）接受H0
对Y方向的方差做的假设为
H0:
=2500,H0:
2500（如果做单侧检验，可不扣分）接受H0
2）相关系数的点估计为0.313（用r=corrcoef（x,y）命令）
3）大约0.76，结果具有随机性
[附]主要程序示例：
%1）~2）
x=[-6.3-71.665.6-79.2-49.7-81.974.6-47.6-120.8（%B-100.8）56.9100.9479.7-60.1-52.78680.6-42.656.415.2];
h1=ztest（x,0,70）,%x方向均值检验
y=[28.91.661.7-68-41.3-30.58717.3-17.81.2-12.639.18532.728.1-9.3-4.55.1-32-9.5];
h2=ztest（y,0,50）,%Y方向均值检验
r=corrcoef（x,y）%相关系数的点估计
pause
n=20;
alpha=0.05;
sx2=var（x）,sx0=70;
chi2=（n-1）*sx2/（sx0^2）
chi2alpha=chi2inv（1-alpha,n-1）
ifchi2<=chi2alphaH0=0
elseH0=1
end
sy2=var（y）,sy0=50;
chi2=（n-1）*sy2/（sy0^2）
chi2alpha=chi2inv（1-alpha,n-1）
ifchi2<=chi2alphaH0=0
elseH0=1
end
%3）
a=0.7;b=0.5;m=0;z=0;
p=0.313;c=1.1;d=0.9;%A
%p=0.311;c=1;d=0.8;%B
n=10000;
fori=1:
n
x=2*rand（1,2）-1;
y=0;
ifx
（1）^2+x
（2）^2<=1
y=exp（-0.5/（1-p*p）*（c^2*x
（1）^2/a^2+d^2*x
（2）^2/b^2-2*p*c*d*x
（1）*x
（2）/a/b））;
z=z+y;
m=m+1;
end
end
P=4*c*d*z/2/pi/a/b/sqrt（1-p*p）/n,m
B卷（班级-学号-姓名）
一、1）
2）3160元；2元；不变
二、1.53077（或1.53073）,龙格-库塔方法ode23（或ode45），3级2阶（或5级4阶）
三、0.0826,[2.3416,0.5359,-0.2961,-0.0095]’,收敛,谱半径为0.3968<1
四、1）同A卷。
2）0.311;3）大约0.69，结果具有随机性
评分标准：
一、1）5分（每个式子一分），2）前两空每空2分，最后一空1分。
二、第一空6分，后两空每空2分。
三、第一空3分，中间一空5分，最后两空每空1分。
四、1）8分（4个检验每个2分）；2）2分；
3）10分（看程序及结果分析：
积分表达式2分，程序5分，结果3分）
升级会员 