,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,*,单击此处编辑母版文本样式,第二级,“,”,常微分方程的数值解法,绪论,在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。,常微分方程描写的物理现象,镭的衰变规律,单摆的运动,伯努利方程,RLC,振荡电路,物理场计算,常微分方程,常微分方程数值解基本思想,8.1 Euler,方法,8.1 Euler,方法,8.1 Euler,方法,-,梯形公式,8.2,改进,Euler,方法,8.2,改进,Euler,方法,1,2,K,C=0.01,法,R=10,欧,E=10,伏,function,main,global E R C;,E=10;R=10;C=0.01;,Q0=E/C;,h=0.01;,t=0:h:1;,QE(1)=Q0;,QEG(1)=Q0;,for i=2:length(t),QE(i,)=QE(i-1)+h*f(QE(i-1);%,欧拉法,k1=h*f(QEG(i-1);%,改进欧拉法,k2=h*f(QEG(i-1)+k1);,QEG(i,)=QEG(i-1)+1/2*(k1+k2);,end,plot(t,QE,r,);%,改进欧拉法的曲线,hold on,plot(t,QEG,b,);%,改进欧拉法的曲线,plot(t,Q0*exp(-1*t./(R*C),.);%,理论曲线,function,y=,f(x,),global E R C;,y=-1*x./(R*C);,K,C=10mF,R=100,欧,E=10V,L=10H,function,main,%,欧拉法：二阶常微分方程,global E R C L;,E=10;R=100;C=0.01;L=10;,Q0=0;I0=0;h=0.01;t=0:h:10;,Q(1)=Q0;,I(1)=I0;,for i=2:length(t),Q(i,)=Q(i-1)+h*I(i-1);%Q,I(i,)=I(i-1)+h*f(Q(i-1),I(i-1);%I,end,plot(t,Q,r,t,I,m,);%,改进欧拉法的曲线,function,y=,f(Q,I,),global E R C L;,y=(E-Q/C-I*R)/L;,（用改进的,Euler,法解）：,function,main,%,改进欧拉法：二阶常微分方程,global E R C L;,E=10;R=100;C=0.01;L=10;,Q0=0;I0=0;h=0.01;t=0:h:10;Q(1)=Q0;I(1)=I0;,for i=2:length(t),%1,、预报,Q(i,)=Q(i-1)+h*I(i-1);,I(i,)=I(i-1)+h*f(Q(i-1),I(i-1);,%2,、计算,Q,k1=h*I(i-1);k2=h*,I(i,);,Q(i,)=Q(i-1)+1/2*(k1+k2);,%3,、计算,I,k1=h*f(Q(i-1),I(i-1);k2=h*,f(Q(i),I(i,);,I(i,)=I(i-1)+1/2*(k1+k2);,end,plot(t,Q,r,t,I,m,);%,改进欧拉法的曲线,function,y=,f(Q,I,),global E R C L;,y=(E-Q/C-I*R)/L;,8.3,龙格,-,库塔（,R-K,）方法,思想：取多点处斜率的加权平均为平均斜率，从而减小误差。,四阶公式：,公式推导,见,P37-44,例：用,R-K,方法解例题,8.1.1,function,main,%,龙格库塔：一阶常微分方程,global E R C;,E=10;R=10;C=0.01;h=0.01;t=0:h:1;Q(1)=E/C;,for i=2:length(t),k1=h*f(Q(i-1);,k2=h*f(Q(i-1)+k1/2);,k3=h*f(Q(i-1)+k2/2);,k4=h*f(Q(i-1)+k3);,Q(i,)=Q(i-1)+1/6*(k1+2*k2+2*k3+k4);,end,plot(t,Q,b);hold,on%,龙格库塔方法的曲线,plot(t,Q(1)*exp(-1*t./(R*C),.);%,理论曲线,function,y=,f(x,),global E R C;,y=-1*x./(R*C);,例：用,R-K,方法解例题,8.1.2,和,function main,%,龙格库塔：二阶常微分方程,global E R C L;,E=10;R=100;C=0.01;L=10;,Q0=0;I0=0;h=0.01;t=0:h:10;Q(1)=Q0;I(1)=I0;,for i=2:length(t),k1=h*I(i-1);m1=h*f(Q(i-1),I(i-1);,k2=h*(I(i-1)+m1/2);m2=h*f(Q(i-1)+k1/2,I(i-1)+m1/2);,k3=h*(I(i-1)+m2/2);m3=h*f(Q(i-1)+k2/2,I(i-1)+m2/2);,k4=h*(I(i-1)+m3);m4=h*f(Q(i-1)+k3,I(i-1)+m3);,Q(i,)=Q(i-1)+1/6*(k1+2*k2+2*k3+k4);,I(i,)=I(i-1)+1/6*(m1+2*m2+2*m3+m4);,end,plot(t,Q,r,t,I,m,);%,龙格库塔的曲线,function y=,f(Q,I,),global E R C L;,y=(E-Q/C-I*R)/L;,练习：,分别用,Euler,法、改进,Euler,法和四阶,R-K,法求解阻尼振动方程：,已知质量,m=10,倔强系数,k=10,阻尼系数,c=2,初始速度,v=0,，初始位置,x=10.,误差概述,误差概述,误差概述,误差概述,8.1.3,数值稳定性分析,数值稳定性分析,定义,8.1.3,若某数值算法的绝对稳定性区域包含,h,平面上的左半平面,Re(h,),0,，则称该方法是,A,稳定的。,隐式,Euler,法是,A,稳定的。,8.2,Runge-Kutta,方法,Runge-Kutta,方法,Runge-Kutta,方法,Runge-Kutta,方法,6.2.2,四阶,Runge-Kutta,方法,四阶,Runge-Kutta,方法,6.2.3 R-K,法的稳定性,R-K,法的稳定性,R-K,法的稳定性,6.2.5,隐式,R-K,法,隐式,R-K,法,隐式,R-K,法,隐式,R-K,法,隐式,R-K,法,