《機(jī)器人技術(shù)-建模、仿真及應(yīng)用》課件 第四章 機(jī)器人動(dòng)力學(xué)

ROBOT機(jī)器人技術(shù)——建模、仿真及應(yīng)用機(jī)器人動(dòng)力學(xué)第四章

目錄剛體動(dòng)力學(xué)基礎(chǔ)PART.1牛頓-歐拉迭代動(dòng)力學(xué)方程PART.2歐拉-拉格朗日方程PART.3引言

動(dòng)力學(xué)正問(wèn)題是已知機(jī)械臂各關(guān)節(jié)的作用力或力矩，求各關(guān)節(jié)的位移、速度和加速度，即機(jī)器人的運(yùn)動(dòng)軌跡

，這可以用于對(duì)機(jī)械臂的仿真。

動(dòng)力學(xué)逆問(wèn)題是已知機(jī)械臂的運(yùn)動(dòng)軌跡，即各關(guān)節(jié)的位移、速度和加速度，求各關(guān)節(jié)所需要的驅(qū)動(dòng)力或力矩

，這可以用于對(duì)機(jī)械臂的控制。

動(dòng)力學(xué)主要研究產(chǎn)生運(yùn)動(dòng)所需要的力。對(duì)于機(jī)器人動(dòng)力學(xué)分析，有兩種經(jīng)典的方法：一種是牛頓-歐拉法，另一種是拉格朗日法。與機(jī)器人運(yùn)動(dòng)學(xué)相似，機(jī)器人動(dòng)力學(xué)也有兩個(gè)相反的問(wèn)題：剛體動(dòng)力學(xué)基礎(chǔ)PART.1質(zhì)量分布剛體的加速度剛體動(dòng)力學(xué)基礎(chǔ)質(zhì)量分布如圖表示一個(gè)剛體，坐標(biāo)系建立在剛體上，

表示單元體

的位置矢量。坐標(biāo)系

中的慣性張量可用3×3矩陣表示如下：矩陣中的各元素為：式中剛體由單元體

組成，單元體的密度為

。每個(gè)單元體的位置由矢量

確定，如圖所示。慣性張量剛體動(dòng)力學(xué)基礎(chǔ)質(zhì)量分布

,

和

為慣量矩，是單元體質(zhì)量

乘以單元體到相應(yīng)轉(zhuǎn)軸垂直距離的平方在整個(gè)剛體上的積分。

交叉項(xiàng)稱為慣量積。對(duì)于一個(gè)剛體來(lái)說(shuō)，上述六個(gè)相互獨(dú)立參量取決于所在坐標(biāo)系的位姿。慣性張量剛體動(dòng)力學(xué)基礎(chǔ)剛體的加速度有關(guān)剛體加速度問(wèn)題：任一瞬時(shí)，對(duì)剛體的線速度和角速度求導(dǎo)。例如，加速度可以通過(guò)計(jì)算空間一點(diǎn)

相對(duì)于坐標(biāo)系

的速度的微分進(jìn)行描述，即：同速度一樣，當(dāng)微分的參考坐標(biāo)系為世界坐標(biāo)系時(shí)，可用下列符號(hào)表示剛體的速度，即：線速度剛體動(dòng)力學(xué)基礎(chǔ)剛體的加速度把坐標(biāo)系

固連在一剛體上，要求描述相對(duì)于坐標(biāo)系

的速度矢量，如圖，這里假設(shè)坐標(biāo)系

是固定的。坐標(biāo)系

相對(duì)于坐標(biāo)系

的位置矢量

和旋轉(zhuǎn)矩陣

來(lái)描述，假設(shè)方位

不隨時(shí)間變化，則

點(diǎn)相對(duì)于坐標(biāo)

的運(yùn)動(dòng)是由于

或

隨時(shí)間的變化引起的。求解坐標(biāo)系

中

的線速度只要寫出坐標(biāo)系

中的兩個(gè)速度分量，求其和為：公式（4.7）只適用于坐標(biāo)系

和坐標(biāo)系

的相對(duì)方位保持不變的情況。坐標(biāo)系

以速度

相對(duì)于坐標(biāo)系

平移坐標(biāo)系

相對(duì)于坐標(biāo)系

的方位是隨時(shí)間變化的，

相對(duì)于

的旋轉(zhuǎn)速度用矢量

來(lái)表示。已知矢量

確定了坐標(biāo)系

