2023高教社杯全國大學生數學建模競賽承諾書我們仔細閱讀了?全國大學生數學建模競賽章程?和?全國大學生數學建模競賽參賽規那么?〔以下簡稱為“競賽章程和參賽規那么〞，可從全國大學生數學建模競賽網站下載〕。我們完全明白，在競賽開始后參賽隊員不能以任何方式〔包括、電子郵件、網上咨詢等〕與隊外的任何人〔包括指導教師〕研究、討論與賽題有關的問題。我們知道，抄襲別人的成果是違反競賽章程和參賽規那么的，如果引用別人的成果或其他公開的資料〔包括網上查到的資料〕，必須按照規定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽章程和參賽規那么，以保證競賽的公正、公平性。如有違反競賽章程和參賽規那么的行為，我們將受到嚴肅處理。我們授權全國大學生數學建模競賽組委會，可將我們的論文以任何形式進行公開展示〔包括進行網上公示，在書籍、期刊和其他媒體進行正式或非正式發表等〕。我們參賽選擇的題號是〔從A/B/C/D中選擇一項填寫〕：C我們的參賽報名號為〔如果賽區設置報名號的話〕：所屬學校〔請填寫完整的全名〕：長春工業大學參賽隊員(打印并簽名)：1.武太彬2.賈光芒3.牛文正指導教師或指導教師組負責人(打印并簽名)：李純潔 〔論文紙質版與電子版中的以上信息必須一致，只是電子版中無需簽名。以上內容請仔細核對，提交后將不再允許做任何修改。如填寫錯誤，論文可能被取消評獎資格。〕日期：2023年9月_16_日賽區評閱編號〔由賽區組委會評閱前進行編號〕：2023高教社杯全國大學生數學建模競賽編號專用頁賽區評閱編號〔由賽區組委會評閱前進行編號〕：賽區評閱記錄〔可供賽區評閱時使用〕：評閱人評分備注全國統一編號〔由賽區組委會送交全國前編號〕：全國評閱編號〔由全國組委會評閱前進行編號〕：古塔的變形摘要本文對古塔的變形問題建立數學模型，它實質上是一個空間解析幾何問題。首先建立空間解析幾何模型，并利用這個模型對問題1進行求解，然后對模型進行數據處理和圖形分析，精確地給出了中心位置坐標；之后對于問題2，在問題1的根底上我們計算出中心坐標的擬合直線，計算出曲率的值；最后對于問題3，運用AR自回歸模型給出古塔的變形趨勢。在問題1中，通過分析這四年古塔的每一層中的8個離散點，運用MATLAB數學軟件建立觀測數據模擬圖，依據此圖作出每一層相應散點的投影平面，經過計算根本近似正八邊形，得出正八邊形的中心，并運用空間直線擬合模型，求出古塔的一條軸心擬合直線；并且運用空間平面擬合模型，求出每一層的8個散點擬合的平面，那么直線與每一平面相交的點，即為每一層的中心坐標，以所求第一層中心坐標為例，1986年第一層中心坐標為，1996年第一層中心坐標為，2023年第一層中心坐標為，2023年第一層中心坐標，其它〔算上塔尖〕14層見正文表5。在問題2中，首先運用MATLAB軟件畫出每一年的俯視圖，即各年的平面圖，可以看出這四年的古塔是逐漸傾斜和彎曲的，且在第五層開始發生了扭曲。建立曲率數學模型，運用軟件求解出相應的曲率值。在問題3中，首先建立帶季節項的AR自回歸模型，運用問題2中所求的曲率作為自變量，代入到模型中，運用時間序列分析的統計軟件SPSS，得出古塔變形的趨勢，即隨年代的增長，古塔的傾斜和彎曲將更加嚴重，且在第五層由于空間中心的改變，將更加扭曲。在本文最后，對模型的優缺點及改良之處進行分析。關鍵詞：空間解析幾何最小二乘法擬合曲率AR預測模型MATLAB7.0軟件1問題的重述與分析由于許多外在原因，古塔會產生各種變形，傾斜，彎曲，扭曲等。為了了解各種變形量，測繪公司先后于1986年7月，1996年8月，2023年3月和2023年3月對該塔進行了四次觀測。討論3個問題，問題1，給出確定古塔各層中心位置的通用方法，列表給出各次測量的古塔各層中心坐標。本文首先對古塔的變形情況進行分析，可以獲取位置信息，且只有四次的觀測數據信息。建立適當的坐標系，先研究一層塔的八個點的投影所得的八個點的坐標點，然后再確定各中心的坐標，從而找到各層的中心點，可以擬合成一條直線。并且運用空間擬合平面模型，求出每一層的8個散點擬合的平面，那么直線與每一平面相交的點，即為每一層的中心坐標。問題2，分析該塔傾斜，彎曲，扭曲等變形情況。首先運用MATLAB軟件畫出每一年的俯視圖，即各年的平面圖，可以看出四年的古塔是傾斜還是彎曲或是扭曲，建立相應的數學模型。問題3，分析該塔的變形趨勢。首先建立帶季節項的AR自回歸模型，運用問題2中所求的結果作為自變量，代入到模型中，得出古塔變形的趨勢。2根本假設1.觀測的數據準確；

