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

常微分方程的数值解法 绪论 在工程和科学计算中，所建立的各 种常微分方程的初值或边值问题，除很 少几类的特殊方程能给出解析解，绝大 多数的方程是很难甚至不可能给出解析 解的，其主要原因在于积分工具的局限 性。因此，人们转向用数值方法去解常 微分方程，并获得相当大的成功，讨论 和研究常微分方程的数值解法是有重要 意义的。 常微分方程描写的物理现象 n镭的衰变规律 n单摆的运动 nRLC振荡电路 n物理场计算 常微分方程 常微分方程数值解基本思想 8.1 Euler方法 8.1 Euler方法 8.1 Euler方法-梯形公式 8.2 改进Euler方法 8.2 改进Euler方法 12K 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),.);%理论曲线 legend(欧拉法的曲线,改进欧拉法的曲线, 理论曲线) 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）方法 n思想：取多点处斜率的加权平均为平均斜 率，从而减小误差。 n四阶公式： 公式推导 见P37-44 n例：用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); n例：用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; 练习： n分别用Euler法、改进Euler法和四阶R-K法求 解阻尼振动方程： 已知质量m=10,倔强系数k=10,阻尼系数c=2, 初始速度v=0，初始位置x=10. 误差概述 误差概述 误差概述 误差概述 8.1.3 数值稳定性分析 数值稳定性分析 n定义8.1.3 若某数值算法的绝对稳定性区 域包含h平面上的左半平面Re(h)0， 则称该方法是A稳定的。 n隐式Euler法是A稳定的。 8.2 Runge-Kutta方法 Runge-Kutta方法 Runge