2013年全国大学生高教杯数学建模C题古塔变形论文.doc

2013高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了全国大学生数学建模竞赛章程和全国大学生数学建模竞赛参赛规则（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 长春工业大学 参赛队员 (打印并签名) ：1. 武太彬 2. 贾光辉 3. 牛文正 指导教师或指导教师组负责人 (打印并签名)： 李纯净 （论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。） 日期： 2013 年 9 月_16_日赛区评阅编号（由赛区组委会评阅前进行编号）：2013高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：古塔的变形摘 要本文对古塔的变形问题建立数学模型，它实质上是一个空间解析几何问题。首先建立空间解析几何模型，并利用这个模型对问题1进行求解，然后对模型进行数据处理和图形分析，精确地给出了中心位置坐标；之后对于问题2，在问题1的基础上我们计算出中心坐标的拟合直线，计算出曲率的值；最后对于问题3，运用AR自回归模型给出古塔的变形趋势。在问题1中，通过分析这四年古塔的每一层中的8 个离散点，运用MATLAB数学软件建立观测数据模拟图，依据此图作出每一层相应散点的投影平面，经过计算基本近似正八边形，得出正八边形的中心，并运用空间直线拟合模型，求出古塔的一条轴心拟合直线；并且运用空间平面拟合模型，求出每一层的8个散点拟合的平面，那么直线与每一平面相交的点，即为每一层的中心坐标，以所求第一层中心坐标为例，1986年第一层中心坐标为，1996年第一层中心坐标为，2009年第一层中心坐标为，2011年第一层中心坐标，其它（算上塔尖）14层见正文表5。在问题2中，首先运用MATLAB软件画出每一年的俯视图，即各年的平面图，可以看出这四年的古塔是逐渐倾斜和弯曲的，且在第五层开始发生了扭曲。建立曲率数学模型，运用软件求解出相应的曲率值。在问题3中，首先建立带季节项的AR自回归模型，运用问题2中所求的曲率作为自变量，代入到模型中，运用时间序列分析的统计软件SPSS，得出古塔变形的趋势，即随年代的增长，古塔的倾斜和弯曲将更加严重，且在第五层由于空间中心的改变，将更加扭曲。在本文最后，对模型的优缺点及改进之处进行分析。关键词：空间解析几何 最小二乘法拟合 曲率 AR预测模型 MATLAB7.0软件1问题的重述与分析由于许多外在原因，古塔会产生各种变形，倾斜，弯曲，扭曲等。为了了解各种变形量，测绘公司先后于1986年7月，1996年8月，2009年3月和2011年3月对该塔进行了四次观测。讨论3个问题，问题1，给出确定古塔各层中心位置的通用方法，列表给出各次测量的古塔各层中心坐标。本文首先对古塔的变形情况进行分析，可以获取位置信息，且只有四次的观测数据信息。建立适当的坐标系，先研究一层塔的八个点的投影所得的八个点的坐标点，然后再确定各中心的坐标，从而找到各层的中心点，可以拟合成一条直线。并且运用空间拟合平面模型，求出每一层的8个散点拟合的平面，那么直线与每一平面相交的点，即为每一层的中心坐标。 问题2，分析该塔倾斜，弯曲，扭曲等变形情况。首先运用MATLAB软件画出每一年的俯视图，即各年的平面图，可以看出四年的古塔是倾斜还是弯曲或是扭曲，建立相应的数学模型。问题3，分析该塔的变形趋势。首先建立带季节项的AR自回归模型，运用问题2中所求的结果作为自变量，代入到模型中，得出古塔变形的趋势。2 基本假设1. 观测的数据准确；2. 以古塔的地面为水平面；3. 假设人在观测时是围绕塔的中心进行八个不同的方面进行大量的测量；4. 两点的距离作到小数点后百分位即为等长，后面忽略；5. 不考虑异常点；3 符号说明已知每个点的横坐标；已知每个点的纵坐标；已知每个点的竖坐标；直线方程的待估系数；空间直线的对应的方向向量；：分别是边和边由其原始位置旋转的角度；: 变化率；4 模型的建立与求解4.1 问题1的求解 首先根据已知每一年数据中的每一层的8个散点坐标运用MATLAB软件，得到如下每年观测数据的模拟图。图1 1986年观测数据模拟图图2 1996年观测数据模拟图图3 2009年观测数据模拟图图4 2011年观测数据模拟图由图1到图4的模拟图可以直观地看到，四次测量古塔的显著变化，相同坐标上塔的上部越来越小，把图2和图3叠加到一个图。图5 古塔扭曲变化对比图图5中，红线代表1996年的塔图，黄线代表2009年的塔图。两个塔图叠加在一起对比更加显明，而在第五层出现扭曲，从第五层往上越来越小，到2011年成了一个点的塔尖。 经过以上的图形分析，本文以空间解析几何为背景，运用最小二乘拟合理论给出空间直线与平面的拟合。下面建立空间内曲线的最小二乘问题，以空间直线为例研究曲线拟合的方法。 已知空间直线的标准方程整理得直线射影式方程 ，其中;这样直线可以看作是用这 2 个方程表示的平面相交的直线，所以可以分别对 2 个方程进行数据拟合。 设表示按拟合方程求得的近似值。一般地，它不同于实测值两者之差，同理可得 当Q取最小值时的值即为方程的系数，即满足下列方程时Q值最小有令 (1)方程组(1)可写成 其中根据m组数据点解方程组就可以求得的值。图6 1996年古塔各层中心点拟合直线分析图图7 2009年古塔各层中心点拟合函数分析图从图6，图7中可以得出：由1996年古塔各层中心点拟合直线与2009年古塔各层中心点拟合直线，方程为，由一次方程的系数k看出1996年为0.0105,2009年为0.0112，可以得出古塔的中心发生了一次倾斜，且2009年比1986年倾斜角度变大。运用附件中MATLAB的程序经过调试，得到如下结果： 1986年 1996年 2009年 2011年接下来运用最小二乘法拟合求空间平面方程：设空间平面方程为，其中为待估参数。设古塔的每层8个点，设点坐标，考虑到数据在三个方向均存在误差，得矩阵形式：即。通常采用矩阵奇异值分解解算待定参数的整体最小二乘解。其中则参数的整体最小二乘估计为：残差矩阵为：。运用MATLAB软件得到古塔各层平面方程，列表为：表1 1986年的拟合平面方程z=(0.003417)*x+(-0.000831)*y+(0.471956) z=(0.003629)*x+(-0.000818)*y+(5.887189) z=(0.003726)*x+(-0.000883)*y+(11.308475) z=(0.003838)*x+(-0.000935)*y+(15.602646) z=(0.003941)*x+(-0.000874)*y+(20.156847) z=(0.005443)*x+(-0.017498)*y+(33.310426) z=(0.005668)*x+(-0.018744)*y+(37.502066) z=(0.005887)*x+(-0.020009)*y+(41.620122) z=(0.006115)*x+(-0.021455)*y+(45.825845) z=(0.000887)*x+(-0.021682)*y+(52.003938) z=(0.001868)*x+(-0.022189)*y+(56.048295) z=(0.002175)*x+(-0.023884)*y+(61.121714) z=(0.004664)*x+(-0.027610)*y+(66.050359)表2 1996年的拟合平面方程z=(0.003715)*x+(-0.001488)*y+(0.684400)z=(0.003782)*x+(-0.000493)*y+(5.617203)z=(0.003735)*x+(-0.001266)*y+(11.515906)z=(0.004045)*x+(0.000715)*y+(14.555943)z=(0.003941)*x+(-0.001286)*y+(20.385763)z=(0.005639)*x+(-0.017135)*y+(32.997144)z=(0.005880)*x+(-0.018358)*y+(37.168214)z=(0.005877)*x+(-0.020489)*y+(41.891314)z=(0.006366)*x+(-0.021014)*y+(45.437916)z=(0.000873)*x+(-0.022225)*y+(52.314084)z=(0.001853)*x+(-0.022770)*y+(56.380505)z=(0.002159)*x+(-0.024511)*y+(61.481252)z=(0.005234)*x+(-0.027989)*y+(65.963401)表3 2009年的拟合平面方程z=(-0.003770)*x+(-0.002373)*y+(5.079985)z=(-0.000194)*x+(-0.003972)*y+(9.661222)z=(-0.003697)*x+(-0.002317)*y+(15.977842)z=(-0.000065)*x+(-0.004433)*y+(19.616420)z=(-0.000946)*x+(-0.004248)*y+(24.611974)z=(-0.020936)*x+(-0.004708)*y+(39.819633)z=(-0.018595)*x+(-0.006811)*y+(43.402884)z=(-0.021604)*x+(-0.006530)*y+(48.330004)z=(-0.022140)*x+(-0.006879)*y+(52.311658)z=(-0.022761)*x+(-0.001523)*y+(52.915196)z=(-0.022023)*x+(-0.002515)*y+(57.363039)z=(-0.023948)*x+(-0.002804)*y+(62.798782)z=(-0.027663)*x+(-0.005448)*y+(70.356472)表4 2011年的拟合平面方程z=(-0.003740)*x+(-0.002333)*y+(5.040576)z=(-0.002530)*x+(-0.002555)*y+(10.061014)z=(-0.004408)*x+(-0.002633)*y+(16.522721)z=(-0.002736)*x+(-0.003593)*y+(20.518390)z=(-0.001351)*x+(-0.004259)*y+(24.823896)z=(-0.021377)*x+(-0.004712)*y+(40.045928)z=(-0.018205)*x+(-0.007014)*y+(43.305695)z=(-0.021560)*x+(-0.006476)*y+(48.273096)z=(-0.025300)*x+(-0.004743)*y+(52.729697)z=(-0.026175)*x+(0.000449)*y+(53.562983)z=(-0.024393)*x+(-0.002167)*y+(58.395736)z=(-0.018681)*x+(-0.004368)*y+(60.918705)z=(-0.028368)*x+(-0.005161)*y+(70.556460)联立直线与平面方程，即得到表5的各年各层的中心坐标。表5 各年各层的中心坐标1986年1996年2009年2011年1层2层3层4层5层6层7层8层9层10层11层12层13层4.2 问题2的求解 问题2是为了分析古塔的变形情况，根据问题1的图1-图5，运用MATLAB软件，得到如下关于轴的俯视图：图8 古塔模拟俯视图图9 1986年观测数据俯视模拟图图10 1996年观测数据模拟俯视图图11 2009年观测数据模拟俯视图 图8至图11分别给出了整体和各个年份的俯视图，从四个图中可以看到并不是均匀的八边形环，有些边环密，有些边环稀，这就直观地说明了古塔发生了倾斜和弯曲，加上问题1的第五层的中心坐标的改变，说明古塔在第五层上也发生了扭曲现象。本文在第一问最后得出每层各中心坐标的拟合直线，运用到第二问中古塔的倾斜问题上，倾斜是指基础两端点倾斜方向的沉降差与其距离的比值建筑中心线或其墙、柱等，在不同高度的点对其相应底部点的偏移现象。运用公式两条直线夹角公式：弯曲，即不直。可以分为形变弯曲及空间弯曲。当杆件受到与杆轴线垂直的外力或在轴线平面内的力偶作用时，杆的轴线由原来的直线变成曲线，这种变形叫弯曲变形。曲率处处不为零的空间称为弯曲空间。物体发生弯曲时产生的形变叫做“弯曲形变”。物体弯曲得越厉害，产生的弹力就越大。例如，将弓拉得越满，箭就射得越远。把一个物体放在支持物上，物体越重，支持物被压弯曲得越厉害，支持力就越大。运用曲率公式：扭曲物体因外力作用而扭转变形,也用于比喻数学的一种改变物体形状的方法，扭曲变形是在弯曲变形的基础上，旋转，扭成螺旋状态，叫扭曲变形。问题1中得出扭曲的发生主要在第五层。首先,计算方向的扭曲变形。由于边和的旋转角度一般都很小(对于一个实际问题,一般因此可以认为约等于,因此:式中:分别是边和边由其原始位置旋转的角度。如果，则说明微元的一边产生均匀竖向位移,则变形后各点仍保持在一个平面内,意味着这个微元发生的是刚性的旋转。否则,当,即在两条边之间存在一个旋转角的差，就说明产生了扭曲变形。则在变形后的面上,扭曲变形的大小可以采用这两个角度之差随着两个对边之间距离的变化率来表示:或者:如果应用以上分析过程考虑y方向的扭曲变形,即首先计算微元另外两条对边和的旋转角度的差,将会得到相同的结果。最终运用MATLAB软件得出4.3 问题3的求解在问题1的基础上，运用MATLAB软件，给出1996年与2009年的中心拟合直线图来。时间序列分析的目标就是通过分析要素（变量）随时间变化的历史过程，揭示其变化发展规律，并对未来状态进行分析预测。如在变形测量中，可以采用时间序列分析方法对观测数据进行分析，以便建立变形体的动态变形预测模型，并对其变形趋势进行预测。而自回归( AR) 模型的参数估计是时间序列分析的基本问题，是在模型结构及阶次已确定的条件下，对模型参数进行估计，使所建立的模型是实际时间序列的“最佳”拟合模型。但在实际的观测中，观测值是由一定观测手段得到的，不可避免地含有随机误差。理论上讲AR( p) 参数的最小二乘估计也就是线性最小二乘估计，所以可以用现在广泛应用的整体最小二乘估计方法进行估计，即不仅考虑自身观测值的误差，同时考虑与其有关的自身前一个或前几个时刻的观测值的误差，从而进行参数估计，得到的预测观测数据更准确。但AR 模型是典型的预测模型，而最小二乘估计是常用的数据分析工具，在预测数据上还是常应用最小二乘估计比较准确，反而最小二乘估计的数据不够准确和稳定。自回归模型，子样观察值，白噪声序列表示为，回归系数用表示，则可以得到模型： ，模型参数的最小二乘估计。设样本观测值，记,则模型可表示为有最小二乘原理可得到模型参数的估计为。那么根据最小二乘估计值可以得到噪声的估计值为噪声方差为，由此模型可以求出未来的变形趋势。5模型的评价及推广5.1模型的优点(1). 在问题的求解中，充分运用了数据和图形，即表格，使结果明了清晰。(2). 本文采用了多种专业软件对模型进行求解，如Matlab，SPSS等提高了模型 的精确度。(3). 本文所有模型建立均完全基于实际的统计数据，拟合成近似的线进行处理， 科学合理。(4). 本文建立多种模型，多种模型的对比更能体现不同模型的特点及优势。(5). 应用相关几何知识，将空间复杂的点分布简化成直线和面，应用几何问题求 空间点。5.2模型的缺点(1). 背景资料的筛选方法有待进一步优化和改进；(2). 模型建立中采用的数据不太准确，利用拟合估计值，对模型产生了相对较大 的误差。(3). 建立的模型不够完善，忽略了许多的因素，过于简单。(4). 由于古塔在现实生活中受到各种影响，理论值难免与实际情况有所偏差。5.3模型的改进 由于都是求出每个点都是一个估计值，从而有较大的误差，需要在模型的建立中进行改进，使估计的所有的数据更加准确。 5.4模型的推广用中心位置的模型方法，不难预测全国所有古塔或高楼由于外力和内力的作用在以后的变化趋势。6.参考文献1 姚泽清,郑旭东,全国大学生数学建模赛题优秀论文评析,北京:国防工业出版社，2012: 198-218.2 周培德，计算几何M. 北京：清华大学出版社，2011：100-105.3 袭杨，空间直线拟合的一种方法J.齐齐哈尔大学学报，2009,25(2):64-68.4 王能超. 数值分析简明教程M. 北京：高等教育出版社，2003：36-37. 5 吕林根，许子道. 解析几何M. 北京：高等教育出版社，2004：156-157.6 于晓秋，李欣，张宏礼. MATLAB 数值实验M. 沈阳：辽宁科学技术出版社，2004：34-40. 7 王正东. 数学软件与数学实验M. 北京：科学出版社，2004：88-89. 8 吴怀宇 时间序列分析与综合M 武汉: 武汉大学出版社， 20049 俞锦成 关于整体最小二乘的可解性J 南京师范大学学报( 自然科学版) ， 1996， 19( 1) : 13-167 附录直线拟合： x=566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9843 567.0218 567.0569 567.1045 567.1518; y=522.7105 522.6684 522.6273 522.5944 522.5591 522.5244 522.5081 522.4924 522.4764 522.4624 522.423 522.3836; z=1.7874 7.3202 12.7552 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A = 0.0104 566.6411B = -0.0066 522.7125 /1986年 x=566.665 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.062 567.1102 567.1578; y=522.7102 522.6674 522.6256 522.5922 522.5563 522.521 522.5042 522.4881 522.4714 522.4572 522.4173 522.3775; z=1.783 7.3146 12.7507 17.0751 21.716 26.2295 29.8323 33.3454 36.8483 40.1676 44.4354 48.7074; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A = 0.0105 566.6411B = -0.0067 522.7124 /1996年 x=566.727 566.7642 566.8004 566.8297 566.861 566.9478 566.98 567.0313 567.0825 567.1381 567.181 567.2238; y=522.7014 522.669 522.6387 522.6127 522.586 522.5335 522.5115 522.4788 522.4457 522.3926 522.3535 522.3147; z=1.7632 7.2905 12.7269 17.052 21.7039 26.2045 29.817 33.3366 36.8222 40.1441 44.4249 48.6839; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A =0.0112 566.6643B =-0.0084 522.7436 /2011年 x=566.7268 566.764 566.8001 566.8293 566.8603 566.9471 566.9792 567.0305 567.0816 567.137 567.1799 567.2225; y=522.7015 522.6693 522.6384 522.6132 522.5866 522.5342 522.5123 522.4797 522.4466 522.3937 522.3547 522.316; z=1.7645 7.309 12.7323 17.0697 21.7094 26.211 29.8246 33.3399 36.8438 40.1611 44.4326 48.6998; F=z;1 1 1 1 1 1 1 1 1 1 1 1; M=F*F; N=F*x; O=F*y; A=(MN) B=(MO) A =0.0112 566.6642B = -0.0084 522.7436 /2009年平面拟合公式：86.1 z=(0.003417)*x+(-0.000831)*y+(0.471956) z=(0.003629)*x+(-0.000818)*y+(5.887189) z=(0.003726)*x+(-0.000883)*y+(11.308475) z=(0.003838)*x+(-0.000935)*y+(15.602646) z=(0.003941)*x+(-0.000874)*y+(20.156847) z=(0.005443)*x+(-0.017498)*y+(33.310426) z=(0.005668)*x+(-0.018744)*y+(37.502066) z=(0.005887)*x+(-0.020009)*y+(41.620122) z=(0.006115)*x+(-0.021455)*y+(45.825845) z=(0.000887)*x+(-0.021682)*y+(52.003938) z=(0.001868)*x+(-0.022189)*y+(56.048295) z=(0.002175)*x+(-0.023884)*y+(61.121714) z=(0.004664)*x+(-0.027610)*y+(66.050359)96.1 z=(0.003715)*x+(-0.001488)*y+(0.684400) z=(0.003782)*x+(-0.000493)*y+(5.617203)z=(0.003735)*x+(-0.001266)*y+(11.515906) z=(0.004045)*x+(0.000715)*y+(14.555943) z=(0.003941)*x+(-0.001286)*y+(20.385763) z=(0.005639)*x+(-0.017135)*y+(32.997144) z=(0.005880)*x+(-0.018358)*y+(37.168214) z=(0.005877)*x+(-0.020489)*y+(41.891314) z=(0.006366)*x+(-0.021014)*y+(45.437916)z=(0.000873)*x+(-0.022225)*y+(52.314084) z=(0.001853)*x+(-0.022770)*y+(56.380505) z=(0.002159)*x+(-0.024511)*y+(61.481252) z=(0.005234)*x+(-0.027989)*y+(65.963401)99.1 z=(-0.003770)*x+(-0.002373)*y+(5.079985) z=(-0.000194)*x+(-0.003972)*y+(9.661222) z=(-0.003697)*x+(-0.002317)*y+(15.977842) z=(-0.000065)*x+(-0.004433)*y+(19.616420) z=(-0.000946)*x+(-0.004248)*y+(24.611974) z=(-0.020936)*x+(-0.004708)*y+(39.819633) z=(-0.018595)*x+(-0.006811)*y+(43.402884) z=(-0.021604)*x+(-0.006530)*y+(48.330004) z=(-0.022140)*x+(-0.006879)*y+(52.311658) z=(-0.022761)*x+(-0.001523)*y+(52.915196) z=(-0.022023)*x+(-0.002515)*y+(57.363039) z=(-0.023948)*x+(-0.002804)*y+(62.798782) z=(-0.027663)*x+(-0.005448)*y+(70.356472)11.1 z=(-0.003740)*x+(-0.002333)*y+(5.040576) z=(-0.002530)*x+(-0.002555)*y+(10.061014) z=(-0.004408)*x+(-0.002633)*y+(16.522721)z=(-0.002736)*x+(-0.003593)*y+(20.518390) z=(-0.001351)*x+(-0.004259)*y+(24.823896) z=(-0.021377)*x+(-0.004712)*y+(40.045928) z=(-0.018205)*x+(-0.007014)*y+(43.305695) z=(-0.021560)*x+(-0.006476)*y+(48.273096)z=(-0.025300)*x+(-0.004743)*y+(52.729697) z=(-0.026175)*x+(0.000449)*y+(53.562983) z=(-0.024393)*x+(-0.002167)*y+(58.395736) z=(-0.018681)*x+(-0.004368)*y+(60.918705) z=(-0.028368)*x+(-0.005161)*y+(70.556460)平面拟合： x=xlsread(D:data.xls,C4:C11);y=xlsread(D:data.xls,D4:D11);z=xlsread(D:data.xls,E4:E11);plot3(x,y,z);scatter3(x,y,z,filled);hold on;X = ones(8,1) x y; b = regress(z,X);hold on;xfit = min(x):0.1:max(x); yfit = min(y):0.1:max(y);XFIT,YFIT= meshgrid (xfit,yfit); ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT;mesh (XFIT,YFIT,ZFIT);title(sprintf( z=(%f)*x+(%f)*y+(%f),b(3), b(2),b(1);两直线夹角：function a=JiaJiao(x,y)% 求两条直线夹角% x,y 是已知三点的横坐标和纵坐标% eg: x=1 2 3;y=4 1 5;if x(2)=x(1) k1=(y(2)-y(1)/(x(2)-x(1);endif x(3)=x(2)k2=(y(3)-y(2)/(x(3)-x(2);endif x(2)=x(1) & x(3)=x(2)a=0;elseif x(3)=x(2) a=pi/2-atan(abs(k1);elseif x(1)=x(2)a=pi/2-atan(abs(k2);elseif 1+k1*k2=0a=pi/2;elsea=atan(abs(k2-k1)/(1+k2*k1); % 夹角enda=a*360/(2*pi); % 转化为角度制向量的夹角：A=1 2 9;B=0 1 0;acos(dot(A,B)/(norm(A)*norm(B)*180/pi曲率:clc; clear all; close all;x0 = linspace(0, 1);y0 = sin(x0).*cos(x0);h = abs(diff(x0(2), x0(1);% 模拟一阶导figure; box on; hold on;ythe1 = cos(x0).2 - sin(x0).2; %理论一阶导yapp1 = gradient(y0, h); %matlab数值近似plot(x0, ythe1, .);plot(x0, yapp1, r);legend(理论值, 模拟值);title(模拟一阶导);% 模拟二阶导figure; box on; hold on;ythe2 = (-4)*cos(x0).*sin(x0); %理论二阶导yapp2 = 2*2*del2(y0, h); %matlab数值近似plot(x0, ythe2,.);plot(x0, yapp2,r);legend(理论值, 模拟值);title(模拟二阶导);% 模拟曲率syms x yy = sin(x)*cos(x);yd2 = diff(y, 2);yd1 = diff(y, 1);k = abs(yd2)/(1+yd12)(3/2);k1 = subs(k, x, x0);k2 = abs(yapp2)./(1+