第八章RLC电路与常微分方程的解法8.1RC电路与常微分方程的欧拉解法RC电路:K21RC先把开关K接通“1”端,电容C充满电后再把开关K接通“2”端,则这时电容C放电过程满足方程:即电容C上的电量是时间t的函数,满足以上微分方程.0dQQRdtC如果设:=RC,t=0时刻电容所带电量为Q0则有:考虑数值微分问题:已知:求f(x)在xn点的导数.可以:或:0(0),0dQQdtQQt11()()()nnnnnfxfxfxxx11()()()nnnnnfxfxfxxx123,,,...xxx123(),(),(),...fxfxfx微分方程化为一般形式:把时间t等间隔离散化:其中:做如下近似:由方程得:000(,)(),dQfQtdtQtQtt012,,,......ttt10,ttt202,......ttt11()()()nnnnttnnQtQtdQQtdttt11()()((),)nnnnnnQtQtfQtttt即:记:则得到解微分方程的欧拉法递推公式:对于RC电路:令1()()((),)nnnnQtQtfQttt11()nnQtQ()nnQtQ10000(,)(),.nnnnQQfQttQQttt100(0),0.nnnQQQtQQt01,10,1QRCt得到:方程的解析解:100101,0.nnnQQQQt0()tQtQe微分方程化为一般形式:把时间t等间隔离散化:其中:欧拉(Euler)差分公式:由方程得:000(,)(),dQfQtdtQtQtt012,,,......ttt10,ttt202,......ttt1()()nnnnttttQtQtdQQdttt1()()((),)nnnnQtQtfQttt即:记:则得到解微分方程的欧拉法递推公式:对于RC电路:例如:1()()((),)nnnnQtQtfQttt11()nnQtQ()nnQtQ10000(,)(),.nnnnQQfQttQQttt100(0),0.nnnQQQtQQt01,10,1QRCt得到:100101,0.nnnQQQQtrc(1,6,1,10);01234560.50.550.60.650.70.750.80.850.90.951欧拉法也可解释为Q(t)在tn处的泰勒展开:取线性部分:欧拉方法的截断误差:21()()()()2nnnnQttQtQttQtt1()()()nnnQtQtQtt2()()OtOt局部:整体:例:写出解如下一阶常微分方程的欧拉公式:得:2(0)1,0,0.1,dyxydxyyxx10020.11,0.nnnnnxyyyyyx8.2RLC电路和改进的欧拉近似法RLC电路图:LRCVaK根据基尔霍夫定律:由于:RLCaVVVV,,RLCdIQVIRVLVdtC得:由于:所有:欧拉法:把二阶微分方程化成一阶微分方程组:22adQdQQLRVdtdtC00000,,(),()1()adQIdtttItIQtQdIQVIRdtLCdQIdtadIQLIRVdtC其中t是自变量,Q和I随着t的改变而改变.1100000()(),(),nnnnnnanQQtIQtIIVIRLCIItQQtttfunction[Q,I,tt]=rlc(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;forn=1:length(tt)-1Q(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-10123456rlc(1,0,[1,5,1,5],50,0.1);05101520253035404550-10123456782.向后的欧拉方法方法分为两步:预估:(一步)校正:00000,,(),()1()adQIdtttItIQtQdIQVIRdtLC0101()nnnnnnanQQtIQtIIVIRLC01100111()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;forn=1:length(tt)-1Q1=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);holdonrlc(1,0,[1,1,1,5],15,0.1);051015-101234563.改进的欧拉法方法分两步:预估:(一步)校正:00000,,(),()1()adQIdtttItIQtQdIQVIRdtLC0101()nnnnnnanQQtIQtIIVIRLC01100111()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;forn=1:length(tt)-1Q1=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:');051015-10123456RC电路:向后的欧拉法:预估:校正:改进的欧拉法:预估:校正:0(0),0dQQdtQQt01nnnQQQt011nnnQQQt01nnnQQQt011()2nnnnQQQQtfunction[Q1,Q2,Q3,tt]=rc3(Q0,T,dt,tao)%RC电路欧拉解法Q1(1)=Q0;Q2(1)=Q0;Q3(1)=Q0;tt=0:dt:T;forn=1:length(tt)-1Q1(n+1)=Q1(n)-dt*Q1(n)/tao;endforn=1:length(tt)-1Q=Q2(n)-dt*Q2(n)/tao;Q2(n+1)=Q2(n)-dt*Q/tao;endforn=1:length(tt)-1Q=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.80.850.90.951一般微分方程:向后的欧拉法:改进的欧拉法:000(,)(),dQfQtdtQtQtt010111000(,)(,)(),nnnnnnnnQQtfQtQQtfQtQQttt010111000(,)[(,)(,)]2(),nnnnnnnnnnQQtfQtfQtfQtQQtQQttt8.3龙格-库塔(R-K)方法对于微分方程:根据微分中值定理:即:000(,)(),dQfQtdtQtQtt111()()(),[,]nnnnnnQtQtQtttt111()()()(),[,]nnnnnnQtQtttQtt11((),),[,]nnnnQQtfQttQ(t)tntn+1用tn处Q(t)的导数代替处导数f(,Q()),则为欧拉法:用tn+1处Q(t)的导数的估计值代替处导数f(,Q()),则为向后的欧拉法:1(,)nnnnQQtfQt即:用tn和tn+1处Q(t)的导数的估计值的平均代替处导数f(,Q()),则为改进的欧拉法:若取多点处斜率(即导数)的加权平均会使误差更小,称为龙格-库塔法010111(,)(,)nnnnnnnnQQtfQtQQtfQt010111(,)[(,)(,)]2nnnnnnnnnnQQtfQtfQtfQtQQt最常用的四阶龙格-库塔法:1123412132413(22)6(,)(,)22(,)22(,)nnnnnnnnnntQQkkkkkftQttkftQkttkftQkkftQtk000(,)(),dQfQtdtQtQtt例1:求解方程:梯形法:四阶龙格-库塔法:(0)1dyydxy010110[]21,1,0nnnnnnnyyxyyyyyxyxx有:1123412132430(22)6221,1,0nnnnnnxyykkkkkyxkykxkykkyxkyxxfunction[y1,y2,xx]=rk1(y0,X,dx)%矩形法和四阶龙格--库塔法y1(1)=y0;y2(1)=y0;xx=0:dx:X;forn=1:length(xx)-1y=y1(n)+dx*y1(n);y1(n+1)=y1(n)+dx*(y1(n)+y)/2;endforn=1:length(xx)-1k1=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(2'yyxyy2.0h