系統(tǒng)建模與仿真習題答案(forstudents)

第一章習題1-1什么是仿真？它所遵循的根本原則是什么？答：仿真是建立在控制理論，相似理論，信息處理技術(shù)和計算技術(shù)等理論根底之上的，以計算機和其他專用物理效應設(shè)備為工具，利用系統(tǒng)模型對真實或假想的系統(tǒng)進行試驗，并借助專家經(jīng)驗知識，統(tǒng)計數(shù)據(jù)和信息資料對試驗結(jié)果進行分析和研究，進而做出決策的一門綜合性的試驗性科學。它所遵循的根本原則是相似原理。1-2在系統(tǒng)分析與設(shè)計中仿真法與解析法有何區(qū)別？各有什么特點？答：解析法就是運用已掌握的理論知識對控制系統(tǒng)進行理論上的分析，計算。它是一種純物理意義上的實驗分析方法，在對系統(tǒng)的認識過程中具有普遍意義。由于受到理論的不完善性以及對事物認識的不全面性等因素的影響，其應用往往有很大局限性。仿真法基于相似原理，是在模型上所進行的系統(tǒng)性能分析與研究的實驗方法。1-3數(shù)字仿真包括那幾個要素？其關(guān)系如何？答:通常情況下，數(shù)字仿真實驗包括三個根本要素，即實際系統(tǒng)，數(shù)學模型與計算機。由圖可見，將實際系統(tǒng)抽象為數(shù)學模型，稱之為一次模型化，它還涉及到系統(tǒng)辨識技術(shù)問題，統(tǒng)稱為建模問題；將數(shù)學模型轉(zhuǎn)化為可在計算機上運行的仿真模型，稱之為二次模型化，這涉及到仿真技術(shù)問題，統(tǒng)稱為仿真實驗。1-4為什么說模擬仿真較數(shù)字仿真精度低？其優(yōu)點如何？。答：由于受到電路元件精度的制約和容易受到外界的干擾，模擬仿真較數(shù)字仿真精度低但模擬仿真具有如下優(yōu)點：描述連續(xù)的物理系統(tǒng)的動態(tài)過程比擬自然和逼真。仿真速度極快，失真小，結(jié)果可信度高。能快速求解微分方程。模擬計算機運行時各運算器是并行工作的，模擬機的解題速度與原系統(tǒng)的復雜程度無關(guān)?？梢造`活設(shè)置仿真試驗的時間標尺，既可以進行實時仿真，也可以進行非實時仿真。易于和實物相連。1-5什么是CAD技術(shù)？控制系統(tǒng)CAD可解決那些問題？答：CAD技術(shù)，即計算機輔助設(shè)計〔ComputerAidedDesign〕,是將計算機高速而精確的計算能力，大容量存儲和處理數(shù)據(jù)的能力與設(shè)計者的綜合分析，邏輯判斷以及創(chuàng)造性思維結(jié)合起來，用以加快設(shè)計進程，縮短設(shè)計周期，提高設(shè)計質(zhì)量的技術(shù)。控制系統(tǒng)CAD可以解決以頻域法為主要內(nèi)容的經(jīng)典控制理論和以時域法為主要內(nèi)容的現(xiàn)代控制理論。此外，自適應控制，自校正控制以及最優(yōu)控制等現(xiàn)代控制測策略都可利用CAD技術(shù)實現(xiàn)有效的分析與設(shè)計。1-6什么是虛擬現(xiàn)實技術(shù)？它與仿真技術(shù)的關(guān)系如何？答：虛擬現(xiàn)實技術(shù)是一種綜合了計算機圖形技術(shù)，多媒體技術(shù)，傳感器技術(shù)，顯示技術(shù)以及仿真技術(shù)等多種學科而開展起來的高新技術(shù)。虛擬現(xiàn)實技術(shù)不斷完善，為控制系統(tǒng)數(shù)字仿真與CAD開辟了一個新時代。1-7什么是離散系統(tǒng)？什么是離散事件系統(tǒng)？如何用數(shù)學的方法描述它們？答：本書所講的“離散系統(tǒng)〞指的是離散時間系統(tǒng)，即系統(tǒng)中狀態(tài)變量的變化僅發(fā)生在一組離散時刻上的系統(tǒng)。它一般采用差分方程，離散狀態(tài)方程和脈沖傳遞函數(shù)來描述。離散事件系統(tǒng)是系統(tǒng)中狀態(tài)變量的改變是由離散時刻上所發(fā)生的事件所驅(qū)動的系統(tǒng)。這種系統(tǒng)的輸入輸出是隨機發(fā)生的，一般采用概率模型來描述。1-8如圖1-16所示某衛(wèi)星姿態(tài)控制仿真實驗系統(tǒng)，試說明：假設(shè)按模型分類，該系統(tǒng)屬于那一類仿真系統(tǒng)？圖中“混合計算機〞局部在系統(tǒng)中起什么作用？與數(shù)字仿真相比該系統(tǒng)有什么優(yōu)缺點？答：〔1〕按模型分類，該系統(tǒng)屬于物理仿真系統(tǒng)。〔2〕混合計算機集中了模擬仿真和數(shù)字仿真的優(yōu)點，它既可以與實物連接進行實時仿真，計算一些復雜函數(shù)，又可以對控制系統(tǒng)進行反復迭代計算。其數(shù)字局部用來模擬系統(tǒng)中的控制器，而模擬局部用于模擬控制對象。與數(shù)字仿真相比，物理仿真總是有實物介入，效果逼真，精度高，具有實時性與在線性的特點，但其構(gòu)成復雜，造價較高，耗時過長，通用性不強。第二章習題2-1思考題：〔1〕數(shù)學模型的微分方程，狀態(tài)方程，傳遞函數(shù)，零極點增益和局部分式五種形式，各有什么特點？答：微分方程是直接描述系統(tǒng)輸入和輸出量之間的制約關(guān)系，是連續(xù)控制系統(tǒng)其他數(shù)學模型表達式的根底。狀態(tài)方程能夠反映系統(tǒng)內(nèi)部各狀態(tài)之間的相互關(guān)系，適用于多輸入多輸出系統(tǒng)。傳遞函數(shù)是零極點形式和局部分式形式的根底。零極點增益形式可用于分析系統(tǒng)的穩(wěn)定性和快速性。利用局部分式形式可直接分析系統(tǒng)的動態(tài)過程?！?〕數(shù)學模型各種形式之間為什么要互相轉(zhuǎn)換？答：不同的控制系統(tǒng)的分析和設(shè)計方法，只適用于特定的數(shù)學模型形式?！?〕控制系統(tǒng)建模的根本方法有哪些？他們的區(qū)別和特點是什么？答：控制系統(tǒng)的建模方法大體有三種：機理建模法，實驗建模法和綜合建模法。機理建模法就是對結(jié)構(gòu)，參數(shù)的物理系統(tǒng)運用相應的物理定律或定理，經(jīng)過合理的分析簡化建立起來的各物理量間的關(guān)系。該方法需要對系統(tǒng)的內(nèi)部結(jié)構(gòu)和特性完全的了解，精度高。實驗建模法是采用歸納的方法，根據(jù)系統(tǒng)實測的數(shù)據(jù)，運用統(tǒng)計規(guī)律和系統(tǒng)辨識等理論建立的系統(tǒng)模型。該方法建立的數(shù)學模型受數(shù)據(jù)量不充分，數(shù)據(jù)精度不一致，數(shù)據(jù)處理方法的不完善，很難在精度上到達更高的要求。綜合建模法是上述兩種方法的結(jié)合。〔4〕控制系統(tǒng)計算機仿真中的“實現(xiàn)問題〞是什么含意？答：“實現(xiàn)問題〞就是根據(jù)建立的數(shù)學模型和精度，采用某種數(shù)值計算方法，將模型方程轉(zhuǎn)換為適合在計算機上運行的公式和方程，通過計算來使之正確的反映系統(tǒng)各變量動態(tài)性能，得到可靠的仿真結(jié)果。〔5〕數(shù)值積分法的選用應遵循哪幾條原則？答：數(shù)值積分法應該遵循的原則是在滿足系統(tǒng)精度的前提下，提高數(shù)值運算的速度和并保證計算結(jié)果的穩(wěn)定。2-2.用matlab語言求以下系統(tǒng)的狀態(tài)方程、傳遞函數(shù)、零極點增益、和局部分式形式的模型參數(shù)，并分別寫出其相應的數(shù)學模型表達式：(1)G〔s〕=(2)MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\h=y=[0202]X解：〔1〕狀態(tài)方程模型參數(shù):編寫matlab程序如下>>num=[172424];>>den=[110355024];>>[ABCD]=tf2ss(num,den)得到結(jié)果：A=,B=,C=,D=[0]所以模型為:=X+u,y=X(2)零極點增益:編寫程序>>num=[172424];>>den=[110355024];>>[ZPK]=tf2zp(num,den)得到結(jié)果Z=-2.7306+2.8531i,-2.7306-2.8531i,-1.5388P=-4,-3,-2,-1K=1(3)局部分式形式:編寫程序>>num=[172424];>>den=[110355024];>>[RPH]=residue(num,den)得到結(jié)果R=4.0000,-6.0000,2.0000,1.0000P=-4.0000,-3.0000,-2.0000,-1.0000H=[]G(s)=〔2〕解：〔1〕傳遞函數(shù)模型參數(shù)：編寫程序>>A=[2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75]；>>B=[4220]';>>C=[0202];>>D=[0];>>[numden]=ss2tf(A,B,C,D)得到結(jié)果num=04.000014.000022.000015.0000den=1.00004.00006.25005.25002.2500(2)零極點增益模型參數(shù):編寫程序>>A=[2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75]；>>B=[4220]';>>C=[0202];>>D=[0];>>[Z,P,K]=ss2zp(A,B,C,D)得到結(jié)果Z=-1.0000+1.2247i-1.0000-1.2247i-1.5000P=-0.5000+0.8660i-0.5000-0.8660i-1.5000-1.5000K=4.0000表達式〔3〕局部分式形式的模型參數(shù)：編寫程序>>A=[2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75]；>>B=[4220]';>>C=[0202];>>D=[0];>>[numden]=ss2tf(A,B,C,D)>>[R,P,H]=residue(num,den)得到結(jié)果R=4.0000-0.00000.0000-2.3094i0.0000+2.3094iP=-1.5000-1.5000-0.5000+0.8660i-0.5000-0.8660iH=[]2-3.用歐拉法求下面系統(tǒng)的輸出響應y(t)在0≤t≤1上，h=0.1時的數(shù)值。要求保存4位小數(shù)，并將結(jié)果與真解比擬。解：歐拉法〔前向歐拉法,可以自啟動〕其幾何意義：把f(t,y)在[]區(qū)間內(nèi)的曲邊面積用矩形面積近似代替。利用matlab提供的m文件編程，得到算法公式。如下所示〔1〕m文件程序為h=0.1;disp('函數(shù)的數(shù)值解為');%顯示‘’中間的文字%disp('y=');%同上%y=1;fort=0:h:1m=y;disp(y);%顯示y的當前值%y=m-m*h;end保存文件q2.m在matalb命令行中鍵入>>q2得到結(jié)果函數(shù)的數(shù)值解為y=10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487〔2〕另建一個m文件求解在t[0,1]的數(shù)值〔%是的真解%〕程序為h=0.1;disp('函數(shù)的離散時刻解為');disp('y=');fort=0:h:1y=exp(-t);disp(y);end保存文件q3.m在matalb命令行中鍵入>>q3函數(shù)的離散時刻解為y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679比擬歐拉方法求解與真值的差異歐拉10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679誤差0-0.0048-0.0007–0.0118–0.0142–0.0160–0.0174–0.0183–0.0188-0.0192-0.0192顯然誤差與為同階無窮小，歐拉法具有一階計算精度，精度較低，但算法簡單。2-4用二階龍格庫塔法求解2-3的數(shù)值解，并于歐拉法求得的結(jié)果比擬。解：我們經(jīng)常用到預報-校正法的二階龍-格庫塔法，此方法可以自啟動，具有二階計算精度，幾何意義：把f(t,y)在[]區(qū)間內(nèi)的曲邊面積用上下底為和、高為h的梯形面積近似代替。利用matlab提供的m文件編程，得到算法公式。如下所示〔1〕m文件程序為h=0.1;disp('函數(shù)的數(shù)值解為');disp('y=');y=1;fort=0:h:1disp(y);k1=-y;k2=-(y+k1*h);y=y+(k1+k2)*h/2;end保存文件q4.m在matlab的命令行中鍵入>>q4顯示結(jié)果為函數(shù)的數(shù)值解為y=10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685比擬歐拉法與二階龍格-庫塔法求解.〔誤差為絕對值〕真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679龍庫10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685誤差00.00020.00030.00040.00050.00060.00060.00060.00070.00060.0006明顯誤差為得同階無窮小，具有二階計算精度，而歐拉法具有以階計算精度，二階龍格-庫塔法比歐拉法計算精度高。2-5．用四階龍格-庫塔法求解題2-3數(shù)值解，并與前兩題結(jié)果相比擬。解：四階龍格-庫塔法表達式，其截斷誤差為同階無窮小，當h步距取得較小時，誤差是很小的.編輯m文件程序h=0.1;disp('四階龍格-庫塔方法求解函數(shù)數(shù)值解為');disp('y=');y=1;fort=0:h:1disp(y=);k1=-y;k2=-(y+k1*h/2);k3=-(y+k2*h/2);k4=-(y+k3*h);y=y+(k1+2*k2+2*k3+k4)*h/6;end保存文件q5.m在matlab命令行里鍵入>>q5得到結(jié)果四階龍格-庫塔方法求解函數(shù)數(shù)值解為y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679(2)比擬這幾種方法：對于四階龍格-庫塔方法真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679龍庫10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679誤差00000000000顯然四階龍格-庫塔法求解精度很高，根本接近真值。三種方法比擬可以得到精度〔四階〕〉精度〔二階〕〉精度〔歐拉〕2-6．二階系統(tǒng)狀態(tài)方程為寫出取計算步長為h時，該系統(tǒng)狀態(tài)變量X=[]的四階龍格-庫塔法遞推關(guān)系式。解：四階龍格-庫塔法表達式所以狀態(tài)變量的遞推公式可以寫作：A=,B=,可以寫成則遞推形式2-7單位反應系統(tǒng)的開環(huán)傳遞函數(shù)如下用matlab語句、函數(shù)求取系統(tǒng)閉環(huán)零極點，并求取系統(tǒng)閉環(huán)狀態(tài)方程的可控標準型實現(xiàn)。解：開環(huán)傳遞函數(shù)，求得閉環(huán)傳遞函數(shù)為在matlab命令行里鍵入>>a=[10];>>b=[14.6];>>c=[13.416.35];>>d=conv(a,b)；>>e=conv(d,c)e=1.00008.000031.990075.21000>>f=[0005100];>>g=e+fg=1.00008.000031.990080.2100100.0000%以上是計算閉環(huán)傳遞函數(shù)的特征多項式%>>p=roots(g)%計算特征多項式的根，就是閉環(huán)傳遞函數(shù)的極點%p=-0.9987+3.0091i-0.9987-3.0091i-3.0013+0.9697i-3.0013-0.9697i>>m=[5100];>>z=roots(m)z=-20%計算零點%綜上：當閉環(huán)傳函形如時，可控標準型為：；所以可控標準型是2-8用matlab語言編制單變量系統(tǒng)三階龍格-庫塔法求解程序，程序入口要求能接收狀態(tài)方程各系數(shù)陣〔A,B,C,D〕,和輸入階躍函數(shù)r(t)=R*1(t);程序出口應給出輸出量y〔t〕的動態(tài)響應數(shù)值解序列。解：m文件為：functiony=hs(A,B,C,D,R,T,h)%T為觀測時間,h為計算步長，R為輸入信號幅值%disp('數(shù)值解為');y=0;r=R;x=[0;0;0;0];N=T/h;fort=1:N;k1=A*x+B*R;k2=A*(x+h*k1/3)+B*R;k3=A*(x+2*h*k2/3)+B*R;x=x+h*(k1+3*k3)/4;y(t)=C*x+D*R;end在命令行里鍵入A=B=C=D=R=T=h=y=hs(A,B,C,D,R,T,h)得到結(jié)果。2-9．用題2-8仿真程序求解題2-7系統(tǒng)的閉環(huán)輸出響應y(t).解：A=,B=,C=,D=[0]在命令行里鍵入>>A=[010000100001-100-80.21-31.99-8];>>B=[0001]';>>C=[-100500];>>D=[0];>>T=1;>>R=1;>>h=0.01;>>y=hs(A,B,C,D,R,T,h)數(shù)值解為08.3333e-0075.8659e-0061.8115e-0053.9384e-0057.0346e-005。。。。%僅取一局部%2-10.用式〔2-34〕梯形法求解試驗方程，分析對計算步長h有何限制，說明h對數(shù)值穩(wěn)定性的影響。解：編寫梯形法程序為得到穩(wěn)定系統(tǒng)最終漸進收斂。系統(tǒng)穩(wěn)定則計算得。h的選取不能超出上述范圍，否則系統(tǒng)不穩(wěn)定。2-11如圖2-27所示斜梁滾球系統(tǒng)，假設(shè)要研究滾球在梁上的位置可控性，需首先建立其數(shù)學模型，力矩電機的輸出轉(zhuǎn)矩M與其電流i成正比，橫梁為均勻可自平衡梁〔即當電機不通電且無滾球時，橫梁可處于=0的水平狀態(tài)〕，是建立系統(tǒng)的數(shù)學模型，并給出簡化后系統(tǒng)的動態(tài)結(jié)構(gòu)圖。解：設(shè)球的質(zhì)心到桿的距離為0，該系統(tǒng)為特殊情況下的球棒系統(tǒng)。另令分別表示棒的慣量、球的質(zhì)量和球的慣量。則球質(zhì)心的位置和速度為其中，。因而動能的移動局部為因而動能的移動局部為球棒系統(tǒng)的旋轉(zhuǎn)動能為因而，系統(tǒng)總的動能等于其中為常數(shù)。此系統(tǒng)的拉格朗日方程組為綜合以上公式的系統(tǒng)的方程組為設(shè)系統(tǒng)在平衡點附近，，，則系統(tǒng)方程可化為對上式進行拉普拉斯變換并化簡后可得到。參考文獻：[1]Hauser,S.Sestry,andP.Kokotovic.“Nonlinearcontrolviaapproximateinput-outputlinearization〞.IEEETrans.onAutomaticControl,vol.37:pp.392-398,1992.[2]R.Sepulchre.“Slowpeakingandlow-gaindesignsforglobalstabilizationofnonlinearsystems〞.submittedforIEEETAC1999.[3]R.Sepulchre,M.Jankovic,andP.KokotovicConstructiveNonlinearControl.Springer-Verlag,1997.[4]R.Teel.“UsingSaturationtostabilizeaclassofsingle-inputpartiallylinearcompositesystems〞.IFACNOLCOS'92Symposium,pages369-374,June1992.2-12如圖2-28所示雙水箱系統(tǒng)中，為流入水箱1的液體流量，為流出水箱2的液體流量，試依據(jù)液容與液阻的概念，建立的系統(tǒng)動態(tài)結(jié)構(gòu)圖。解：根據(jù)液容和液阻的概念，可分別列出兩個水箱的數(shù)學模型對上式進行在零初始條件下進行拉普拉斯變換得化簡后可得第三章習題3-2設(shè)典型閉環(huán)結(jié)構(gòu)控制系統(tǒng)如圖4-47所示，當階躍輸入幅值時，用sp4_1.m求取輸出的響應。解：用sp4_1.m求解過程如下：在MATLAB語言環(huán)境下，輸入以下命令語句>>a=[0.0160.8643.273.421];>>b=[3025];>>X0=[0000];%系統(tǒng)狀態(tài)向量初值為零>>V=2;%反應系數(shù)>>n=4;>>T0=0;Tf=10;>>h=0.01;R=20;%仿真步長h=0.01，階躍輸入幅值>>sp4_1%調(diào)用sp4_1.m函數(shù)>>plot(t,y)運行結(jié)果為：附：sp4_1.m函數(shù)為b=b/a(1);a=a/a(1);A=a(2:n+1);A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];B=[zeros(1,n-1),1]';m1=length(b);C=[fliplr(b),zeros(1,n-m1)];Ab=A-B*C*V;X=X0';y=0;t=T0;N=round((Tf-T0)/h);fori=1:NK1=Ab*X+B*R;K2=Ab*(X+h*K1/2)+B*R;K3=Ab*(X+h*K2/2)+B*R;K4=Ab*(X+h*K3)+B*R;X=X+h*(K1+2*K2+2*K3+K4)/6;y=[y,C*X];t=[t,t(i)+h];end3-4系統(tǒng)結(jié)構(gòu)圖如圖3-55，寫出該系統(tǒng)的聯(lián)結(jié)矩陣和，并寫出聯(lián)結(jié)矩陣非零元素陣。解：根據(jù)圖3-55中,拓撲連結(jié)關(guān)系，可寫出每個環(huán)節(jié)輸入受哪些環(huán)節(jié)輸出的影響，現(xiàn)列如入下:把環(huán)節(jié)之間的關(guān)系和環(huán)節(jié)與參考輸入的關(guān)系分別用矩陣表示出來，即=，=，3-6假設(shè)系統(tǒng)為圖4-5b雙輸入-雙輸出結(jié)構(gòu)，試寫出該系統(tǒng)的聯(lián)接矩陣，，說明應注意什么？解：根據(jù)圖4-5b中,拓撲連結(jié)關(guān)系，可列寫如下關(guān)系式:轉(zhuǎn)換成矩陣形式為所以聯(lián)接矩陣=，=此時應注意輸入聯(lián)接矩陣變?yōu)樾汀?--8求圖3-56非線性系統(tǒng)的輸出響應y(t),并與無非線性環(huán)節(jié)情況進行比擬。解：〔1〕不考慮非線性環(huán)節(jié)影響時，求解過程如下：先將環(huán)節(jié)編號標入圖中。2)在MATLAB命令窗口下，按編號依次將環(huán)節(jié)參數(shù)輸入P陣；>>P=[0.110.51;01200;2110;10110];3)按各環(huán)節(jié)相對位置和聯(lián)接關(guān)系，有聯(lián)接矩陣如下：，，所以非零元素矩陣>>WIJ=[101;14-1;211;321;431];4〕由于不考慮非線性影響，則非線性標志向量和參數(shù)向量均應賦零值；>>Z=[0000];S=[0000];5〕輸入運行參數(shù)：開環(huán)截至頻率約為1，故計算步長h取經(jīng)驗公式值，即，取h=0.01；每0.25秒輸出一點。故取=25。>>h=0.01;>>L1=25;>>n=4;>>T0=0>>Tf=20;>>nout=4;>>Y0=10;>>sp4_4;>>plot(t,y,'r')>>holdon運行結(jié)果如圖中紅色實線所示。(2)考慮非線性環(huán)節(jié)N影響時，只需將非線性標志向量Z和參數(shù)向量S的相應分量正確輸入即可。在MATLAB命令窗口中輸入以下語句：>>Z=[4000];S=[5000];%第一個線性環(huán)節(jié)后有飽和非線性，參數(shù)值為5。>>sp4_4;>>plot(t,y,'--')運行結(jié)果如圖中藍色虛線所示。從圖中可以清楚的地看出，飽和非線性環(huán)節(jié)對線性系統(tǒng)輸出響應的影響。附：sp4_4函數(shù)為:A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1));W0=zeros(n,1);W=zeros(n,n);fork=1:mif(WIJ(k,2)==0);W0(WIJ(k,1))=WIJ(k,3);elseW(WIJ(k,1),WIJ(k,2))=WIJ(k,3);end;end;fori=1:nif(A(i)==0);FI(i)=1;FIM(i)=h*C(i)/B(i);FIJ(i)=h*h*(C(i)/B(i))/2;FIC(i)=1;FID(i)=0;if(D(i)~=0);FID(i)=D(i)/B(i);elseendelseFI(i)=exp(-h*A(i)/B(i));FIM(i)=(1-FI(i))*C(i)/A(i);FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i);FIC(i)=1;FID(i)=0;if(D(i)~=0);FIC(i)=C(i)/D(i)-A(i)/B(i);FID(i)=D(i)/B(i);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);fork=1:N-1fori=1:L1Ub=Uk;Uk=W*Y+W0*Y0;fori=1:nif(Z(i)~=0)if(Z(i)==1)Uk(i)=satu(Uk(i),S(i));endif(Z(i)==2)Uk(i)=dead(Uk(i),S(i));endif(Z(i)==3)[Uk(i),Ubb(i)]=backlash(Ubb(i),Uk(i),Ub(i),S(i));endendendUdot=(Uk-Ub)/h;Uf=2*Uk-Ub;X=FI'.*X+FIM'.*Uk+FIJ'.*Udot;