新編-第八章RLC電路與常微分方程的-精品課件

1、第八章 RLC電路與常微分方程的解法 8.1 RC電路與常微分方程的歐拉解法RC電路: K 2 1 R C 先把開關(guān)K接通“1” 端,電容C充滿電后再把開關(guān)K接通“2”端,則這時電容C放電過程滿足方程: 即電容C上的電量是時間t的函數(shù),滿足以上微分方程.如果設(shè): =RC, t=0時刻電容所帶電量為Q0 則有:考慮數(shù)值微分問題:已知: 求f(x) 在xn 點的導(dǎo)數(shù).可以:或:微分方程化為一般形式:把時間t等間隔離散化:其中:做如下近似: 由方程得:即:記:則得到解微分方程的歐拉法遞推公式:對于RC電路:令得到:方程的解析解:微分方程化為一般形式:把時間 t 等間隔離散化:其中:歐拉(Euler)

2、 差分公式: 由方程得:即:記:則得到解微分方程的歐拉法遞推公式:對于RC電路:例如：得到: rc(1,6,1,10);歐拉法也可解釋為Q(t)在tn處的泰勒展開:取線性部分:歐拉方法的截斷誤差:例: 寫出解如下一階常微分方程的歐拉公式:得: 8.2 RLC電路和改進的歐拉近似法RLC 電路圖: L R C Va K 根據(jù)基爾霍夫定律:由于:得:由于:所有:歐拉法:把二階微分方程化成一階微分方程組:其中t是自變量,Q和I隨著t的改變而改變.function Q,I,tt=rlc(Q0,I0,con,T,dt)% RLC電路歐拉解法Q(1)=Q0;I(1)=I0;R=con(1);L=con(2

3、);C=con(3);V=con(4);tt=0:dt:T;for n=1:length(tt)-1 Q(n+1)=Q(n)+dt*I(n); I(n+1)=I(n)+dt*(V-R*I(n)-Q(n)/C)/L;endplot(tt,Q,r,tt,I,b); rlc(1,0,1,1,1,5,15,0.1); rlc(1,0,1,5,1,5,50,0.1);2. 向后的歐拉方法方法分為兩步:預(yù)估:(一步)校正:或者(k+1步)校正:function Q,I,tt=rlc1(Q0,I0,con,T,dt)% RLC電路向后歐拉解法Q(1)=Q0;I(1)=I0;R=con(1);L=con(2)

4、;C=con(3);V=con(4);tt=0:dt:T;for n=1:length(tt)-1 Q1=Q(n)+dt*I(n); I1=I(n)+dt*(V-R*I(n)-Q(n)/C)/L; Q(n+1)=Q(n)+dt*I1; I(n+1)=I(n)+dt*(V-R*I1-Q1/C)/L;endplot(tt,Q,r-,tt,I,b-); rlc1(1,0,1,1,1,5,15,0.1); hold on rlc(1,0,1,1,1,5,15,0.1);3. 改進的歐拉法方法分兩步:預(yù)估:(一步)校正:或(k+1步)校正:function Q,I,tt=rlc2(Q0,I0,con,T

5、,dt)% RLC電路改進歐拉解法Q(1)=Q0;I(1)=I0;R=con(1);L=con(2);C=con(3);V=con(4);tt=0:dt:T;for n=1:length(tt)-1 Q1=Q(n)+dt*I(n); I1=I(n)+dt*(V-R*I(n)-Q(n)/C)/L; Q(n+1)=Q(n)+dt*(I1+I(n)/2; I(n+1)=I(n)+dt*(V-R*(I1+I(n)/2- (Q1+Q(n)/2/C)/L;endplot(tt,Q,r:,tt,I,b:);RC電路:向后的歐拉法:預(yù)估:校正:改進的歐拉法:預(yù)估:校正:function Q1,Q2,Q3,tt

6、=rc3(Q0,T,dt,tao)% RC電路歐拉解法Q1(1)=Q0;Q2(1)=Q0;Q3(1)=Q0;tt=0:dt:T;for n=1:length(tt)-1 Q1(n+1)=Q1(n)-dt*Q1(n)/tao;endfor n=1:length(tt)-1 Q=Q2(n)-dt*Q2(n)/tao; Q2(n+1)=Q2(n)-dt*Q/tao;endfor n=1:length(tt)-1 Q=Q3(n)-dt*Q3(n)/tao; Q3(n+1)=Q3(n)-dt*(Q+Q3(n)/2/tao;endQa=Q0*exp(-tt/tao);plot(tt,Qa,b,tt,Q1,

7、r-,tt,Q2,r-,tt,Q3,r:); rc3(1,6,1,10)一般微分方程:向后的歐拉法:改進的歐拉法: 8.3 龍格-庫塔(R-K)方法對于微分方程:根據(jù)微分中值定理:即: Q(t) tn tn+1 用tn處Q(t)的導(dǎo)數(shù)代替處導(dǎo)數(shù) f(,Q(),則為歐拉法: 用tn+1處Q(t)的導(dǎo)數(shù)的估計值代替處導(dǎo)數(shù) f(,Q(),則為向后的歐拉法:即: 用tn和tn+1處Q(t)的導(dǎo)數(shù)的估計值的平均代替處導(dǎo)數(shù) f(,Q(),則為改進的歐拉法: 若取多點處斜率(即導(dǎo)數(shù))的加權(quán)平均會使誤差更小,稱為龍格-庫塔法最常用的四階龍格-庫塔法:例1: 求解方程:梯形法:四階龍格-庫塔法:有:functi

8、on y1,y2,xx=rk1(y0,X,dx)% 矩形法和四階龍格-庫塔法y1(1)=y0;y2(1)=y0;xx=0:dx:X;for n=1:length(xx)-1 y=y1(n)+dx*y1(n); y1(n+1)=y1(n)+dx*(y1(n)+y)/2;endfor n=1:length(xx)-1 k1=y2(n); k2=y2(n)+dx*k1/2; k3=y2(n)+dx*k2/2; k4=y2(n)+dx*k3; y2(n+1)=y2(n)+dx*(k1+2*k2+2*k3+k4)/6;endy=exp(xx);plot(xx,y,b,xx,y1,r-,xx,y2,r:)

9、; rk1(1,5,1);例2：四階的龍格-庫塔公式例3：四階的龍格-庫塔公式微分方程組:龍格-庫塔公式:例4: 求解阻尼振動方程首先把它轉(zhuǎn)化為一階微分方程組:四階龍格-庫塔公式:function x,v,tt=rk2(m,k,c,x0,v0,T,dt)% 四階龍格-庫塔法x(1)=x0;v(1)=v0;tt=0:dt:T;for n=1:length(tt)-1 k1=v(n); l1=-(c*v(n)+k*x(n)/m; k2=v(n)+dt*l1/2; l2=-(c*(v(n)+dt*l1/2)+k*(x(n)+dt*k1/2)/m; k3=v(n)+dt*k2/2; l3=-(c*(v

10、(n)+dt*l2/2)+k*(x(n)+dt*k2/2)/m; k4=v(n)+dt*k3; l4=-(c*(v(n)+dt*l3)+k*(x(n)+dt*k3)/m; x(n+1)=x(n)+dt*(k1+2*k2+2*k3+k4)/6; v(n+1)=v(n)+dt*(l1+2*l2+2*l3+l4)/6;endplot(tt,x);h=line(Color,1 0 0, Marker,.,MarkerSize,20,EraseMode,xor);for i=1:length(tt) set(h,Xdata,tt(i),Ydata,x(i); pause(dt); end x,v,tt=rk2(10,10,2,10