第八章RLC电路与常微分方程的01574

1、第八章 rlc电路与常微分方程的解法 8.1 rc电路与常微分方程的欧拉解法rc电路: k 2 1 r c 先把开关k接通“1” 端,电容c充满电后再把开关k接通“2”端,则这时电容c放电过程满足方程: 即电容c上的电量是时间t的函数,满足以上微分方程.0dqqrdtc如果设: =rc, t=0时刻电容所带电量为q0 则有:考虑数值微分问题:已知: 求f(x) 在xn 点的导数.可以:或:0(0),0dqqdtqq t 11()()()nnnnnf xf xfxxx11()()()nnnnnf xf xfxxx123,.x x x123(),(),(),.f xf xf x微分方程化为一般形式

2、:把时间t等间隔离散化:其中:做如下近似: 由方程得:000( , )( ),dqf q tdtq tq tt012, , ,.t t t10,ttt202,.ttt 11( )()( )nnnnt tnnq tq tdqq tdttt11( )()( ( ),)nnnnnnq tq tf q tttt即即:记记:则得到解微分方程的欧拉法递推公式则得到解微分方程的欧拉法递推公式:对于对于rc电路电路:例如：例如：1()( )( ( ),)nnnnq tq tf q ttt11()nnq tq( )nnq tq10000(,)( ),.nnnnqqf q ttqq ttt100(0),0.nnn

3、qqqtqqt01, 10, 1qrct 得到:方程的解析解:100101,0.nnnqqqqt0( )tq tq e rc(1,6,1,10);01234560.50.550.60.650.70.750.80.850.90.951欧拉法也可解释为欧拉法也可解释为q(t)在在tn处的处的泰勒展开泰勒展开:取线性部分取线性部分:欧拉方法的截断误差欧拉方法的截断误差:21()( )( )( )2nnnnq ttq tq ttq tt 1()( )( )nnnq tq tq tt2()()otot局部：整体：例: 写出解如下一阶常微分方程的欧拉公式:得:2(0)1,0,0.1,dyxydxyyxx

4、10020.11,0.nnnnnxyyyyyx 8.2 rlc电路和改进的欧拉近似法rlc 电路图: l r c va k 根据基尔霍夫定律:由于:rlcavvvv,rlcdiqvir vlvdtc得:由于:所有:w欧拉法:把二阶微分方程化成一阶微分方程组:22ad qdqqlrvdtdtc00000, ( ),( )1()adqidttti ti q tqdiqvirdtlcdqidtadiqlirvdtc其中t是自变量,q和i随着t的改变而改变.1100000()( ),( ),nnnnnnanqqtiqtiivi rlcii tqq tttfunction q,i,tt=rlc(q0,i

5、0,con,t,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 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);051015-10123456 rlc(1,0,1,5,1,5,50,0.1);05101520253035404550-10123456782. 向后的欧拉方法方法分为两步:预估:

6、(一步)校正:00000, ( ),( )1()adqidttti ti q tqdiqvirdtlc0101()nnnnnnanqqtiqtiivi rlc01100111()nnnnnnanqqtiqtiivirlc或者(k+1步)校正:1111111111111(),kknnnkkknnnankknnnnqqtiqtiivirlcqqiifunction q,i,tt=rlc1(q0,i0,con,t,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(t

7、t)-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);051015-101234563. 改进的欧拉法方法分两步:预估:(一步)校正:00000, ( ),( )1()adqidttti ti q tqdiqvirdtlc0101()nnnnnnanqqtiqtiivi rlc

8、01100111()2()()()22nnnnnnnnnnaiiqqtiiqqtiivrlc或(k+1步)校正:1111111111111()2()()()22,kknnnnkkknnnnnnakknnnniiqqtiiqqtiivrlcqqiifunction q,i,tt=rlc2(q0,i0,con,t,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)

9、/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:);051015-10123456rc电路:向后的欧拉法:预估:校正:改进的欧拉法:预估:校正:0(0),0dqqdtqq t 01nnnqqqt011nnnqqqt01nnnqqqt011()2nnnnqqqqtfunction q1,q2,q3,tt=rc3(q0,t,dt,tao)% rc电路欧拉解法电路欧拉解法q1(1)=q0;q2(1)=q0;q3(1)=q0;tt=0:dt:t;f

