2015 高教社杯全国大学生数学建模竞赛 C 题论文.pdf

2015 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 C 题论文题论文 月上柳梢头，人约黄昏后月上柳梢头，人约黄昏后 摘摘 要要 “月上柳梢头，人约黄昏后”里面所提到的约会时间，就是月出与黄昏后同 时出现的时间段。本文通过建立数学模型，并以北京为例,计算了北京的各参数 值，与现实数据作比较来验证模型，再判断什么条件下会出现“月上柳梢头，人 约黄昏后”这一现象，从而给出黄昏后的定义以及发生这一情景的条件。然后根 据条件，利用 Excel 表格对哈尔滨、上海、广州、昆明、成都、和乌鲁木齐这六 个城市的经纬度、日落时间、月出时间、以及日落月出时间差等进行了计算和统 计，通过分析比较，从而判断出各个城市地区能否出现“月上柳梢头，人约黄昏 后”这一现象。 关键词关键词：地心天顶距，日落，月出，月亮高度，黄昏后 问题重述问题重述 “月上柳梢头，人约黄昏后”是北宋学者欧阳修的名句，写的是与佳人相约 的情景。请用天文学的观点赏析该名句，并进行如下的讨论： 1. 定义“月上柳梢头”时月亮在空中的角度和什么时间称为“黄昏后” 。根据天 文学的基本知识，在适当简化的基础上，建立数学模型，分别确定“月上柳 梢头”和“人约黄昏后”发生的日期与时间。并根据已有的天文资料（如太 阳和月亮在天空中的位置、日出日没时刻、月出月没时刻）验证所建模型的 合理性。 2. 根据所建立的模型，分析 2016 年北京地区“月上柳梢头，人约黄昏后”发生 的日期与时间。根据模型判断 2016 年在哈尔滨、上海、广州、昆明、成都、 乌鲁木齐是否能发生这一情景？如果能，请给出相应的日期与时间；如果不 能，请给出原因。 问题分析问题分析 “月上柳梢头，人约黄昏后”是北宋学者欧阳修的名句，释义：去年元宵夜 之时，花市上灯光明亮如同白昼。佳人相约，在月上柳梢头之时，黄昏之后的情 景。显然出现这一情景的时间比较短，否则佳人相约就没有了确切的时间，必会 造成某一方长时间的空等。 经过人们的长期观察， 出现这一现象， 必定是在太阳、 月亮和地球三者共线之日，也就是只会发生在农历每月的十五日。 然后根据这一现象以及我们的生活中的实际观察，做出适当的假设，经过计 算，定义出“月上柳梢头”时月亮在空中的角度范围。 再通过分析与计算，得出某地农历每月的日落与月出时间，看其交叉范围是 否在定义的角度范围内，从而确定“月上柳梢头” 、 “人约黄昏后”这一现象发生 的时间条件。再根据日落与月出时间的比较与分析，确定黄昏后的时间。 通过建立数学模型，并对模型进行验证分析和修正，再根据现象发生的时间 条件进而算出 2016 年北京地区“月上柳梢头，人约黄昏后”发生的日期与时间。 通过对哈尔滨、上海、广州、昆明、成都、乌鲁木齐这六个城市的经纬度， 日落时间，月出时间的统计，再计算出日落月出的时间差以及月亮与地平面的夹 角，从而判定这些城市是否会发生“月上柳梢头，人约黄昏后”的现象。 模型假设模型假设 1. 假设柳树高度为 5m，人距柳树的距离 15 米,人的身高为 1.6m,根据三角 函数和相似三角形基本数学知识求出月亮在空中的角度为 12.77。 2.假设当时诗人是在现在的北京，假设当时的月亮与地平面的夹角是 0 20。 3.假设没有雾霾、台风以及各种天气因素的影响。 4.假设把观测点当作一个理想的点来验算。 5.假设云层对太阳光没有散射效应。 符号说明符号说明 h 为太阳平黄经 Ps 为太阳近地点的平黄经 为赤道夹角(黄赤交角) T 表示儒略世纪数 p 为月亮在近地点的平黄经 N 为月亮升交点的平黄经 为观测点的地方恒星时 模型模型一一的建立与求解的建立与求解 我们观测到的太阳高度和方位角与太阳在天球中的坐标和观测点所处的位 置有关。太阳在天球中的坐标是时间的函数。我们首先介绍计算太阳在天球中位 置的基本单元随时间变化的天文参数， 然后介绍随天文参数变化的太阳坐标 以及太阳相对于观测点的地心天顶距的计算方法。 根据这些基本参数给出日落时 刻、太阳高度的计算方法。 天文参数天文参数 1895 年纽康姆给出某时刻太阳天文参数h，Ps、的计算公式。如下： 2 279.6966836000.768920.00030hTT 23 281.220831.1719020.004500.000003PsTTT 23 23 278.261 46.8450.00590.00183TTT 其中，T表示儒略世纪数，由儒略日数计算，其计算公式为： 2451545 36525 JD T 其中，JD为儒略日数，为自 1900 年 1 月 0 日 12 时起至计算时刻之间的天 数。可从天文年历中查出，本文运用下列公式计算: 设Y为给定年份，M为月份，D为该月日期（可以带小数） 对格里高利历，有A=INT() 100 Y ,2() 4 A BAINT, 则 (365.25*(Y 4716)INT(30.6001*(M 1)D B 1524.52JDINT 其中， 第一项为所考虑年份， 每年按 365 天计算， 第二项为考虑闰年的情况， 第三项考虑 1900 年 1 月 0 日的 12 个小时。 太阳的坐标太阳的坐标 已知上述天文参数，可求出计算时刻太阳的黄经 ： 0.033501*sin(h)0.000351*sin(2(h) S hPsPs 0.000005*sin(3*(h)Ps 因为太阳在黄道上，所以黄纬为零。 由于天赤道、黄道和太阳所在的时圈构成直角三角形；其中直角就是天赤道 和时圈的交角。根据球面正弦定律，这个三角形的边和角有如下关系： sin:sinsin :sin90 s 式中，为太阳赤纬；为按（3）式计算黄赤交角。上式可写为： sinsin *sin s （1） （2） （3） （4） （6） （7） （5） （8） 太阳相对于观测点的位置太阳相对于观测点的位置 本研究采用观测点的地心天顶距来表示太阳在某时刻相对于某一点的位 置，设观测点的经纬度、，则可由下式计算某时的地心天顶距： cossin*sin *sincos*cos*cossin *cos *sin sss 其中为观测点的地方恒星时，由下式计算： (14) 15180小算时刻（小时） 太阳太阳落下时刻的计算落下时刻的计算 根据(9)式可求出某时刻太阳相对于某一点的地心天顶距。若地心天顶距为 90 ，则该点在晨昏线上。因此，求解方程 即可得到儒略世纪数T，然后转化为某天的时分秒，即可得到某天的太阳落 山时间。 求解日落时间求解日落时间 经查询可知，北京的地理位置是东经 116 度 24 分，北纬 39 度 54 分，时区 是东八区，将北京时间 2015 年 11 月 7 日代入模型计算，见附录一，可以算的北 京的日落时间为：17.3572165434923 时，即 17:21，与我们实际查询到的当日北 京的日落时间 17:06 基本吻合， 见附录二， 故该模型求解日落时间是比较精确的。 假设作者欧阳修的创作时间是 1999 年，经网上资料查询这首诗是在元夜所 作，即农历的正月十五左右，即 3 月 1 日，假设创作地点是北京。 将数据代入该模型，见附录三，计算可得到： 日落时间为：17.0406013165627 时，即：17：02 （9） （10） （11） 模型模型二二的建立与求解的建立与求解 已知确定月亮位置的太阳参数表示如下： 2 279.6966836000.768920.00030hTT 23 23 278.261 46.8450.00590.00183TTT 根据布朗在 1919 年给出计算月亮位置的天文参数s， p ，N，其中s为月 亮的平黄经，其角速度为每小时 0.5490165 ， p 为月亮在近地点的平黄经，其 角速度为每小时 0.0046418 ，N为月亮升交点的平黄经，其角速度为每小时 0.0000020 。其计算公式如下： 23 270.69668481267.890570.001980.000002sTTT 23 334.329564069.34030.010320.00001pTTT 23 259.183281934.142010.002080.00002NTTT 月亮的坐标月亮的坐标 由以上的天文参数，便可求解出时刻月球的黄经 m 和黄纬 m ： 0.089504*sin(s N)0.004897*sin(2*s p)0.004847*sin(p) m NN 0.003024*sin(2*)shN 0.109760*sin(s p)0.022236*sin(22*h p)0.011490*sin(2*(s h) m s 0.003728*sin(2*s 2*p) 月球相对于观测点的位置月球相对于观测点的位置 用观测点的地心天顶距来表示月亮在某时刻相对于某一点的位置， 设观测 点的经纬度为和，则可由下式计算某时某地的地心天顶距： cossin *sin *cos*sincos *sin mmm （12） （13） （14） （15） （16） （17） （18） cos *cos*cos*cos * sin*cos*cossin *sin mmmmm 其中为观测点的地方恒星时，可由以下式子来计算： 11 *15180小算时刻 小时 月出时刻的计算月出时刻的计算 根据（19）式可求出某时刻月亮相对于某一点的地心天顶距。若地心天顶距 为90，则该点为月出或月落的时刻。因此，只要求解以下方程即可： 1sin *sin *cos*sincos *sin mmm cos *cos*cos*cos * sin*cos*cossin *sin mmmmm 根据所得到的儒略世纪数，然后转化为某天的时分秒，即可得到某天的月球 出升时间和落山时间。 模型求解模型求解 同样把时间 1999 年 3 月 1 日，地点是北京（东经 116 度 24 分，北纬 39 度 54 分）代入模型，见附录四，经计算求得月出时间为：16.4513879494116 时， 即 16:27。 显然，在某一地点某一个人的视角范围内，月亮从升起到落下是 180 度，用 时 12 小时，则此时在该点的水平面上，月亮的角速度： 180 15 / 12 wh h 由此可见 “月上柳梢头， 人约黄昏后” 发生这一情景， 月亮与水平线的夹角： *w日落时间-月出时间 代入可知是 8.838 。刚好符合月上柳梢头的 0 20 ，所以这个时间地 点假设成立。 经查询资料知，黄昏时候，月球高度十度左右，因此可知以上计算结果的正 确性。并且，资料显示，黄昏后有确切的定义，即太阳入地十度。所以根据地心 天顶距，代入模型一，可以算出 1999 年 3 月 1 日北京黄昏后的时间为： 17.3588536507063 时，即 17:21.并将其代入（22） 、 （23）式，算的月亮与水平 （19） （20） （21） （22） （23） 线的夹角为：13.6，这与假设一基本吻合，所以该定义是十分准确的，即黄昏 后为太阳入地十度。 现设Rl为日落时间，Yc为月出时间，则 *wRlYc 。又因为出现“月 上柳梢头，人约黄昏后”这一现象，月亮与水平线的夹角得在 0 20 之间，即 0*20wRlYc，也就是： 20 0, 15 RlYc 。 针对问题二，要计算 2016 年北京地区“月上柳梢头，人约黄昏后”发生的 日期与时间，以及判断 2016 年在哈尔滨、上海、广州、昆明、成都、乌鲁木齐 是否能发生这一情景？只需要算出各地区农历每月十五日的日落月出时间， 并检 验是否满足（24）式，倘若满足，则会发生这一情景，否则，将不会出现“月上 柳梢头，人约黄昏后”这一现象。 模型的验算及预测模型的验算及预测 运用模型可以计算出 2016 年北京地区农历每月十五的日落月出时间（见附 录五） ，并将日落与月出的时间差随月份的变化可做出以下折线图： 北京 -3.99 -2.66 -1.33 0 1.33 2.66 123456789101112 月份 日落-日出 由图像可以看出， 只有在 6 月 19 日和 7 月 18 日北京日落月出的时间差才满 足（24）式，即为 2016 年北京地区“月上柳梢头，人约黄昏后”发生的日期与 时间。 同时，我们还可以通过模型，计算出 2016 年哈尔滨、上海、广州、昆明、 成都、乌鲁木齐各地区的日落时间、月出时间、以及它们的时间差（见附录六） ， （24） 并将其均做成了以下折线图： 哈尔滨 -5.32 -3.99 -2.66 -1.33 0 1.33 2.66 12345678910 11 12 月份 日落-月出 上海 -3.99 -2.66 -1.33 0 1.33 123456789101112 月份 日落-月出 广州 -2.66 -1.33 0 1.33 2.66 123456789101112 月份 日落-月出 昆明 -3.99 -2.66 -1.33 0 1.33 12345678910 11 12 月份 日落-月出 成都 -3.99 -2.66 -1.33 0 1.33 12345678910 11 12 月份 日落-月出 乌鲁木齐 -5.32 -3.99 -2.66 -1.33 0 1.33 123456789101112 月份 日落-月出 根据以上的折线图与（24）式，可以预测出 2016 年哈尔滨、上海、广州、 昆明、成都、乌鲁木齐这六个地区出现“月上柳梢头，人约黄昏后”这一情景的 日期。其中，上海、昆明以及成都这三个地区出现这情景的日期均为 7 月 18 日， 广州则会在 4 月 21 日和 9 月 15 日发生这个情景，乌鲁木齐出现的日期为 7 月 18 日、8 月 17 日和 9 月 15 日，而在哈尔滨出现这个场景的次数最多，分别是在 4 月 21 日、5 月 21 日、6 月 19 日、8 月 17 日和 9 月 15 日。 模型评价模型评价 本文在正确、清晰地分析了题意地基础上，建立了合理科学的数学模型。这 两个模型的思路都比较清晰，一目了然，无论是模型的表达，还是后面的求解过 程都比较明了，每个表达式的含义都很明确，而且易于应用于现实的生活中。 不足的是，模型过于理想化，考虑的影响因素较少，忽略了各种天气问题， 以及云层对太阳光的散射作用等，所以求解的结果可能存在一些误差。 参考文献参考文献 1万永革、庄献华，防专校园太阳高度及升降方位的计算J，防灾技术高 等专科学校学报，2003 年 3 月 2万永革、孟晓春、黄猛，月亮高度及升降时刻与方位的计算，防灾技术 高等专科学校学报，2003 年 9 月 3 附录附录 附录一：附录一： JD=365.25*(2015+4716)+367+7+(2-20+5)-1524.52 JD = 2457334.23 T=(JD-2451545)/36525 T = 0.158500479123887 h=4.881627973+628.3319509*T+0.000005235988*T2 h = 104.472543371037 Ps=4.908229409+0.020453548*T+0.000078539816*T2+0.00000005235988*T3 Ps = 4.9114732794751 e=0.409319678+0.00022712*T-0.000000028604*T2+0.00000000887209*T3 e = 0.409355675945545 r=h+0.033501*sin(h-Ps)+0.000351*sin(2*(h-Ps)+0.000005*sin(3*(h-Ps) r = 104.444582343072 x=solve(sin(39.95)*sin(e)*sin(r)+cos(39.95)*cos(r)*x+sqrt(1-x2)*cos(e)*sin(r)=0,x) x = (2.5*(100000000000000000000000000000000.0*cos(e)*sin(r)*(3.9547046624944954147383792 959034e65*cos(r)2 + 1.0e66*cos(e)2*sin(r)2 - 6.0452953375055045852616207040966e65*sin(e)2*sin(r)2)(1/2) + 4.8895150738482495508826868606319e64*cos(r)*sin(e)*sin(r)/(9.8867616562362385368459 482397586e64*cos(r)2 + 2.5e65*cos(e)2*sin(r)2) -(2.5*(100000000000000000000000000000000.0*cos(e)*sin(r)*(3.954704662494495414738379 2959034e65*cos(r)2 + 1.0e66*cos(e)2*sin(r)2 - 6.0452953375055045852616207040966e65*sin(e)2*sin(r)2)(1/2) - 4.8895150738482495508826868606319e64*cos(r)*sin(e)*sin(r)/(9.8867616562362385368459 482397586e64*cos(r)2 + 2.5e65*cos(e)2*sin(r)2) x2= -(2.5*(100000000000000000000000000000000.0*cos(e)*sin(r)*(3.954704662494495414738379 2959034e65*cos(r)2 + 1.0e66*cos(e)2*sin(r)2 - 6.0452953375055045852616207040966e65*sin(e)2*sin(r)2)(1/2) - 4.8895150738482495508826868606319e64*cos(r)*sin(e)*sin(r)/(9.8867616562362385368459 482397586e64*cos(r)2 + 2.5e65*cos(e)2*sin(r)2) x2 = 0.944733937415299 time=(x2-e-2.798053679+3.141592654)/0.261799388+14 time = 17.3572165434923 附录二：附录二： 附录三：附录三： JD=2452653.75+122+1+B-1524.52 JD = 2451734.23 T=(JD-2451545)/36525 T = 0.00518083504448957 h=4.881627973+628.3319509*T+0.000005235988*T2 h = 8.13691216393576 Ps=4.908229409+0.020453548*T+0.000078539816*T2+0.00000005235988*T3 Ps = 4.90833537756636 e=0.409319678+0.00022712*T-0.000000028604*T2+0.00000000887209*T3 e = 0.409320854670489 r=h+0.033501*sin(h-Ps)+0.000351*sin(2*(h-Ps)+0.000005*sin(3*(h-Ps) r = 8.13406124716575 x=solve(sin(39.95)*sin(e)*sin(r)+cos(39.95)*cos(r)*x+sqrt(1-x2)*cos(e)*sin(r)=0,x) x = (2.5*(100000000000000000000000000000000.0*cos(e)*sin(r)*(3.9547046624944954147383792 959034e65*cos(r)2 + 1.0e66*cos(e)2*sin(r)2 - 6.0452953375055045852616207040966e65*sin(e)2*sin(r)2)(1/2) + 4.8895150738482495508826868606319e64*cos(r)*sin(e)*sin(r)/(9.8867616562362385368459 482397586e64*cos(r)2 + 2.5e65*cos(e)2*sin(r)2) -(2.5*(100000000000000000000000000000000.0*cos(e)*sin(r)*(3.954704662494495414738379 2959034e65*cos(r)2 + 1.0e66*cos(e)2*sin(r)2 - 6.0452953375055045852616207040966e65*sin(e)2*sin(r)2)(1/2) - 4.8895150738482495508826868606319e64*cos(r)*sin(e)*sin(r)/(9.8867616562362385368459 482397586e64*cos(r)2 + 2.5e65*cos(e)2*sin(r)2) x1=(2.5*(100000000000000000000000000000000.0*cos(e)*sin(r)*(3.9547046624944954147383 792959034e65*cos(r)2 + 1.0e66*cos(e)2*sin(r)2 - 6.0452953375055045852616207040966e65*sin(e)2*sin(r)2)(1/2) + 4.8895150738482495508826868606319e64*cos(r)*sin(e)*sin(r)/(9.8867616562362385368459 482397586e64*cos(r)2 + 2.5e65*cos(e)2*sin(r)2) x1 = 0.861809443498602 xita=acos(x1) xita = 0.531970103985351 time=(x1-e-2.798053679+3.141592654)/0.261799388+14 time = 17.0406013165627 附录四：附录四： 1999/100 ans = 19.99 1999/4 ans = 499.75 A=19 A = 19 B=2-A+499 B = 482 365.25*(1999+4716) ans = 2452653.75 30.6001*(3+1) ans = 122.4004 JD=2452653.75+122+1+B-1524.52 JD = 2451734.23 T=(JD-2451545)/36525 T = 0.00518083504448957 s=4.72454834 +8399.70927457*T+0.00003456*T2+0.00000003*T3 s = 48.2420565141439 p=5.83515161 +71.02338662*T-0.00018012*T2+0.00000017*T3 p = 6.20311205554464 N=4.52360160 +33.