《数值与算法》第八讲-常微分方程初值问题

1、数数 值值 分分 析析 (8)Numerical AnalysisWenjian Yu2第八章第八章 常微分方程初值问题常微分方程初值问题n(主要是前主要是前3节节)常微分方程初值问题常微分方程初值问题Wenjian Yu3常微分方程基本概念常微分方程基本概念 常微分方程常微分方程nWenjian Yu4常微分方程常微分方程nWenjian Yu5初值问题初值问题:常微分方程常微分方程 例子例子1n含电容元件的电路问题含电容元件的电路问题n通过通过节点分析法节点分析法得到微分方程组得到微分方程组Wenjian Yu6问题的解反映了电容充问题的解反映了电容充/放电过程放电过程C1R2R1R3C2

2、1Ri1Ci2Ri2Ci3Ri1v2v电流电流/ /电压关系电压关系ddCCviCt1RRivR节点电流方程节点电流方程1120CRRiii2320CRRiii112211222232111d0d0111ddvRRRCvtCvvRRRt 1(0)1v2(0)0vt常微分方程常微分方程 例子例子2n双联摆的运动双联摆的运动n两个摆锤两个摆锤(重物重物), 刚性杆的重量可忽略刚性杆的重量可忽略n不考虑摩擦力不考虑摩擦力,运动不会停止运动不会停止, 且在初始角且在初始角度较大时摆锤的轨迹呈现混沌现象度较大时摆锤的轨迹呈现混沌现象Wenjian Yu7求解微分方程初值问题求解微分方程初值问题, 得到得

3、到摆锤的运动规律摆锤的运动规律 Matlab演示演示swinger常微分方程常微分方程nWenjian Yu8线性齐次常系数微分方程线性齐次常系数微分方程实际的问题基本上都是稳定的实际的问题基本上都是稳定的!(由于历史原因由于历史原因)常微分方程常微分方程nWenjian Yu9局部稳定局部稳定简单方法与有关概念简单方法与有关概念Wenjian Yu10简单的初值问题数值解法简单的初值问题数值解法 初值问题的数值解法初值问题的数值解法nWenjian Yu11否则为否则为多步法多步法否则为否则为隐格式方法隐格式方法欧拉法欧拉法nWenjian Yu12“左矩形左矩形”求积公式求积公式h=0.1

4、h=0.050.1 1.000000 1.004837 0.05 1.000000 0.3 1.0350920.2 1.010000 1.018731 0.1 1.002500 0.35 1.0483370.3 1.029000 1.040818 0.15 1.007375 0.4 1.0634200.4 1.056100 1.070320 0.2 1.014506 0.45 1.0802490.5 1.090490 1.106531 0.25 1.023781 0.5 1.098737步长步长h=0.1, 和和0.05步长步长小的更准小的更准数值解法的稳定性数值解法的稳定性nWenjian

5、Yu13-1 0欧拉法解模型问题欧拉法解模型问题的的稳定区域稳定区域数值解法的稳定性数值解法的稳定性nWenjian Yu14-1 0欧拉法稳定欧拉法稳定数值解法的稳定性数值解法的稳定性nWenjian Yu1500.0250.050.0750.10.1250.151-1.52.25-3.3755.0625-7.59375 11.39061 0.082085 0.006738 0.000553 4.5410-53.7310-63.0610-7这里设的这里设的h太大太大!计算结果如下表计算结果如下表:数值解法的局部截断误差数值解法的局部截断误差nWenjian Yu16整体误差整体误差稳定的问题

6、，整体误差小稳定的问题，整体误差小于局部误差之和于局部误差之和不稳定的问题呢？不稳定的问题呢？一般仅能控制局部误差一般仅能控制局部误差整体误差？整体误差？ 局部误差局部误差数值解法的局部截断误差数值解法的局部截断误差nWenjian Yu17欧拉法是欧拉法是一阶方法一阶方法我们讨论的所有方法都至少有我们讨论的所有方法都至少有1阶准确度阶准确度数值解法的数值解法的收敛性收敛性：随着：随着h0，误差，误差0向后欧拉法与梯形法向后欧拉法与梯形法n从数值积分的角度推导从数值积分的角度推导n向后欧拉法向后欧拉法:n梯形法梯形法:n两者均为单步、两者均为单步、隐格式隐格式方法方法, 每每步步计算要求解计算