中一固定點(diǎn)的位置，則可得點(diǎn)

的角速度為：角速度剛體動(dòng)力學(xué)基礎(chǔ)剛體的加速度固定坐標(biāo)系

中的矢量

以角速度

相對(duì)于坐標(biāo)系

旋轉(zhuǎn)首先討論兩坐標(biāo)系的原點(diǎn)重合、相對(duì)相速度為零的情況，而且它們的原點(diǎn)始終保持重合，其中一個(gè)或兩個(gè)坐標(biāo)系固連在剛體上，如圖所示。

線速度和角速度同時(shí)存在時(shí)，且兩坐標(biāo)系原點(diǎn)不重合，把線速度帶入上式，可以得到從坐標(biāo)系

觀測(cè)坐標(biāo)系

中固定速度矢量的普遍公式：線加速度剛體動(dòng)力學(xué)基礎(chǔ)剛體的加速度式

描述了當(dāng)坐標(biāo)系

的原點(diǎn)與坐標(biāo)系

的原點(diǎn)重合時(shí)，坐標(biāo)系

下的速度矢量

，方程左邊描述的是矢量

隨時(shí)間變化的情況，由于兩個(gè)坐標(biāo)系的原點(diǎn)重合，因此可以把改寫成如下形式：對(duì)式

求導(dǎo)，當(dāng)坐標(biāo)系

的原點(diǎn)與坐標(biāo)系

的原點(diǎn)重合時(shí)，可得到

的加速度在坐標(biāo)系

的表達(dá)式：上式第一項(xiàng)和最后一項(xiàng)應(yīng)用式（1），則式（2）變?yōu)椋海?）（2）線加速度剛體動(dòng)力學(xué)基礎(chǔ)剛體的加速度將

同類項(xiàng)合并，整理得：為了將結(jié)論推廣到兩個(gè)坐標(biāo)系原點(diǎn)不重合的一般情況，附加一個(gè)表示坐標(biāo)系

原點(diǎn)線加速度的項(xiàng)，最終得到一般表達(dá)式：值得指出的是，當(dāng)

是常量時(shí)，即：此時(shí)，上式簡(jiǎn)化為：角加速度剛體動(dòng)力學(xué)基礎(chǔ)剛體的加速度假設(shè)坐標(biāo)系

以角速度

相對(duì)于坐標(biāo)系

轉(zhuǎn)動(dòng)，同時(shí)坐標(biāo)系

以角速度

相對(duì)于坐標(biāo)系

轉(zhuǎn)動(dòng)。為求

在坐標(biāo)系

中進(jìn)行矢量相加，即：對(duì)上式求導(dǎo)，得：將式

代入上式右側(cè)最后一項(xiàng)中，得：上式用于計(jì)算操作臂連桿的角加速度。牛頓-歐拉迭代動(dòng)力學(xué)方程PART.2歐拉方程牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓方程牛頓-歐拉法應(yīng)用實(shí)例牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓方程作用于剛體質(zhì)心的力F引起剛體運(yùn)動(dòng)加速度如圖所示的剛體質(zhì)心正以加速度

做加速運(yùn)動(dòng)。此時(shí)，由牛頓方程可得作用在質(zhì)心上的力F引起剛體加速度為：式中，m代表剛體總質(zhì)量。牛頓-歐拉迭代動(dòng)力學(xué)方程歐拉方程如圖所示為一個(gè)旋轉(zhuǎn)剛體，其角速度和角加速度分別為

、

。此時(shí)，由歐拉方程可得，作用在剛體上的力矩N引起剛體的轉(zhuǎn)動(dòng)為：作用在剛體上的力矩N，剛體旋轉(zhuǎn)角速度

和角加速度式中

是剛體在坐標(biāo)系{C}中的慣性張量。剛體的質(zhì)心在坐標(biāo)系{C}的原點(diǎn)上。牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程有關(guān)操作臂給定運(yùn)動(dòng)軌跡的力矩計(jì)算的問(wèn)題。（1）計(jì)算速度和加速度的向外迭代法：為了計(jì)算作用在連桿上的慣性力，需要計(jì)算操作臂每個(gè)連桿在某一時(shí)刻的角速度、線加速度和角加速度?？蓱?yīng)用迭代方法完成這些計(jì)算。首先對(duì)連桿1進(jìn)行計(jì)算，接著計(jì)算下一個(gè)連桿，這樣一直向外迭代到連桿n。角速度在連桿之間的傳遞如圖所示，連桿