10、or 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,r-,tt,q2,r-,tt,q3,r:); rc3(1,6,1,10)01234560.50.550.60.650.70.750

11、.80.850.90.951一般微分方程:向后的欧拉法:改进的欧拉法:000( , )( ),dqf q tdtq tq tt010111000(,)(,)( ),nnnnnnnnqqtf q tqqtf qtqq ttt010111000(,) (,)(,)2( ),nnnnnnnnnnqqtf q tf q tf qtqqtqq ttt 8.3 龙格-库塔(r-k)方法对于微分方程:根据微分中值定理:即:000( , )( ),dqf q tdtq tq tt111()( )( ), ,nnnnnnq tq tqt ttt111()( )()( ), ,nnnnnnq tq ttt qt

12、t11( ( ), ), ,nnnnqqtf qt t q(t) tn tn+1 用tn处q(t)的导数代替处导数 f(,q(),则为欧拉法: 用tn+1处q(t)的导数的估计值代替处导数 f(,q(),则为向后的欧拉法:1(,)nnnnqqtf q t即: 用tn和tn+1处q(t)的导数的估计值的平均代替处导数 f(,q(),则为改进的欧拉法: 若取多点处斜率(即导数)的加权平均会使误差更小,称为龙格-库塔法010111(,)(,)nnnnnnnnqqtf q tqqtf qt010111(,) (,)(,)2nnnnnnnnnnqqtf q tf q tf qtqqt最常用的四阶龙格-库

13、塔法:1123412132413(22)6( ,)(,)22(,)22(,)nnnnnnnnnntqqkkkkkf t qttkf tqkttkf tqkkf tqtk例1: 求解方程:梯形法:四阶龙格-库塔法:(0)1dyydxy01011021,1,0nnnnnnnyyx yyyyyxyxx 有:1123412132430(22)6221,1,0nnnnnnxyykkkkkyxkykxkykkyxkyxx function y1,y2,xx=rk1(y0,x,dx)% 矩形法和四阶龙格矩形法和四阶龙格-库塔法库塔法y1(1)=y0;y2(1)=y0;xx=0:dx:x;for n=1:le

14、ngth(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:); rk1(1,5,1);00.511.522.533.544.55050100150例2：四阶的龙格-库塔公式1)0(2yyxyy2 . 0h3342

15、231121432112 . 0)2 . 0(22 . 01 . 0) 1 . 0(21 . 01 . 0) 1 . 0(21 . 02)22(31 . 0kyxkykkyxkykkyxkykyxykkkkkyynnnnnnnnnnnnnn例3：四阶的龙格-库塔公式0)0(20314sin622iitdtdi)001. 0(20)001. 0(314sin622)0005. 0(20)0005. 0(314sin622)0005. 0(20)0005. 0(314sin62220314sin622)22(6001. 0342312143211kitkkitkkitkitkkkkkyinnnnn

16、nnnnn微分方程组:( , , )( , , )dxf x y tdtdyg x y tdt龙格-库塔公式:112341123411211211322322(22)6(22)6(,);(,)(,);(,)222222(,);(,22222nnnnnnnnnnnnnnnnnnnnnntxxkkkktyyllllkf xy tlg xy tttttttkf xk yl tlg xk yl ttttttkf xkyl tlg xkyl t43314331)2(,);(,)nnnnnntkf xtkytl tlg xtkytl t例4: 求解阻尼振动方程首先把它转化为一阶微分方程组:22d xdxm

17、ckxdtdt 1()dxvdtdvcvkxdtm四阶龙格-库塔公式:112341123411212113232243433(22)6(22)61;()1;()()2221;()()2221;()()nnnnnnnnnnnnnnnntxxkkkktvvllllkv lcvkxmtttkvl lc vlk xkmtttkvl lc vlk xkmkvtk lc vtlk xtkmm0010,10,2,10,0,0kcxvtfunction 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(

18、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(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