7、要求解(非线性非线性)方程方程n例例8.6:用向后欧拉法求解用向后欧拉法求解Wenjian Yu18右矩形右矩形梯形梯形00.0250.050.0750.10.1250.1510.006663 0.001904 0.0005441 0.082085 0.006738 0.000553 4.5410-53.7310-63.0610-7向后欧拉法向后欧拉法nWenjian Yu19准确解准确解0 1稳定区域稳定区域无条件稳定无条件稳定(unconditionally stable)!向后欧拉法向后欧拉法nWenjian Yu20具有具有1阶准确度阶准确度!向后欧拉法与梯形法向后欧拉法与梯形法nWe

8、njian Yu21稳定的条件是稳定的条件是: 思考思考无条件稳定无条件稳定!具有具有2阶准确度阶准确度简单方法与有关概念简单方法与有关概念Wenjian Yu22Runge-Kutta方法方法 n在欧拉法基础上改进在欧拉法基础上改进再增加一次函数求值再增加一次函数求值: 数值积分的数值积分的中矩形中矩形或或梯形公式梯形公式利用欧拉法算半个步长的结果利用欧拉法算半个步长的结果, 估算中点处被积函数值估算中点处被积函数值先用欧拉法估计区间终点处斜率先用欧拉法估计区间终点处斜率, 再用它与再用它与起始点斜率的平均值算一整步起始点斜率的平均值算一整步Runge-Kutta方法方法( , ( )n n

9、n nththnnnnt tyyf s y s dsyyf s y s ds1 1中中矩形矩形梯梯形公式形公式(中中点公式点公式)(Heun方法方法/改进的欧拉法改进的欧拉法)都比欧拉法准确都比欧拉法准确均为均为2级级R-K公式公式Wenjian Yu23nRunge-Kutta方法方法( , ( )n nn nththnnnnt tyyf s y s dsyyf s y s ds1 1(, ()r rnninininnininii iyyhc f th y thyyhc f th y th1 11 1()niniy thy th只能估算只能估算() nnnny tyy ty(,) nnnnk

10、f t ykf t y1 1, 令令, , , r rkkkkkk2323r rnni inni ii iyyhc kyyhc k1 11 1Wenjian Yu24nRunge-Kutta方法方法r rnni inni ii iyyhc kyyhc k1 11 1()niniy thy th其他其他, 用所有前面点的信息用所有前面点的信息Wenjian Yu25积分节点积分节点n几种显式几种显式R-K公式公式参数的值不按具体数值参数的值不按具体数值积分公式设置积分公式设置, 而根据而根据准确度阶数准确度阶数要求设置要求设置2级公式：改进欧拉法、级公式：改进欧拉法、中点公式中点公式经典经典4级

11、、级、4阶阶Runge-Kutta法法 (1905)Runge-Kutta方法方法Wenjian Yu26nRunge-Kutta方法方法此时此时局部截局部截断误差断误差只要只要对非模型问题也有相同结论！对非模型问题也有相同结论！Wenjian Yu27如如n例例8.7: 用用2阶改进欧拉、阶改进欧拉、3阶阶Ralston、4阶经典阶经典R-K解问题解问题, h=0.1, 算到算到y(2)精确解为精确解为Runge-Kutta方法方法r4对应的对应的r级级R-K公式公式有有r阶准确度阶准确度高于高于4阶的阶的公式很少单公式很少单独使用独使用Wenjian Yu28r4对应的对应的r级级R-K公

12、式公式达不到达不到r阶阶准确度准确度二阶二阶Heun 三阶三阶Ralston 四阶四阶Runge-Kutta准确值准确值10.40.40.40.41.10.4756410.4746260.47463830.47463821.20.5834080.5813640.58138680.58138671.30.7281350.7250340.72506630.72506621.40.9153290.9111370.91117730.91117711.51.1511101.1457851.14583361.14583331.61.4421691.4356641.43572031.43572001.71.

13、7957381.7880041.78806741.78806711.82.2195782.2105612.21063152.21063111.92.7219612.7116062.71168362.71168322.03.3116653.2999163.30000043.3000000nRunge-Kutta方法方法(此时此时 )一般显式单步法一般显式单步法恰好与欧拉法一样恰好与欧拉法一样显格式显格式, 都都不是不是无条件无条件稳定的稳定的局部局部截断截断误差误差判断单步法判断单步法收敛性收敛性的简便方法的简便方法Wenjian Yu29简单方法与有关概念简单方法与有关概念Wenjian Yu

14、30多步法多步法 n多步法多步法Wenjian Yu31(线性线性m步法步法)固定步长固定步长hnTaylor展开法求展开法求线性线性m步法步法的系数的系数 这些系数应该等于这些系数应该等于0例例8.8: 求两步法公式求两步法公式 中参数值中参数值多步法多步法Wenjian Yu32满足它们才可能满足它们才可能收敛收敛 (相容性相容性), 有有二阶准确度二阶准确度n多步法多步法Wenjian Yu33同同例例8.8的结果的结果!多步法公式中包含多步法公式中包含p个待定参数，至少可达到个待定参数，至少可达到p-1阶准确度阶准确度n多步法多步法Wenjian Yu34Vandermonde阵阵T,

15、 非奇异非奇异插值节点插值节点函数值函数值n常用的多步法公式常用的多步法公式Adams公式公式的推导：的推导：用插值多项式近似被积函数用插值多项式近似被积函数例例8.10: 推导推导m=4对应的显式对应的显式Adams公式公式多步法多步法Wenjian Yu35可证明满足最高准确度阶数可证明满足最高准确度阶数类似地算其类似地算其他系数，得他系数，得显式四阶显式四阶Adams-Bashforth公式(单项式函数代入法单项式函数代入法)nAdams公式公式几种显式公式几种显式公式几种隐式公式几种隐式公式多步法多步法Wenjian Yu36阶数阶数稳定阈值稳定阈值 误差常数误差常数11 -21/22

16、3/2-1/2 -15/12323/12 -16/12 5/12 -6/113/8455/24 -59/24 37/24 -9/24-3/10251/720欧拉法欧拉法阶数阶数稳定阈值稳定阈值 误差常数误差常数11 -1/221/21/2 -1/1235/128/12-1/12 -6-1/2449/2419/24 -5/241/24-3-19/720向后欧拉法向后欧拉法梯形法梯形法并非无条件并非无条件稳定稳定!n多步法多步法Wenjian Yu37用用Matlab求解初值问题求解初值问题Wenjian Yu38用用Matlab解解ODE-IVP Matlab中的中的ODE-IVP求解器求解器n

17、Wenjian Yu39T, Y, TE, YE, IE = solver(odefun, tspan, y0, options) 求解单个求解单个ODEn火焰燃烧问题火焰燃烧问题当当点燃一根火柴点燃一根火柴时时, 火焰火焰迅速增大直到一个迅速增大直到一个临界体积临界体积, 然然后后维持这一体积维持这一体积不变不变, 此时此时火焰内部燃烧耗费的氧气和火焰内部燃烧耗费的氧气和其表面现存的氧气达到了一种其表面现存的氧气达到了一种平衡平衡.火焰火焰(近似为球近似为球)半径半径y满足满足ODE f= (t,y) y2-y3; ode23(f, 0, 2.0e4, 1e-4)Wenjian Yu4023

18、-(0)yyyy 设初始半径设初始半径=0.000100.511.52x 10400.20.40.60.811.21.4例例8.15尝试用尝试用ode23s求解它求解它00.10.20.30.40.50.60.70.80.9100.511.522.533.544.5求解求解ODE方程组方程组nWenjian Yu41例例8.1612122 , (0)0,(0)16yyyyyt function ydot = myode2(t, y);ydot= y(2); 6*t; %列向量列向量双联摆问题的求解双联摆问题的求解Wenjian Yu42隐格式非线性常微分方程初值隐格式非线性常微分方程初值问题问题, 初值为初值为, , 0, 0Todeset可以设置质量矩阵可以设置质量矩阵质量矩阵质量矩阵swinger_solve.m很大时混沌现象很大时混沌现象: