范玲玲，樂騰勝，童智能，項炳泉，張金鐘

（1.南昌市京東學(xué)校，江西 南昌 330029；2.安徽省建筑科學(xué)研究設(shè)計院，安徽 合肥 230001；3.江西科技師范大學(xué)建筑工程學(xué)院，江西 南昌 330013）

1 引言

數(shù)學(xué)建模就是構(gòu)造數(shù)學(xué)模型的過程，即用數(shù)學(xué)的語言-公式、符號、圖表等刻畫和描述一個實際問題，然后經(jīng)過數(shù)學(xué)的處理-計算、迭代等得到定量的結(jié)果，以供人們做分析、預(yù)報、決策和控制[1]。隨著計算機廣泛使用與科學(xué)技術(shù)迅速發(fā)展，科學(xué)計算已是科學(xué)研究、工程設(shè)計中的一個重要的手段。合理的利用MATLAB數(shù)學(xué)軟件來輔助數(shù)學(xué)演算和繪圖，已成為與理論分析、科學(xué)實驗并駕齊驅(qū)的科學(xué)研究方法[2～5]。懸臂梁在各種荷載及不同支撐下的擾度分析是工民建中常見的問題，主要涉及到了常微分方程初值問題的數(shù)值求解。對于一般的初值問題，通常采用改進(jìn)的Euler公式，就能保證二階精度[6-7]。如果方程右端項f（x，y）足夠光滑，計算精度較高，經(jīng)典四階RK方法是不錯的方法，這一點在分析懸臂梁的擾度問題中得到了很好的應(yīng)用[8-9]。

本文針對懸臂梁的擾度隨著梁的長度的變化而變化，在MATLAB中運用經(jīng)典四階RK公式計算，作出懸臂梁擾度圖（z-y圖），確定出懸臂梁中擾度最大的位置以及求出相對應(yīng)的最大擾度近似解，并與撓度微分方程計算得出的精確解相對比，驗證了運用MATLAB軟件近似計算懸臂梁擾度的可行性，促進(jìn)了懸臂梁擾度計算的研究與分析。

2 經(jīng)典四階龍格-庫塔算法的介紹

經(jīng)典龍格-庫塔（RK）方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進(jìn)行抑制，所以其實現(xiàn)原理也較復(fù)雜。同Euler等算法一樣，該算法也是構(gòu)建在數(shù)學(xué)支持的基礎(chǔ)之上的。對于一階精度的歐拉公式有：

當(dāng)用點Xn處的率近似值K1與右端點Xn+1處的斜率K2的算術(shù)平均值作為平均斜率K*的近似值，那么就會得到二階精度的改進(jìn)歐拉公式：

依次類推，如果在區(qū)間[Xn，Xn+1]內(nèi)多預(yù)估幾個點上的斜率值K1、K2、……Km，并用他們的加權(quán)平均數(shù)作為平均斜率K*的近似值，顯然能構(gòu)造出具有很高精度的高階計算公式。經(jīng)數(shù)學(xué)推導(dǎo)、求解，可以得出四階龍格－庫塔公式，也就是在工程中應(yīng)用廣泛的經(jīng)典四階龍格－庫塔算法[10]：

3 問題描述

如圖1所示的一段被嵌入墻內(nèi)的懸臂梁，在固定點A，位移y和斜度dz皆為零。設(shè)梁是一根均勻細(xì)桿，其長度為L。梁關(guān)于垂直方向y的擾度滿足微分方程：

其中：I是梁的橫截面關(guān)于其主軸的慣性力矩；E是彈性模量；常數(shù)ρ是梁的線密度；g是重力加速度。

選取參數(shù)L=2m，線密度ρ=10kg/m，慣性力矩與彈性模量的乘積 IE=2400kg·m3/s2[11]。

根據(jù)以上條件，確定出懸臂梁中擾度最大的位置以及求出相對應(yīng)的最大擾度，作出懸臂梁擾度圖（z-y圖）。

圖1 懸臂梁示意圖

4 MATLAB近似求解

將區(qū)間[0,2]進(jìn)行100等分，采用經(jīng)典四階RK公式計算。

用RK求解微分方程的MATLAB程序[12]如下：function RK4=RK4(a,b,m)

運行MATLAB程序得出懸臂梁擾度圖（z-y圖）如圖2所示。

圖2 懸臂梁繞度圖

再運行程序得出近似解：最大擾度在端點B出現(xiàn)，即ymax1（2）=16.48cm。

5 解析法精確求解

將已知參數(shù)L=2m，線密度ρ=10kg/m，慣性力矩與彈性模量的乘積 IE=2400kg·m3/s2代入（11）式，化簡求解得：

對（15）式求最值得出：當(dāng)z=2時，即在懸臂梁最右端，撓度最大，ymax2（2）=16.67cm。

6 結(jié)論