试题八(蒙特卡罗方法_龙格-库塔方法)

考试课程数学实验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=[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.00000000400584844.999999993870908z=-3.069999999925403e+003ans=33.9999999989193572.5000000001404410.0000000002784052）每天的最大净利润是___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)x=9.99999999999940045.000000000000625z=-3.209999999999987e+003二、(10分)已知常微分方程组初值问题221'''()0,4xyxyxy()2,2y2'().2y试用数值方法求()6y__1.73205____(保留小数点后5位数字)。你用的MATLAB命令是______ode45(@f,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)]输出结果：0.5235987755982991.732050795523993三、(10分)已知线性代数方程组Ax=b,其中5701322625131121023A,4321xxxxx,6347b,若方程组右端项有小扰动]1.0,0,0,0[b，试根据误差估计式估计11xx___0.0743___（xx,分别表示原问题的解和右端项小扰动后对应的解的变化量）；若取初值]0,0,0,0[)0(x，则用高斯-赛德尔迭代法求解Ax=b时，)5(x_(1.7160,0.3926,-0.1306,0.1381)_；对本题而言，此迭代方法是否收敛___是__，原因是__谱半径ρ(B)=0.3971__。线性代数方程组解的误差分析：1()xbbAAcondAxbb故其误差上限为：A=[5-701;-32262;5-131-1;21023];b=[6347];db=[0000.1];d=cond(A,1)*norm(db,1)/norm(b,1)输出结果:d=0.074339065208930A=[5-701;-32262;5-131-1;21023];D=diag(diag(A));%从稀疏矩阵A中提取DL=-tril(A,-1);%从稀疏矩阵A中提取LU=-triu(A,1);%从稀疏矩阵A中提取Ub=[6347]';%设定方程组右端项向量x=zeros(4,1);%设定方程组初始向量m=inv(D-L)*U;n=inv(D-L)*b;%高斯-赛德尔迭代法forj2=1:5y=m*(x(:,j2));fori=1:4x(i,j2+1)=y(i,:)+n(i,:);endendt2=x(:,end)%输出迭代法最终结果j2输出结果:t2=1.7159723472264450.392646824062879-0.1305711006230470.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.665.6-79.2-49.7-81.974.6-47.6-120.856.9Y28.91.661.7-68-41.3-30.58717.3-17.81.2X100.9479.7-60.1-52.78680.6-42.656.415.2Y-12.639.18532.728.1-9.3-4.55.1-32-9.51)根据这组数据对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=0h2=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==0k1=chi2inv(alpha/2,n-1)k2=chi2inv(1-alpha/2,n-1)ifk=k1&k=k2h=0;elseh=1;endendiftail==1k0=chi2inv(1-alpha,n-1)ifk=k0h=0;elseh=1;endendiftail==-1k0=chi2inv(alpha,n-1)ifk=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.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=0.3134123051021973)用蒙特卡罗方法求炮弹落在椭圆形区域内的概率（取10000个数据点；请附程序）。22222211(,)exp22(1)21xxyxxyxxyxpxyrrr%炮弹命中椭圆形区域概率源程序：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: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位数字。（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分)已知常微分方程组初值问题2212'''()0,()2,'()422xyxyxyyy试用数值方法求()8y_______________________(保留小数点后5位数字)。你用的MATLAB命令是__________________________________________，其精度为_____________。三、(10分)已知线性代数方程组Ax=b,其中5701322625131121023A,4321xxxxx,8325b,若方程组右端项有小扰动]1.0,0,0,0[b，试根据误差估计式估计11xx_____________（xx,分别表示原问题的解和右端项小扰动后对应的解的变化量）；若取初值]0,0,0,0[)0(x，则用高斯-赛德尔迭代法求解Ax=b时，)5(x____________________________________________；对本题而言，此迭代方法是否收敛_________________，原因是__________________________________。四、（20分）炮弹射击的目标为一椭圆形区域，在X方向半轴长100m，Y方向半轴长80m.当瞄准目标的中心发射炮弹时，在众多随机因素的影响下，弹着点服从以目标中心为均值的二