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

1、常微分方程的数值解法绪论 在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因而，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。常微分方程描写的物理现象n镭的衰变规律n单摆的运动nRLC振荡电路n物理场计算sinmgdtdmdtdlECIRILIQdtddtdQ常微分方程常微分方程数值解基本思想8.1 Euler方法。算出出发，逐次公式，它可以从这就是显式的）（的近似值，则有表示以于是该式可离散为：代替散化，并用方法首先将

2、微分算子离为实现这一目标，.,Euler) 1 (,.2 , 1 , 0),()()(,()()(|,Euler321010yyyynyxhfyyxyyxyxfhxyhxydxdyxxnnnnnnnnnnxnnnnnnxxhxxxxnxyxxxyxxxyxyxy000)(200000.)(!)()(! 2)( )( )()(8.1 Euler方法。才能得到步要解函数方程的不同在于，它每算一方法，它与显式公式或向后这就是隐式的）（的近似值，则有表示以于是该式可离散为：代替如果用11111101)2(EulerEuler)2(,.2 , 1 , 0),()()(,()()(,nnnnnnnnnnn

3、nynyxhfyyxyyxyxfhxyhxyxxhxxxxxxnxyxxxyxxxyxyxynnnn000)(200000.)(!)()(! 2)( )( )()()()()()(1)()(11,|3,.,),(),(),(,.,),(),(knnknknknnnnnknnnnknnnnnnnnyyyykyxfyxfhyyyxhfyyynyxfyxfhyy取时当）（）（时常用以下迭代式：计算称为梯形公式。）（平均值的结果，则得如果取以上两式的算术8.1 Euler方法-梯形公式。并取即校正算法构成一类预估迭代一次一般只由于迭代工作量较大计算对于)1(11)0(11)1( 1)0(11),(),

4、(2),(,)2(nnnnnnnnnnnnnyyyxfyxfhyyyxhfyyy8.2 改进Euler方法),(,()(,Euler,.),(),(),()(nnnnnnnnnnnnnnyxhfyxfyxfhyynkyhxhfkyxhfkkkyy亦可写成方法该式称为改进上式还常写成8.2 改进Euler方法12KC=0.01法R=10欧E=10伏,.)2 , 1 , 0()()()()()(21Euler,.)2 , 1 , 0()( 1 , 0,01. 0;01. 0/10/)(11212111000nRCkQhkQhfkRCQhQhfkkkQQnRCQhQQEulerthCEQQtQRCQ

5、dtdQnnnnnnnnn法：改进的法：上结果，此时计算取初始状态：下常微分方程解：分析得电量满足如function mainglobal 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);endplot(t,QE,r);%欧拉法的曲线欧拉法的曲线hold onp

6、lot(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);KC=10mFR=100欧E=10VL=10H ,.)2 , 1 , 0( ,/ )/(,.)2 , 1 , 0( ,Euler10, 0,01. 0; 0I ; 01100nLRICQEhIInhIQQthQnnnnnnn法：上结果，此时计算取初始状态：function main%欧

7、拉法：二阶常微分方程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);%Iendplot(t,Q,r,t,I,m);%欧拉法的曲线function y=f(Q,I)global E R C L;y=(E-Q/C-I*R)/L;（用改进的Euler法解）：,.),(/)(,.),()()(EulernLRICQEhdtdIhkLRICQEhdtd

8、IhkkkIInIhdtdQhkhIdtdQhkkkQQnnnnnnnnnnnnnn）（）（和法：、改进的,.),( ,/)/(,.),( ,EulernLRICQEhIInhIQQnnnnnnn法预报：、先用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

9、-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);endplot(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四阶公式：n ),()21,21()21,21(),()22(61

10、342343kyhxhfkkyhxhfkkyhxhfkyxhfkkkkkyynnnnnnnnnn公式推导见P37-44n例：用R-K方法解例题8.1.1)RC(),()RC/()21,21()RC/()21,21()RC(),()22(61342343kQhkQhthfkkQhkQhthfkkQhkQhthfkQhQthfkkkkkQQnnnnnnnnnnnnnnfunction 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(

11、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);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)21(),()21()21,21()21()21,21(),()22(61342343lIhlIhthfklIh

12、lIhthfklIhlIhthfkhIIthfkkkkkQQnnnnnnnnnnnnnn),( )21,21,21( )21,21,21( / )/(),( )22(61342343kQlIhthflkQlIhthflkQlIhthflLCQRIEhQIthflllllIInnnnnnnnnnnnnnnn和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

13、); 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);%龙格库塔的曲线龙格库塔的曲线function y=f(Q,I)

