数模第一章答案(张绍辉给).doc
- 文档编号:18694902
- 上传时间:2023-09-17
- 格式:DOC
- 页数:22
- 大小:368KB
数模第一章答案(张绍辉给).doc
《数模第一章答案(张绍辉给).doc》由会员分享,可在线阅读,更多相关《数模第一章答案(张绍辉给).doc(22页珍藏版)》请在冰点文库上搜索。
第一章习题参考答案
1.请编写绘制以下图形的MATLAB命令,并展示绘得的图形.
(1)、分别是椭圆的内切圆和外切圆.
解答
方法一(显函数和伸缩变换)
的显函数形式为,并利用伸缩变换:
的横、纵坐标都是的两倍,的横、纵坐标分别是的两倍和一倍.
编写程序时运用好MATLAB函数plot的语法格式2(x是向量,y是矩阵),以及格式4,使程序简洁.
使用命令axisequal,才能绘得真正的圆.
程序:
x=-1:
.05:
1;%由40段折线连接成半圆周
y=sqrt(1-x.^2);
plot(x,[y;-y],'k',2.*x,[y;-y;2.*y;-2.*y],'k')
axisequal
title('方法一(显函数)')
绘得的图形:
评价:
方法一绘得的图形在外切圆和椭圆的左右两端看起来明显还是折线,而在其余地方看起来比较光滑,原因在外切圆和椭圆的左右两端,导数趋于无穷大,所以,虽然x的步长是固定的,但是在左右两端,y会比别处有更显著的变化.当然,如果令x的步长更小,例如x=-1:
.01:
1,绘得的图形将会看起来更光滑一些.
方法二(参数方程和伸缩变换)
的参数方程为,关于伸缩变化和MATLAB函数plot的语法的讨论与方法一相同.特意选取参数t的步长,使得半圆周仍然由40段折线连接而成,如同方法一一样.
程序:
t=linspace(0,2*pi,81);%由40段折线连接成半圆周
x=cos(t);
y=sin(t);
plot(x,y,'k',2.*x,[y;2.*y],'k')
axisequal
title('方法二(参数方程)')
绘得的图形:
评价:
虽然半圆周由同样多的折线段连接而成,但是方法二绘得的图形看起来处处一样光滑,事实上,方法二通过等分圆心角来取得圆周上的采样点,并连结线段,所以绘得的“圆形”实际上是正多边形.
(2)指数函数和对数函数的图像关于直线y=x对称.
解答
指数函数和对数函数互为反函数,在MATLAB函数plot的输入当中交换x和y的次序,就实现反函数,而且图像就是关于直线y=x对称.
根据两点确定一条直线的原理,绘制直线段只需给出两端点的坐标.
使用命令axisequal,才能绘得真正的对称图形,加上坐标网格,能增强对称效果.
指数函数的自变量x的取值区间的左端不能太小,否则绘得的图像会在左边有一段与x轴重合.
程序:
x=-3:
.1:
3;
y=exp(x);
plot(x,y,'k',y,x,'k',[-3,20],[-3,20],'k')
axisequal
axis([-3,20,-3,20])
grid
xlabel('x')
ylabel('y')
title('y=e^x和y=lnx的函数图像关于直线y=x对称')
绘得的图形:
(3)黎曼函数
的图像(要求分母q的最大值由键盘输入).
解答
输入的英文单词是input,通过在MATLAB帮助文档检索input这个关键词,查到实现键盘输入的MATLAB函数是input,语法格式为:
(1)user_entry=input('prompt')输入项是一个字符数组,input将该字符数组显示在命令窗口,作为提示语,等待用户从键盘输入数值数组并按回车键;input将用户从键盘输入的数值数组赋值给user_entry所代表的变量名.
(2)user_entry=input('prompt','s')第二输入项是规定的字符“'s'”,表示等待用户从键盘输入的是字符数组.
条件“是既约分数,0 当q=2时,有黎曼函数的图像的一个坐标点(1/2,1/2); 当q=3时,有黎曼函数的图像的两个坐标点(1/3,1/3)和(2/3,1/3); 当q=4时,有黎曼函数的图像的两个坐标点(1/4,1/4)和(3/4,1/4),而x=2/4不是既约分数;…… 可见,在程序中应设置二重循环语句: 第一重(外重)是分母q从2到由键盘输入的最大值的循环,步长为1; 第二重(内重)是分子p从1到q-1的循环,步长为1,如果p和q的最大公约数等于1,就把p/q和1/q分别添加入横坐标数组x和纵坐标数组y. 注意x和y需要在循环语句之前说明为空数组. 当x为有理数(0 程序: n=input('分母的最大值n='); x=[]; y=[]; forq=2: n forp=1: q-1 ifgcd(p,q)==1 x=[x,p./q]; y=[y,1./q]; end end end plot(x,y,'k.') axis([0,1,0,1]) xlabel('x') ylabel('y') s=strcat('分母的最大值n=',num2str(n)); gtext(s); title('第一章习题1(3),黎曼函数的图像') 当分母的最大值n=30时绘得的图形: 当分母的最大值n=100时绘得的图形: 3.两个人玩双骰子游戏,一个人掷骰子,另一个人打赌掷骰子者不能掷出所需点数,输赢的规则如下: 如果第一次掷出3或11点,打赌者赢;如果第一次掷出2、7或12点,打赌者输;如果第一次掷出4,5,6,8,9或10点,记住这个点数,继续掷骰子,如果不能在掷出7点之前再次掷出该点数,则打赌者赢.请模拟双骰子游戏,要求写出算法和程序,估计打赌者赢的概率.你能从理论上计算出打赌者赢的精确概率吗? 请问随着试验次数的增加,这些概率收敛吗? 解答 (一)算法 输入模拟试验的次数n; 输出打赌者赢的概率p. 第1步初始化计数器k=0; 第2步对i=1,2,…,n,循环进行第3~7步; 第3步产生两个在1~6这6个整数中机会均等地取值的随机数,并把这两个随机数之和赋值给x; 第4步如果x是3或11,那么k加1,进入下一步循环;否则,做第5步; 第5步如果x不是2、7和12,那么做第6~8步;否则,直接进入下一步循环; 第6步产生两个在1~6这6个整数中机会均等地取值的随机数,并把这两个随机数之和赋值给y; 第7步如果y不等于x,也不等于7,重复第6步所做的; 第8步如果y等于7,那么k加1,进入下一步循环;否则,直接进入下一步循环; 第9步计算概率p=k./n. (二)MATLAB程序 n=10000; k=0; fori=1: n x=ceil(rand*6)+ceil(rand*6); ifx==3|x==11 k=k+1; elseifx~=2&x~=7&x~=12 y=ceil(rand*6)+ceil(rand*6); whiley~=x&y~=7 y=ceil(rand*6)+ceil(rand*6); end ify==7 k=k+1; end end end p=k/n 注&与&&,|与||分别是短运算符与长运算符,功能有区别,但在上面的程序中效果一样。 (三)理论计算 每个骰子的点数有6种结果: 1,2,3,4,5,6.掷两个骰子,掷出的点数组合共有6×6=36种结果.定义是点数之和的概率,则 ,,, ,,. 玩双骰子游戏,第1次掷出的点数之和为3或11,则打赌者赢,可算得此情形打赌者赢的的概率为 当第1次掷出的点数之和是4,5,6,8,9或10,继续掷骰子,直到掷出的点数之和是7或原来的值为止,如果先得到的点数之和是7,则打赌者赢,利用条件概率,可算得此情形打赌者赢的概率为 所以打赌者赢的理论概率P为 (四)收敛性分析 一次打赌相当于伯努利概型,记为随机变量X,取值为0(表示打赌者输)或1(表示打赌者赢),则X的期望为P,方差为. n次打赌,即相互独立地重复试验n次,试验结果可记作随机变量序列,则打赌者赢的频率为平均值. 弱大数定律: ,都有 . 强大数定律: . 中心极限定理: 当时,随机变量 的分布趋向于标准正态分布(也就是说,当n充分大的时候,随机变量的分布近似于均值为P、方差为的正态分布). 用循环语句实现以下计算: 考虑试验次数n=100、400、1600和6400共四种情况(即),在每一种情况下,重复计算1000次,得到频率的1000个值,用MATLAB函数std计算无偏标准差,用MATLAB函数hist绘制直方图,添加标示理论概率P的铅直线,用asix控制坐标范围. 程序: n=[100,400,1600,6400]; form=1: 4; forj=1: 1000 k=0; fori=1: n(m) x=ceil(rand*6)+ceil(rand*6); ifx==3|x==11 k=k+1; elseifx~=2&x~=7&x~=12 y=ceil(rand*6)+ceil(rand*6); whiley~=x&y~=7 y=ceil(rand*6)+ceil(rand*6); end ify==7 k=k+1; end end end p(j,m)=k/N(m); end end P=251/495%打赌者赢精确概率,极限分布的均值 S=sqrt(P*(1-P)./n)%极限分布的均方差 sigma=std(p)%1000个频率的无偏标准差 form=1: 4 subplot(2,2,m) hist(p(: m),15) holdon plot([P,P],[0,250],'k') axis([min(p(: 1)),max(p(: 1)),0,250])%范围固定 str=strcat('试验次数n=',num2str(n(m))); title(str) holdoff end 命令窗口显示的计算结果为: P= 0.50707 S= 0.0499950.0249980.0124990.0062494 sigma= 0.0517770.0254460.0121440.0063688 绘得的图形(四幅直方图的横坐标范围是固定的): 从计算结果和绘得的图形可以看出: 固定试验次数n,重复计算1000次,得到频率的1000个值,的确如中心极限定理所论断的,只要n足够大,这1000个值的分布就近似于均值为P、方差为的正态分布(可以做正态总体的假设检验,详见课本6.1节); 随着试验次数n的增加,频率的1000个值的分布范围越来越窄,这是因为无偏标准差总是近似于,并随n的增加而越来越小. 这就说明: 随着试验次数n的增加,模拟得到的概率收敛于理论概率P=0.50707. 4.根据表1.14的数据,完成下列数据拟合问题: 表1.14美国人口统计数据(百万人) 年份 1790 1800 1810 1820 1830 1840 1850 1860 人口 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 年份 1870 1880 1890 1900 1910 1920 1930 1940 人口 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 年份 1950 1960 1970 1980 1990 2000 人口 150.7 179.3 204.0 226.5 251.4 281.4 (1)如果用指数增长模型模拟美国人口1790年至2000年的变化过程,请用MATLAB统计工具箱的函数nlinfit计算指数增长模型的以下三个数据拟合问题: (i)取定=3.9,=1790,拟合待定参数; (ii)取定=1790,拟合待定参数和; (iii)拟合待定参数、和. 要求写出程序,给出拟合参数和误差平方和的计算结果,并展示误差平方和最小的拟合效果图. 解答 程序如下: t=1790: 10: 2000; x=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92,... 106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4]; ea=@(b,t)3.9.*exp(b (1).*(t-1790)); ba=nlinfit(t,x,ea,0.03) sea=sum((x-ea(ba,t)).^2) eb=@(b,t)b (2).*exp(b (1).*(t-1790)); bb=nlinfit(t,x,eb,[0.03,3.9]) seb=sum((x-eb(bb,t)).^2) ec=@(b,t)b (2).*exp(b (1).*(t-b(3))); bc=nlinfit(t,x,ec,[0.03,3.9,1790]) sec=sum((x-ec(bc,t)).^2) [bb (1)-bc (1),seb-sec,eb(bb,1790)-ec(bc,1790)] 计算结果为: ba= 0.021194 sea= 17418 bb= 0.01422314.994 seb= 2263.9 Warning: TheJacobianatthesolutionisill-conditioned,andsomemodelparametersmaynotbeestimatedwell(theyarenotidentifiable).Usecautioninmakingpredictions. >Innlinfitat224 bc= 0.0142237.75071743.6 sec= 2263.9 ans= -1.3146e-0084.1567e-0083.5492e-005 结果分析: (i)取定=3.9,=1790,拟合,则得=0.021194,误差平方和等于17418,最大; (ii)取定=1790,拟合和,则得=14.994,=0.014223,误差平方和等于2263.9; (iii)若拟合、和,则得=1743.6,=7.7507,=0.014223,误差平方和等于2263.9,、误差平方和以及的计算结果与方法b的误差都非常微小,但是MATLAB给出警告信息,指出条件存在病态,参数未必能拟合得好. 综上所述,(ii)是最佳拟合方案,其拟合效果图为: 绘图程序: plot(t,x,'k+',t,eb(bb,t),'ko') axis([1780,2010,0,300]) xlabel('年份t'),ylabel('美国人口x(百万)') title('指数增长模型非线性拟合效果图') legend('实际值','理论值',2) 注参考课本73页注1.7.4,如果比较(ii)和(iii)用nlinfit算得的协方差矩阵估计COV,会发现(ii)的COV为: COV= 2.8016e-007-0.00075625 -0.000756252.1129 而(iii)的COV为: COV= 2.9668e-007-8.3191e+005-7.5469e+006 -8.3191e+0053.8741e+0203.5145e+021 -7.5469e+0063.5145e+0213.1883e+022 前者是准确的,而后者却存在很大的误差; 事实上,指数增长模型是微分方程初值问题的经过点的积分曲线的显函数表达式.(i)要拟合的积分曲线的初值点固定在点(1790,3.9),只有待定;(ii)要拟合的积分曲线的初值点(1790,)和参数都待定.(ii)的候选的积分曲线的集合比(i)的大得多,所以(ii)能找到误差平方和更小的结果. 至于(ii)和(iii),在理论上会找到同一条积分曲线,但是在计算上,由于(ii)的待定参数比(iii)的少一个,计算也更精确、更省时;方法c的待定参数和之中有一个冗余,会导致严重的计算误差. (2)通过变量替换,可以将属于非线性模型的指数增长模型转化成线性模型,并用MATLAB函数polyfit进行计算,请说明转化成线性模型的详细过程,然后写出程序,给出拟合参数和误差平方和的计算结果,并展示拟合效果图. 解答 对指数增长模型两边求对数,得 固定,引进变量替换 ,,, 则转化为一次多项式.根据表1.14的数据,计算得,其中 , 然后用MATLAB函数polyfit计算出的最小二乘估计值,进而 程序: t=1790: 10: 2000; x=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92,... 106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4]; X=t-1790; Y=log(x); b=polyfit(X,Y,1) r=b (1) x_0=exp(b (2)) eb=@(a,t)a (2).*exp(a (1).*(t-1790)); xhat=eb([r,x_0],t); sse=sum((x-xhat).^2) plot(t,x,'k+',t,xhat,'ko') axis([1780,2010,0,300]) xlabel('年份t'),ylabel('美国人口x(百万)') title('指数增长模型线性化拟合效果图') legend('实际值','理论值',2) 计算结果为: b= 0.0202191.7992 r= 0.020219 x_0= 6.045 sse= 34892 拟合得=6.045,=0.020219,误差平方和等于34892.拟合效果图为: (3)请分析指数增长模型非线性拟合和线性化拟合的结果有何区别? 原因是什么? 解答 指数增长模型线性化拟合的误差平方和比非线性拟合大很多;观察拟合误差比较图,可以发现: 非线性拟合的拟合误差比较均匀,线性化拟合的拟合误差却随着人口的增加而越开越大. 绘图程序: plot(x,x-eb(bb,t),'kx',x,x-xhat,'ks',[0,300],[0,0],'k') xlabel('美国人口(百万)'),ylabel('拟合误差') title('指数增长模型的拟合误差比较图') legend('非线性拟合','线性化拟合',3) 原因: 因为对于数值越大的数据,由于求对数带来的损失越大,以至于线性化拟合的误差越大. (4)如果用阻滞增长模型模拟美国人口1790年至2000年的变化过程,请用MATLAB统计工具箱的函数nlinfit计算阻滞增长模型的以下三个数据拟合问题: (i)取定=3.9,=1790,拟合待定参数和; (ii)取定=1790,拟合待定参数、和; (iii)拟合待定参数、、和. 要求写出程序,给出拟合参数和误差平方和的计算结果,并展示误差平方和最小的拟合效果图. 解答 程序如下: t=1790: 10: 2000; x=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76,92,... 106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4]; za=@(b,t)b (1).*3.9./(3.9+(b (1)-3.9).*exp(-b (2).*(t-1790))); ba=nlinfit(t,x,za,[300,.02]) sea=sum((x-za(ba,t)).^2) zb=@(b,t)b (1).*b(3)./(b(3)+(b (1)-b(3)).*exp(-b (2).*(t-1790))); bb=nlinfit(t,x,zb,[300,.02,3.9]) seb=sum((x-zb(bb,t)).^2) zc=@(b,t)b (1).*b(3)./(b(3)+(b (1)-b(3)).*exp(-b (2).*(t-b(4)))); bc=nlinfit(t,x,zc,[300,.02,3.9,1790]) sec=sum((x-zc(bc,t)).^2) 计算结果为: ba= 342.440.027353 sea= 1224.9 bb= 446.570.0215477.6981 seb= 457.74 Warning: TheJacobianatthesolutionisill-conditioned,andsomemodelparametersmaynotbeestimatedwell(theyarenotidentifiable).Usecautioninmakingpredictions. >Innlinfitat224 bc= 446.570.0215475.17521771.3 sec= 457.74 结果分析: (i)若取定=3.9,=1790,拟合和,则得=0.027353,=342.44,误差平方和等于1224.9; (ii)若取定=1790,拟合、和,则得=7.6981,=0.021547,=446.57,误差平方和
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数模 第一章 答案 张绍辉