的角速度為：由角加速度的公式可得：牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程當(dāng)?shù)趇+1個(gè)關(guān)節(jié)是移動(dòng)關(guān)節(jié)時(shí)，上式可簡(jiǎn)化為：由線速度公式可以得到每個(gè)連桿坐標(biāo)系原點(diǎn)的線加速度：當(dāng)?shù)趇+1個(gè)關(guān)節(jié)是移動(dòng)關(guān)節(jié)時(shí)，上式可簡(jiǎn)化為：牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程進(jìn)而可以得到每個(gè)連桿質(zhì)心的線加速度：假定坐標(biāo)系{C}連于連桿i上,坐標(biāo)系原點(diǎn)位于連桿質(zhì)心，且各坐標(biāo)軸方位與原連桿坐標(biāo)系{i}方位相同。由于上式與關(guān)節(jié)的運(yùn)動(dòng)無(wú)關(guān)，因此無(wú)論是旋轉(zhuǎn)關(guān)節(jié)還是移動(dòng)關(guān)節(jié)，上式對(duì)于第i+1個(gè)連桿來(lái)說(shuō)都是有效的。計(jì)算每個(gè)連桿質(zhì)心的線加速度和角加速度后，運(yùn)用牛頓-歐拉公式可以分別計(jì)算出作用在連桿質(zhì)心上的慣性力和力矩，即：式中，坐標(biāo)系

的原點(diǎn)位于連桿質(zhì)心，各坐標(biāo)軸方位與原連桿坐標(biāo)系

方位相同。牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程（2）計(jì)算力和力矩的向內(nèi)迭代法：計(jì)算出作用在每個(gè)連桿上的力和力矩之后，需要計(jì)算關(guān)節(jié)力矩，它們是實(shí)際施加在連桿上的力和力矩。根據(jù)典型連桿在無(wú)重力狀態(tài)下的受力圖，如圖所示，列出力平衡方程和力矩平衡方程，每個(gè)連桿都受到相鄰連桿的作用力和作用力矩以及附加的。這里定義了一些專用符號(hào)用來(lái)表示相鄰的作用力和力矩：?jiǎn)蝹€(gè)連桿的力平衡、力矩平衡=連桿

作用在連桿

上的力；②=連桿

作用在連桿

上的力矩。將所有作用在連桿

上的力相加，得到力平衡方程如下：牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程將所有作用在質(zhì)心上的力矩相加，并且令他們的和為零，得到力矩平衡方程如下：將式

的結(jié)果以及附加旋轉(zhuǎn)矩陣的方法帶入上式，可得：最后，重新排列力和力矩方程，形成相鄰連桿從高序號(hào)向低序號(hào)排列的迭代關(guān)系分別為：牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程應(yīng)用這些方程對(duì)連桿依次求解，從連桿n開始向內(nèi)迭代一直到機(jī)器人基座。在靜力學(xué)中，可通過(guò)下式計(jì)算一個(gè)連桿施加于相鄰連桿的力矩在

方向的分量求得關(guān)節(jié)力矩：式中，

表示線性驅(qū)動(dòng)力。注意：對(duì)一個(gè)在自由空間中運(yùn)動(dòng)的機(jī)器人來(lái)說(shuō)，

和

等于零，因此應(yīng)用這些方程首先計(jì)算連桿n時(shí)是很簡(jiǎn)單的；如果機(jī)器人與環(huán)境接觸，

和

不為零，力平衡方程中就包含了接觸力和力矩。牛頓-歐拉迭代動(dòng)力學(xué)方程牛頓-歐拉法迭代動(dòng)力學(xué)方程（3）牛頓-歐拉迭代動(dòng)力學(xué)算法:由關(guān)節(jié)運(yùn)動(dòng)計(jì)算關(guān)節(jié)力矩的完整算法由兩部分組成:第一部分是對(duì)每個(gè)連桿應(yīng)用牛頓-歐拉方程，從連桿1到連桿n向外迭代計(jì)算連桿的速度和加速度；第二部分是對(duì)每個(gè)連桿n到連桿1向內(nèi)迭代計(jì)算連桿間的相互作用力和力矩以及關(guān)節(jié)驅(qū)動(dòng)力矩。對(duì)于轉(zhuǎn)動(dòng)關(guān)節(jié)來(lái)說(shuō)，這個(gè)算法歸納如下：1）外推，

