CP080-计算物理常微分方程解.ppt

常微分方程的数值解法,绪论,在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。,常微分方程描写的物理现象,镭的衰变规律单摆的运动RLC振荡电路物理场计算,常微分方程,常微分方程数值解基本思想,8.1Euler方法,8.1Euler方法,8.1Euler方法-梯形公式,8.2改进Euler方法,8.2改进Euler方法,functionmainglobalERC;E=10;R=10;C=0.01;Q0=E/C;h=0.01;t=0:h:1;QE(1)=Q0;QEG(1)=Q0;fori=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);endplot(t,QE,r);%欧拉法的曲线holdonplot(t,QEG,b);%改进欧拉法的曲线plot(t,Q0*exp(-1*t./(R*C),.);%理论曲线legend(欧拉法的曲线,改进欧拉法的曲线,理论曲线)functiony=f(x)globalERC;y=-1*x./(R*C);,functionmain%欧拉法：二阶常微分方程globalERCL;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;fori=2:length(t)Q(i)=Q(i-1)+h*I(i-1);%QI(i)=I(i-1)+h*f(Q(i-1),I(i-1);%Iendplot(t,Q,r,t,I,m);%欧拉法的曲线functiony=f(Q,I)globalERCL;y=(E-Q/C-I*R)/L;,（用改进的Euler法解）：,functionmain%改进欧拉法：二阶常微分方程globalERCL;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;fori=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、计算Qk1=h*I(i-1);k2=h*I(i);Q(i)=Q(i-1)+1/2*(k1+k2);%3、计算Ik1=h*f(Q(i-1),I(i-1);k2=h*f(Q(i),I(i);I(i)=I(i-1)+1/2*(k1+k2);endplot(t,Q,r,t,I,m);%改进欧拉法的曲线functiony=f(Q,I)globalERCL;y=(E-Q/C-I*R)/L;,8.3龙格-库塔（R-K）方法,思想：取多点处斜率的加权平均为平均斜率，从而减小误差。四阶公式：,公式推导见P37-44,例：用R-K方法解例题8.1.1,functionmain%龙格库塔：一阶常微分方程globalERC;E=10;R=10;C=0.01;h=0.01;t=0:h:1;Q(1)=E/C;fori=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);endplot(t,Q,b);holdon%龙格库塔方法的曲线plot(t,Q(1)*exp(-1*t./(R*C),.);%理论曲线functiony=f(x)globalERC;y=-1*x./(R*C);,例：用R-K方法解例题8.1.2,和,functionmain%龙格库塔：二阶常微分方程globalERCL;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;fori=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);endplot(t,Q,r,t,I,m);%龙格库塔的曲线functiony=f(Q,I)globalERCL;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.2Runge-Kutta方法,Runge-Kutta方法,Runge-Kutta方法,R