14、global E R C L;y=(E-Q/C-I*R)/L;练习：n分别用Euler法、改进Euler法和四阶R-K法求解阻尼振动方程：n n 已知质量m=10,倔强系数k=10,阻尼系数c=2,初始速度v=0，初始位置x=10.kxdtdxcdtxdm误差概述有关，称为增量函数。与函数为而隐式单步法一般形式显式单步法一般形式为),(),(),(1111yxfhyxyxhyyhyxhyynnnnnnnnnn误差概述上的局部截断误差。称为点则由单步法计算出值计算没有误差，即的步在点果第上的整体截断误差。如称为在点则的值或隐式逐步计算，得出发，由单步法显式从初值定义nnnnnnnnnnnnnnx

15、yxyTyxyyxnxyxyeyxyxy,)(,),()(,)(8.1.1误差概述差截断如果给定的算法的局部称为为局部截断误)(,(则非零项)()(阶精度。如果则称该算法具有误差为2 1. . 8定义pnnppnnnpnhxyxghOhx,yxgTp)O(hT误差概述)()(2)O()(y2Euler式对于)O(h)(2Euler对于显式则精度相对的越高。越大，误差阶，一个算法，局部截断一般431321321hOxyhThxhTxyhTnnnnnn ：对于梯形公式：法隐：法说来8.1.3 数值稳定性分析稳定。值不大，则，则称该算法的舍入误的舍入误差。如果算则称该算法是数值稳定|，得值方法又算，

16、如果结果为计算，则时有一舍入误差假设1111数，也成绝对稳定且通过某种数实际计算nnnnnnnnnnyyyyy数值稳定性分析n定义8.1.3 若某数值算法的绝对稳定性区域包含h平面上的左半平面Re(h)0，则称该方法是A稳定的。n隐式Euler法是A稳定的。8.2 Runge-Kutta方法。这里仍假定，使局部截断误差，适当选择参数可设为方法启发，更一般算式受改进的)(),()(,.)2 , 1 , 0(),(),(Euler31112112122111nnnnnnnnnnnxyyhOyxyTnkyhxhfkyxhfkkkyyRunge-Kutta方法)(),(),(),()()(),(),(

17、),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfknnynnnnxnnnynnxnn展开式由二元函数)(),(),(),()()()(3222211hOhyxfyxfyxfxyhxyynnynnnnxnnn于是Runge-Kutta方法一。自由参数，即解答不唯程，因此有一个由于四个参数，三个方展式相比较得与21211:Taylor2221方法。这是改进的此时算式为可得取Euler),(),()(21, 1,21,21) 1 (12121121kyhxhfkyxhfkkkyynnnnnnRunge-Kutta方法.K-R)21,21(),(,21, 1

18、, 0)2(1212121方法这是二阶此时算式为可得取kyhxhfkyxhfkkyynnnnnn方法。这也是二阶（又有算式可得取K-R)32,32(),()341,32,43,41) 3(12121121kyhxhfkyxhfkkkyynnnnnn6.2.2 四阶Runge-Kutta方法4342414323132123432421414323213133121221443322111,.)2 , 1 , 0(),(),(),(),(KR，约定分方程初值问题可设为方法类似，解一阶常微与二阶nkkkyhxhfkkkyhxhfkkyhxhfkyxhfkkkkkyynnnnnnnnnn四阶Runge

19、-Kutta方法稳定。但不是该公式结构整齐算式求得经典四阶类似以上推导过程，可A,),()21,21()21,21(),()22(61KR332312143211kyhxhfkkyhxhfkkyhxhfkyxhfkkkkkyynnnnnnnnnn6.2.3 R-K法的稳定性nnnnyhhhhkyhhhkyhhkyhkyxy)(41)(21)()(41)(21()(21(dd,K-R4324323221有应用于实验方程法为例以经典R-K法的稳定性nnnnnhhhhyhhhhkkkkyy! 4)(! 3)(! 2)()(1 :! 4)(! 3)(! 2)()(1 )22(61432143243211相应的误差为于是R-K法的稳定性., 078. 2:,1|! 4)(!