向外迭代。2）內(nèi)推，

向內(nèi)迭代。已知關(guān)節(jié)位置、速度和加速度，應(yīng)用外推和內(nèi)推公式可以計(jì)算出所需的關(guān)節(jié)力矩。假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例牛頓-歐拉迭代動(dòng)力學(xué)方程首先，確定牛頓-歐拉迭代公式中各參量的值。每個(gè)連桿質(zhì)心的位置矢量分別為:由于假設(shè)為集中質(zhì)量，因此每個(gè)連桿質(zhì)心的慣性張量為零矩陣，即:末端執(zhí)行器上沒(méi)有作用力，因而有:機(jī)器人基座不旋轉(zhuǎn)，因此有:包括重力因素，有:相鄰連桿坐標(biāo)系之間的相對(duì)轉(zhuǎn)動(dòng)由下式給出:牛頓-歐拉法應(yīng)用實(shí)例應(yīng)用方程式外推和內(nèi)推公式，對(duì)連桿1用向外迭代法求解如下:牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例對(duì)連桿2用向外迭代法求解如下:牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例對(duì)連桿2用向內(nèi)迭代法求解如下:（48）牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例對(duì)連桿1用向內(nèi)迭代法求解如下:牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例牛頓-歐拉法應(yīng)用實(shí)例取中的方向分量，得關(guān)節(jié)力矩:上式將驅(qū)動(dòng)力矩表示為關(guān)于關(guān)節(jié)位置、速度和加速度的函數(shù)。牛頓-歐拉迭代動(dòng)力學(xué)方程假設(shè)機(jī)械臂每個(gè)連桿的質(zhì)量都集中在連桿的末端，設(shè)其質(zhì)量分別為

和

；設(shè)其桿長(zhǎng)分別為

、

，關(guān)節(jié)角分別為

、

，關(guān)節(jié)力矩分別為

、

；例歐拉-拉格朗日方程PART.3歐拉-拉格朗日方程歐拉-拉格朗日方程是用廣義坐標(biāo)表示完整工業(yè)機(jī)器人系統(tǒng)的動(dòng)力學(xué)方程，拉格朗日函數(shù)

定義為系統(tǒng)的全部動(dòng)能

和全部勢(shì)能

之差，即：式中，動(dòng)能

取決于機(jī)器人系統(tǒng)中連桿的位姿和速度；而勢(shì)能

取決于連桿的位形。系統(tǒng)動(dòng)力學(xué)方程式，即歐拉-拉格朗日方程如下：式中，

表示與廣義坐標(biāo)

相關(guān)的廣義力；

為相應(yīng)的廣義速度。勢(shì)能計(jì)算運(yùn)動(dòng)方程動(dòng)能計(jì)算拉格朗日法仿真實(shí)例歐拉-拉格朗日方程動(dòng)能計(jì)算在機(jī)械臂中，連桿是運(yùn)動(dòng)部件，連桿

的動(dòng)能

連桿質(zhì)心線速度產(chǎn)生的動(dòng)能和連桿角速度產(chǎn)生的動(dòng)能之和。因此，對(duì)于有

個(gè)連桿的機(jī)器人系統(tǒng)，其總動(dòng)能是每一連桿相關(guān)運(yùn)動(dòng)產(chǎn)生的動(dòng)能之和，表達(dá)如下：式中，

為連桿

的質(zhì)心

的三維線速度矢量；

為連桿

的三維角速度矢量；

為連桿

的質(zhì)量，是標(biāo)量；

為連桿

的慣性張量。由于

和

分別為關(guān)節(jié)變量

和關(guān)節(jié)速度

的函數(shù)，由

可知機(jī)器人的動(dòng)能是關(guān)節(jié)變量

和關(guān)節(jié)速度

的函數(shù)。歐拉-拉格朗日方程勢(shì)能計(jì)算與動(dòng)能計(jì)算類似，機(jī)器人的總勢(shì)能也是各連桿的勢(shì)能之和。假設(shè)連桿是剛性體，勢(shì)能的表達(dá)如下：式中，矢量

為關(guān)節(jié)變量的函數(shù)，且該函數(shù)是非線性的。由此可知，總勢(shì)能

只關(guān)于關(guān)節(jié)變量

的函數(shù)，與關(guān)節(jié)速度

無(wú)關(guān)。歐拉-拉格朗日方程運(yùn)動(dòng)方程按照式

和式