2.以古塔的地面為水平面；

3.假設人在觀測時是圍繞塔的中心進行八個不同的方面進行大量的測量；4.兩點的距離作到小數點后百分位即為等長，后面忽略；5.不考慮異常點；3符號說明每個點的橫坐標；每個點的縱坐標；每個點的豎坐標；直線方程的待估系數；空間直線的對應的方向向量；：分別是邊和邊由其原始位置旋轉的角度；:變化率；4模型的建立與求解4.1問題1的求解首先根據每一年數據中的每一層的8個散點坐標運用MATLAB軟件，得到如下每年觀測數據的模擬圖。圖11986年觀測數據模擬圖圖21996年觀測數據模擬圖圖32023年觀測數據模擬圖圖42023年觀測數據模擬圖由圖1到圖4的模擬圖可以直觀地看到，四次測量古塔的顯著變化，相同坐標上塔的上部越來越小，把圖2和圖3疊加到一個圖。圖5古塔扭曲變化比照圖圖5中，紅線代表1996年的塔圖，黃線代表2023年的塔圖。兩個塔圖疊加在一起比照更加顯明，而在第五層出現扭曲，從第五層往上越來越小，到2023年成了一個點的塔尖。經過以上的圖形分析，本文以空間解析幾何為背景，運用最小二乘擬合理論給出空間直線與平面的擬合。下面建立空間內曲線的最小二乘問題，以空間直線為例研究曲線擬合的方法。空間直線的標準方程整理得直線射影式方程，其中;;這樣直線可以看作是用這2個方程表示的平面相交的直線，所以可以分別對2個方程進行數據擬合。設表示按擬合方程求得的近似值。一般地，它不同于實測值兩者之差，同理可得當Q取最小值時的值即為方程的系數，即滿足以下方程時Q值最小有令(1)方程組(1)可寫成其中根據m組數據點解方程組就可以求得的值。圖61996年古塔各層中心點擬合直線分析圖圖72023年古塔各層中心點擬合函數分析圖從圖6，圖7中可以得出：由1996年古塔各層中心點擬合直線與2023年古塔各層中心點擬合直線，方程為，由一次方程的系數k看出1996年為0.0105,2023年為0.0112，可以得出古塔的中心發生了一次傾斜，且2023年比1986年傾斜角度變大。運用附件中MATLAB的程序經過調試，得到如下結果：1986年1996年2023年2023年接下來運用最小二乘法擬合求空間平面方程：設空間平面方程為，其中為待估參數。設古塔的每層8個點，設點坐標，考慮到數據在三個方向均存在誤差，得矩陣形式：即。通常采用矩陣奇異值分解解算待定參數的整體最小二乘解。其中那么參數的整體最小二乘估計為：殘差矩陣為：。運用MATLAB軟件得到古塔各層平面方程，列表為：表11986年的擬合平面方程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.620232)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)表21996年的擬合平面方程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)表32023年的擬合平面方程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)表42023年的擬合平面方程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年2023年2023年1層2層3層4層5層6層7層8層9層10層11層12層13層4.2問題2的求解問題2是為了分析古塔的變形情況，根據問題1的圖1-圖5，運用MATLAB軟件，得到如下關于軸的俯視圖：圖8古塔模擬俯視圖圖91986年觀測數據俯視模擬圖圖101996年觀測數據模擬俯視圖圖112023年觀測數據模擬俯視圖圖8至圖11分別給出了整體和各個年份的俯視圖，從四個圖中可以看到并不是均勻的八邊形環，有些邊環密，有些邊環稀，這就直觀地說明了古塔發生了傾斜和彎曲，加上問題1的第五層的中心坐標的改變，說明古塔在第五層上也發生了扭曲現象。本文在第一問最后得出每層各中心坐標的擬合直線，運用到第二問中古塔的傾斜問題上，傾斜是指根底兩端點傾斜方向的沉降差與其距離的比值建筑中心線或其墻、柱等，在不同高度的點對其相應底部點的偏移現象。運用公式兩條直線夾角公式：彎曲，即不直。可以分為形變彎曲及空間彎曲。當桿件受到與桿軸線垂直的外力或在軸線平面內的力偶作用時，桿的軸線由原來的直線變成曲線，這種變形叫彎曲變形。曲率處處不為零的空間稱為彎曲空間。物體發生彎曲時產生的形變叫做“彎曲形變〞。物體彎曲得越厲害，產生的彈力就越大。例如，將弓拉得越滿，箭就射得越遠。把一個物體放在支持物上，物體越重，支持物被壓彎曲得越厲害，支持力就越大。運用曲率公式：扭曲物體因外力作用而扭轉變形,也用于比喻數學的一種改變物體形狀的方法，扭曲變形是在彎曲變形的根底上，旋轉，扭成螺旋狀態，叫扭曲變形。問題1中得出扭曲的發生主要在第五層。首先,計算方向的扭曲變形。由于邊和的旋轉角度一般都很小(對于一個實際問題,一般因此可以認為約等于,因此:式中:分別是邊和邊由其原始位置旋轉的角度。如果，那么說明微元的一邊產生均勻豎向位移,那么變形后各點仍保持在一個平面內,意味著這個微元發生的是剛性的旋轉。否那么,當,即在兩條邊之間存在一個旋轉角的差，就說明產生了扭曲變形。那么在變形后的面上,扭曲變形的大小可以采用這兩個角度之差隨著兩個對邊之間距離的變化率來表示:或者:如果應用以上分析過程考慮y方向的扭曲變形,即首先計算微元另外兩條對邊和的旋轉角度的差,將會得到相同的結果。最終運用MATLAB軟件得出4.3問題3的求解在問題1的根底上，運用MATLAB軟件，給出1996年與2023年的中心擬合直線圖來。時間序列分析的目標就是通過分析要素〔變量〕隨時間變化的歷史過程，揭示其變化開展規律，并對未來狀態進行分析預測。如在變形測量中，可以采用時間序列分析方法對觀測數據進行分析，以便建立變形體的動態變形預測模型，并對其變形趨勢進行預測。而自回歸(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]姚澤清,鄭旭東,全國大學生數學建模賽題優秀論文評析<2005年—2023年A題>,北京:國防工業出版社，2023:198-218.[2]周培德，計算幾何[M].北京：清華大學出版社，2023：100-105.[3]襲楊，空間直線擬合的一種方法[J].齊齊哈爾大學學報，2023,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]．武漢:武漢大學出版社，2004．[9]俞錦成．關于整體最小二乘的可解性[J]．南京師范大學學報(自然科學版)，1996，19(1):13-16．7附錄直線擬合：>>x=[566.6648566.7196566.7735566.8161566.8621566.9084566.9467566.9843567.0218567.0569567.1045567.1518];y=[522.7105522.6684522.6273522.5944522.5591522.5244522.5081522.4924522.4764522.4624522.423522.3836];z=[1.78747.320212.755217.078321.720526.235129.836933.350936.854940.172144.440948.7119];F=[z;111111111111];M=F*F';N=F*x';O=F*y';A=(M\N)'B=(M\O)'A=0.0104566.6411B=-0.0066522.7125//1986年>>x=[566.665566.7205566.7751566.8183566.8649566.9118566.9506566.9884567.0265567.062567.1102567.1578];y=[522.7102522.6674522.6256522.5922522.5563522.521522.5042522.4881522.4714522.4572522.4173522.3775];z=[1.7837.314612.750717.075121.71626.229529.832333.345436.848340.167644.435448.7074];F=[z;111111111111];M=F*F';N=F*x';O=F*y';A=(M\N)'B=(M\O)'A=0.0105566.6411B=-0.0067522.7124//1996年>>x=[566.727566.7642566.8004566.8297566.861566.9478566.98567.0313567.0825567.1381567.181567.2238];y=[522.7014522.669522.6387522.6127522.586522.5335522.5115522.4788522.4457522.3926522.3535522.3147];z=[1.76327.290512.726917.05221.703926.204529.81733.336636.822240.144144.424948.6839];F=[z;111111111111];M=F*F';N=F*x';O=F*y';A=(M\N)'B=(M\O)'A=0.0112566.6643B=-0.0084522.7436//2023年>>x=[566.7268566.764566.8001566.8293566.8603566.9471566.9792567.0305567.0816567.137567.1799567.2225];y=[522.7015522.6693522.6384522.6132522.5866522.5342522.5123522.4797522.4466522.3937522.3547522.316];z=[1.76457.30912.732317.069721.709426.21129.824633.339936.843840.161144.432648.6998];F=[z;111111111111];M=F*F';N=F*x';O=F*y';A=(M\N)'B=(M\O)'A=0.0112566.6642B=-0.0084522.7436//2023年平面擬合公式：86.1z=(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.620232)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.1z=(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.1z=(-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.1z=(-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');holdon;X=[ones(8,1)xy];b=regress(z,X);holdon;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)));兩直線夾角：functiona=JiaJiao(x,y)

