建立系統(tǒng)微分方程一般步驟

1、建立系統(tǒng)微分方程一般步驟：(1)將系統(tǒng)劃分為多個(gè)環(huán)節(jié)，確定各環(huán)節(jié)的輸入及輸出信號(hào)，每個(gè)環(huán)節(jié)都可考慮寫一個(gè)方 程；(2)根據(jù)物理定律或通過實(shí)驗(yàn)等方法得出物理規(guī)律，列出各環(huán)節(jié)的原始方程式，并考慮適 當(dāng)簡(jiǎn)化、線性化；(3)將各環(huán)節(jié)方程式聯(lián)立，消去中間變量，最后得出只含有輸入變量、輸出變量以及參量 的系統(tǒng)方程式。建立LRC電路的微分方程式RLo1-Ii(t)1|+ Iui(t)Cu0用MATLAB語言編程實(shí)現(xiàn)仿真的主要步驟是調(diào)用MATLAB 中的 ODE (Ordinary Differential Equation,常微分方程)解函數(shù)。MATLAB 提供的常用ODE解函數(shù)如下：? ode45此算法被

2、推薦為首選算法；? ode23這是一個(gè)比ode45低階的算法；?ode113用于更高階或大的標(biāo)量計(jì)算；?ode23t用于解決難度適中的問題；?ode23s用于解決難度較大的問題；? ode15s 與ode23相同，但要求的精度更高；?ode23tb用于解決難度較大的問題。這些ODE解函數(shù)的調(diào)用格式基本相同。例如，ode45的基本調(diào)用格式為t , x=ode45('方程函數(shù)名',tspan , x0, tol)其中，方程函數(shù)名為描述系統(tǒng)狀態(tài)方程的M函數(shù)名稱，tspan一般為仿真時(shí)間范圍(例如，取tspan=t0,tf , t。為起始計(jì)算時(shí)間，tf為終止計(jì)算時(shí)間)；x0為系統(tǒng)狀態(tài)變

3、量初值，應(yīng)使 該向量元素個(gè)數(shù)等于系統(tǒng)狀態(tài)變量的個(gè)數(shù)；tol用來指定的精度，其默認(rèn)值為10-3 (即0.1%的相對(duì)誤差)，一般應(yīng)用中可以直接采用默認(rèn)值。函數(shù)返回兩個(gè)結(jié)果t向量和x陣。由于計(jì)算中采用了步長(zhǎng)自動(dòng)控制策略，因而t向量不一定是等間隔。但是，仿真結(jié)果可以用plot(t,x)指令繪制出來。例：電路如下圖所示，R=1.6C, L=2.1H, C=0.3F,初始狀態(tài)是電感電流為零，電容電壓為0.2V, t=0時(shí)接入1.5V的電壓，求0 ct <10s時(shí)i(t) , u0(t)的值，并畫出電流和電容電壓的關(guān)系曲線。Example ： Circuit is shown below, R=1.6

4、G, L=2.1H , C=0.3F , Initial condition: inductancecurrent is 0, capacitance voltage is 0.2V. Access 1.5V voltage when t=0, Evaluate i(t) , u0(t)when 0 : t ： 10s , and plot the curve between current and capacitance voltage.RL(t)i(t)+Ui(t)Cu0解：(1)根據(jù)基爾霍夫電壓定律、電流定律得到系統(tǒng)的原始微分方程為L(zhǎng) 則 Ri(t) uO(t) =u«) dt

5、i(t) = CdH dt(2)消去中間變量i(t),得電路微分方程：2LC d u0(t)rc du0(JI u0(t) =u)dtdt(3)劃高階微分方程為一階微分方程：設(shè) x1=u0(t), x2 = u0(t)則有11 = x21卡 x2 =(ui(t) -x1 -RCx2)LC(4)編程系統(tǒng)微分方程描述函數(shù)的程序：(rlcsys .m)function dx=rlcsys(t,x)dx=zeros(2,1);ui=1.5;R=1.6;L=2.1;C=0.3;dx(1)=x(2);dx(2)=1/(L*C)*(ui-x(1)-R*C*x(2);主程序( jierlc.m )clearc

6、lccloset0=0;tfinal=10;x0=0.2 0;t,x=ode45(rlcsys,t0 tfinal,x0);C=0.3;cfu=x(:,1); %C functionifu=C*x(:,2); %i functionfigure(1)subplot(211)plot(t,cfu)title( 'C/V' )xlabel( 'time/s' )gridsubplot(212)plot(t,ifu) title( 'L/A' )xlabel( 'time/s' )gridfigure(2)plot(cfu,ifu)gr

7、idtitle( 'C、L relation' )xlabel( 'C')ylabel( 'L')1、假設(shè)著名的Lorenz模型的狀態(tài)方程表示為xNt) = 8x1(t)/3+x2(t)x3(t)* x2(t) = 10x2(t) +10x3(t)；3(t)=x1(t)x2(t)+28x2(t)-x3(t)若令其初值為Xi(0=X2=0, x3=100,解方程。(采用ode45()函數(shù))Lorenz model state equation is shown below: xjt) = -8x(t)/3+x2(t)x3(t)；2(t)=-10X2(t) 10X3(t) *X3(t) = -X1(t)X2(t) 28X2(t) -X3(t)I，C,0Assume that initial values are x(t) = X2=0 , X3(t) =10 , please solve the equation.(Note:ode45()程序：(lorenz.m )function dx=lorenz(t,x)dx=zeros(3,1);dx(1)=-8*x(1)/3+x(2)*x(3);dx(2)=-10*x(2)+10*