計(jì)算系統(tǒng)總動(dòng)能和總勢(shì)能，牛頓歐拉章節(jié)的例題中機(jī)器人的拉格朗日函數(shù)可寫為：拉格朗日函數(shù)對(duì)關(guān)節(jié)變量

、關(guān)節(jié)速度

和時(shí)間

求導(dǎo)可得動(dòng)力學(xué)運(yùn)動(dòng)方程，其中勢(shì)能與關(guān)節(jié)速度

無(wú)關(guān)，即有：式中，

表示與廣義坐標(biāo)

相關(guān)的廣義力；

為相應(yīng)的廣義速度。實(shí)際計(jì)算機(jī)械臂動(dòng)力學(xué)時(shí)，

和

對(duì)應(yīng)連桿轉(zhuǎn)矩

；而

和

對(duì)應(yīng)連桿推力

。歐拉-拉格朗日方程拉格朗日法仿真實(shí)例機(jī)器人是結(jié)構(gòu)復(fù)雜的連桿系統(tǒng)，一般采用齊次變換的方法，用拉格朗日方程建立其系統(tǒng)動(dòng)力學(xué)方程，對(duì)其位姿和運(yùn)動(dòng)狀態(tài)進(jìn)行描述。機(jī)器人動(dòng)力學(xué)方程的具體推導(dǎo)過(guò)程如下：1.選取坐標(biāo)系，選定完全而且獨(dú)立的廣義關(guān)節(jié)變量

。2.選定相應(yīng)關(guān)節(jié)上的廣義力

：當(dāng)

位移變量時(shí)，

為力；當(dāng)

是角度變量時(shí)，

為力矩。3.求出機(jī)器人各構(gòu)件的動(dòng)能和勢(shì)能，構(gòu)造拉格朗日函數(shù)。4.代入拉格朗日方程求得機(jī)器人系統(tǒng)動(dòng)力學(xué)方程。例歐拉-拉格朗日方程拉格朗日法仿真實(shí)例說(shuō)明下圖機(jī)器人動(dòng)力學(xué)方程的推導(dǎo)過(guò)程。選取笛卡兒標(biāo)系。連桿1和連桿2的關(guān)節(jié)變量分別是轉(zhuǎn)角

和

，連桿1和連桿2的質(zhì)量分別是

和

，桿長(zhǎng)分別為

和

，質(zhì)心分別在

和

處，離關(guān)節(jié)中心的距離分別為

和

，其中底座與大地固定。因此，連桿1質(zhì)心

的位置坐標(biāo)為：連桿1質(zhì)心

速度的平方為：連桿2質(zhì)心

的位置坐標(biāo)為：連桿2質(zhì)心

速度的平方為：歐拉-拉格朗日方程拉格朗日法仿真實(shí)例系統(tǒng)動(dòng)能為：系統(tǒng)勢(shì)能為：拉格朗日函數(shù)為：選取笛卡兒標(biāo)系，設(shè)定連桿1的關(guān)節(jié)變量

，則其角速度連桿2的關(guān)節(jié)變量

，則其角速度

，連桿1的桿長(zhǎng)

，其質(zhì)量

。連桿2的桿長(zhǎng)

，其質(zhì)量

。根據(jù)上述動(dòng)力學(xué)方程的推導(dǎo)過(guò)程對(duì)機(jī)械臂10s內(nèi)的動(dòng)能和勢(shì)能變化進(jìn)行仿真求解，仿真代碼及結(jié)果如下。歐拉-拉格朗日方程拉格朗日法仿真實(shí)例例已知條件：

歐拉-拉格朗日方程拉格朗日法仿真實(shí)例t=0:2:10;theta1=2*pi/180*(t.^2);theta2=pi/180*(t.^2);a=4*pi/180*t;b=2*pi/180*t;m1=20;m2=15;l1=4;l2=3;g=9.8;p1=1/2*l1;p2=1/2*l2;X1=p1*sin(theta1);Y1=-p1*cos(theta1);V12=(p1*a).^2;X2=l1*sin(theta1)+p2*sin(theta1+theta2);Y2=-l1*cos(theta1)-p2*cos(theta1+theta2);V22=(l1^2)*(a.^2)+(p2^2)*((a+b).^2)+2*l1*p2*((a.^2)+a.*b).*cos(theta2);Ek1=1/2*m1*V12;Ek2=1/2*m2*V22;Ek=Ek1+Ek2;subplot(2,1,1);plot(t