%求兩條直線夾角

%x,y是三點的橫坐標和縱坐標

%eg:x=[123];y=[415];

ifx(2)~=x(1)

k1=(y(2)-y(1))/(x(2)-x(1));

end

ifx(3)~=x(2)

k2=(y(3)-y(2))/(x(3)-x(2));

end

ifx(2)==x(1)&x(3)==x(2)

a=0;

elseifx(3)==x(2)

a=pi/2-atan(abs(k1));

elseifx(1)==x(2)

a=pi/2-atan(abs(k2));

elseif1+k1*k2==0

a=pi/2;

else

a=atan(abs((k2-k1)/(1+k2*k1)));%夾角

end

a=a*360/(2*pi);%轉化為角度制向量的夾角：A=[129];

B=[010];

acos(dot(A,B)/(norm(A)*norm(B)))*180/pi曲率:

clc;clearall;closeall;

x0=linspace(0,1);

y0=sin(x0).*cos(x0);

h=abs(diff([x0(2),x0(1)]));

%模擬一階導

figure;boxon;holdon;

ythe1=cos(x0).^2-sin(x0).^2;%理論一階導

yapp1=gradient(y0,h);%matlab數值近似

plot(x0,ythe1,'.');

plot(x0,yapp1,'r');

legend('理論值','模擬值');

title('模擬一階導');

%模擬二階導

figure;boxon;holdon;

ythe2=(-4)*cos(x0).*sin(x0);%理論二階導

yapp2=2*2*del2(y0,h);%matlab數值近似

plot(x0,ythe2,'.');

plot(x0,yapp2,'r');

legend('理論值','模擬值');

title('模擬二階導');

%模擬曲率

symsxy

y=sin(x)*cos(x);

yd2=diff(y,2);

yd1=diff(y,1);

k=abs(yd2)/(1+yd1^2)^(3/2);

k1=subs(k,x,