1、化工应用数学第三章 第三章第三章 常微分方程的数值解法常微分方程的数值解法主主 要要 内内 容容 3.1、引言引言 3.2、初值问题、初值问题 3.3、边值问题、边值问题 12/3/20221化工应用数学第三章 3.1、引引 言言1.1.什么是微分方程什么是微分方程?在化学工程中关于扩散、反应、传质、传热在化学工程中关于扩散、反应、传质、传热和流体流动等问题的数学模型。和流体流动等问题的数学模型。微分方程：微分方程：是包含一个未知函数及其导数关是包含一个未知函数及其导数关系的方程。系的方程。常微分方程：常微分方程：其中只包含一个自变量的导数其中只包含一个自变量的导数的方程。的方程。12/3/2
2、0222化工应用数学第三章 123123123:pPAsSATTAAAAAAAAdck cdtdcrk cdtdcrk cdtdckkkcdtdckkkccdt r r对于 r求解123kkkAPASAT 求解这样的方程一般我们需要给定初值。求解这样的方程一般我们需要给定初值。例题：例题：12/3/202231230()AAkkktcec利用以前所学知识，进行积分可得：利用以前所学知识，进行积分可得：化工应用数学第三章 01d()00,()xdTkdxdxdTT LTdx 由以上两例可以看出，初值问题和边值问题区别在于：前者由以上两例可以看出，初值问题和边值问题区别在于：前者在自变量一端给定附
3、加条件，后者在自变量两端附加条件。在自变量一端给定附加条件，后者在自变量两端附加条件。例题：例题：一维均匀介质稳态导热问题。设其一端绝热，另一端恒温一维均匀介质稳态导热问题。设其一端绝热，另一端恒温为为T1，则此问题模型可用一常微分方程边值问题来描述。则此问题模型可用一常微分方程边值问题来描述。12/3/20224化工应用数学第三章 3.2、初值问题、初值问题1.尤尤拉法（拉法（Euler）2局部截断误差局部截断误差3改进尤拉法改进尤拉法4.龙格龙格-库塔法库塔法5.常微分方程组常微分方程组 6.步长的选择步长的选择7.收敛性和稳定性收敛性和稳定性8.线性多步法线性多步法12/3/20225化
4、工应用数学第三章 00(,)(1.1)()(1.2)yf x yaxby xy一般初值问题可表示为 其中f(x,y)是已知函数，(1.2)是定解条件也称为初值初值条件。条件。12/3/202261.尤拉尤拉(Euler)法法化工应用数学第三章 1.尤拉尤拉(Euler)法法考虑一阶初值问题考虑一阶初值问题：0(,)()yf x yaxby ay 通过欧拉方法的讨论通过欧拉方法的讨论12/3/20227化工应用数学第三章 即等距，要计算出解函数即等距，要计算出解函数y y(x x)在一系列节点在一系列节点()/iiixaihbnhah,一般取节点节点处的近似值处的近似值01 .naxxxb把区间
5、把区间a,b 分为分为n个小区间个小区间步长为1(-)iiihxxn n 等份等份12/3/2022801,nu uu1.尤拉尤拉(Euler)法法化工应用数学第三章 假定是假定是 y(x)初值问题的解，那么，把在节点初值问题的解，那么，把在节点 xi 附近展开附近展开成成Taylor级数，则级数，则)(2)()()(21yhxyhxyxyiii)(2)(,()(2yhxyxhfxyiii 0(,)()yf x yaxby ay 12/3/20229()iiuy x取前两项，并令取前两项，并令1(,)iiiiuuh fxu1.尤拉尤拉(Euler)法法化工应用数学第三章 100021111(,
6、)(,).(,)nnnnuuhf xuuuhf x uuuhf xu依上述公式逐次计算可得：依上述公式逐次计算可得：每步计算每步计算只用到只用到1nunu12/3/2022101(,)iiiiuuhfxu1.尤拉尤拉(Euler)法法化工应用数学第三章 几何意义几何意义0000111222(,)()(,)()PPxyy xxyy xuxxuxx 过 点 的 曲 线 是 解在 作的 切 线(斜 率)与 直 线交 于点再 作 切 线（斜 率）与 直 线交 于点12/3/2022110 x1x2x0P1P2P()yy x化工应用数学第三章.1)0(),10(2yxyxyy解解 取步长h=0.1，欧拉
7、公式的具体形式为)2(1nnnnnyxyhyy 已知y0=1,由此式可得例题：例题：用尤拉法解下列方程组，取h=0.1。.191818.1)1.12.01.1(1.01.1)2(1.11.01)2(1111200001 yxyhyyyxyhyy12/3/2022121.尤拉尤拉(Euler)法法化工应用数学第三章 xy21 与准确解与准确解 相比，可看出欧拉公式的计算结相比，可看出欧拉公式的计算结果误差较大果误差较大.xn 尤拉公式数值解yn准确解y(xn)误差0.20.40.60.81.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.1832
8、16 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.05271912/3/2022131.尤拉尤拉(Euler)法法化工应用数学第三章 定义定义若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。尤尤拉法的局部截断误差：拉法的局部截断误差：222()()hiyxO h因此，尤拉法具有因此，尤拉法具有 1 阶精度。阶精度。),()()(),()()(2112iiiihiiiiiiyxhfxyxyyxhfxyuxyR 定义定义 在假设 ,即第i步计算是精确的前提下，考虑的截断误差 称为
9、局部截断误差。称为局部截断误差。2.局部截断误差局部截断误差12/3/202214()iiuy x11()iiiRy xu化工应用数学第三章 显式尤拉公式显式尤拉公式将将y(xi+1)=y(xi+h)在在 处处Taylor展开展开隐式尤拉公式隐式尤拉公式将将y(xi+1)=y(xi+h)在在 处处Taylor展开展开显式法计算较简单，隐式法计算复杂，可以编程。显式法计算较简单，隐式法计算复杂，可以编程。3.改进的尤拉方法改进的尤拉方法12/3/202215化工应用数学第三章 隐式尤拉公式隐式尤拉公式 显式显式与与隐式隐式两类方法各有特点，考虑到数值稳两类方法各有特点，考虑到数值稳定性等其他因素
10、，人们有时需要选用定性等其他因素，人们有时需要选用隐式隐式方法，但方法，但使用使用显式显式算法远比算法远比隐式隐式方便。方便。隐式方程通常用迭代法求解，而迭代过程的实隐式方程通常用迭代法求解，而迭代过程的实质是质是逐步逐步显式化显式化。12/3/2022163.改进的尤拉方法改进的尤拉方法化工应用数学第三章 1111(,)(,)2(,)iiiiiiiiiihuuf x yf xyuuhf x y尤拉法尤拉法隐式尤拉法隐式尤拉法12/3/202217两式相加可得改进的尤拉公式，通常将其与显示尤拉公式联用两式相加可得改进的尤拉公式，通常将其与显示尤拉公式联用3.改进的尤拉方法改进的尤拉方法化工应用
11、数学第三章 改进尤拉法改进尤拉法/*modified Eulers method*/Step 1:先用先用显式尤显式尤拉公式作拉公式作预测预测，算出，算出Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy)1,.,0(),(,),(211niyxfhuxfyxfhuuiiiiiiii12/3/2022181(,)iiiiuuhf x y111(,)(,)2iiiiiihuf xyuf xy改进的尤拉公式是二阶方法，但计算量增加了一倍。改进的尤拉公式是二阶方法，但计算量增加了一倍。3.改进的尤拉方法改进的尤拉方法化工应用数学第三章 10,0.5(
12、0)10.1yxyxyh 例：110.110.1(1)12iiiiiiiiyyxyxyxy 12/3/202219解：由改进的尤拉公式可得解：由改进的尤拉公式可得111(,)(,)2iiiiiihuf x yuf xy 将数值代入，进行逐步计算。将数值代入，进行逐步计算。一般对于改进的尤拉法可采用计算机编程进行计算。一般对于改进的尤拉法可采用计算机编程进行计算。化工应用数学第三章 考察改进的尤拉法，可以将其改写为：考察改进的尤拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 基本思想基本思想:是从是从(xi,yi)点出发，以点出发，以某一斜某一
13、斜率率沿直线沿直线达到达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为点。欧拉法及其各种变形所能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度的单步递推格式。4.4.龙格龙格-库塔法库塔法12/3/202220化工应用数学第三章 首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyx
14、phfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进尤拉法推广为：将改进尤拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2:将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 12/3/202221化工应用数学第三章 Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点的点的泰勒泰
15、勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式。塔格式。21,121 p注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？12/3/202222化工应用数学第三章 其中其中
16、i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKK
17、KKyyiihihihihiiihii 12/3/202223化工应用数学第三章 11234(22)6iihyyKKKK例题例题1211322433(,)1(,)()()12222(,)()()12222(,)()()1iiiiiiiiiiiiiiiiKf x yxyhhhhKf xyKxyKhhhhKf xyKxyKKf xh yhKxhyhK解：12/3/2022241(0)1yxyy 0,0.50.1xh12340,0.05,0.0475,0.09525KKKK将K值代入式，即可计算出方程组的解。化工应用数学第三章 注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki
18、 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n12/3/202225化工应用数学第三章 0000()(,)()()(,)()y xf x y zy xyz xg x y zz xz，1123411234(22)6(22)6iiiihyykkkkhzzLLLL11(,)(,)iiiiiikf xyzLg xyz 211211(,)2
19、22(,)222iiiiiihhhkf xykzLhhhLg xykzL322322(,)222(,)222iiiiiihhhkf xykzLhhhLg xykzL 433433(,)(,)iiiiiikf xh yhk zhLLg xh yhk zhL5.常微分方程组常微分方程组12/3/202226化工应用数学第三章 5.常微分方程组常微分方程组12/3/202227高阶常微分方程0001()(,)()()yxf x y yy xyyxy 思路：转化降阶，通常转化为常微分方程组（思路：转化降阶，通常转化为常微分方程组（1阶）来求解。阶）来求解。0100()(,)()=()()z xf x
20、y zz xyy xzy xy 令 z=y化工应用数学第三章 222sin(0)0.4(0)0.6xyyyexyy 例 题=z y令 则方程转化为2 22sin(0)0.6=(0)0.4xzzyexzz yy 1123411234(22)6(22)6iiiihyykkkkhzzLLLL12/3/202228222sinxyzzzy ex 化工应用数学第三章 121sin22iixiiikzLexzy212()221121sin()2()2()222iihxiiihkzLhhLexzhLyk322()232221sin()2()2()222iihxiiihkzLhhLexzhLyk432()43
21、32sin()2()2()iixhiiihkzLLexhzhLyhk12/3/202229化工应用数学第三章 6.步长的选择步长的选择利用 p 阶方法可构造出：步长为 h:步长为 h/2:（1）（2）上述两式整理可得：（3）12/3/202230 从局部截断误差上来看，与 有关，故 越小，对降低局部误差有利，但是步长越小，势必引起计算量的加大。1()pO hh121112211()()()2()()2hppiihppiiy xychO hhy xycO h2111221hphiiipyyy化工应用数学第三章 2111(1)(2)=ihhiiyyy 令对于给定的精度，如果，可以将步长反复减半进行
22、计算，直到为止，这时以最终得到的作为结果；如果，可以反复将步长加倍，直到 ,这时将步长再减半一次，即为所求结果。21121(1)21hhhiiiipyyy xy 将(3)式变形为：12/3/2022316.步长的选择步长的选择化工应用数学第三章 7.收敛性和稳定性收敛性和稳定性11y()iihy xh 当步长 趋于零时，方程的数值解趋近于精确解时，则称所采用的数值方法是收敛的。当步长 确定之后，随着步数的增加，计算中累积的误差不会超出所允许的范围。收敛性：收敛性：稳定性：稳定性：12/3/202232 收敛性与稳定性从两个不同的角度描述了微分方程数值解法收敛性与稳定性从两个不同的角度描述了微分
23、方程数值解法的实用价值，只有既收敛又稳定的方法，才可以提供比较可靠的实用价值，只有既收敛又稳定的方法，才可以提供比较可靠的计算结果。的计算结果。化工应用数学第三章 11(,)(,)+1fiixiixdyf x yydxyyf x y dxk对于如果用阶的牛顿后插公式代替被积函数（x,y)则012012-,.,.yy yy y龙格 库塔法和尤拉法均为单步法，所谓单步法，即只要给出初始即可顺序地算出。而多步法则除了给出y 外尚需借助其他单步法提供等若干个点。21210(1)(1).(1)().2!(1)(1).(1).(0,1,2.)2!kiiiikiiiiiip pp ppky xyp yyyk
24、p pp ppkyyhyp yyy dp ikk11110(,)13(,)(,)2iiiiiiiiiihkyyhf x ykyyf x yf xy 8.8.线性多步法线性多步法12/3/202233化工应用数学第三章 n1.打靶法打靶法 n2差分法差分法3.3、边值问题、边值问题12/3/202234化工应用数学第三章 边值问题的数值解边值问题的数值解 /*Boundary-Value Problems*/2 阶常微分方程边值问题阶常微分方程边值问题 )(,)(),(),(byaybaxyyxfy 打靶法打靶法/*shooting method*/先猜测一个初始斜率先猜测一个初始斜率 y (a
25、)=m，通过解初，通过解初值问题值问题 mayaayyyxfy)()(),(y(b)=(s)找出找出s*使得使得(s*)=，即把问，即把问题转化为求方程题转化为求方程 (s)=0 的根。的根。yx0aby x()斜率斜率=s0()s0斜率斜率=s1()s112/3/202235化工应用数学第三章 184(0)0,(10)0yyyy 例题：例题：用打靶法解。用打靶法解。12/3/2022361.打靶法打靶法 化工应用数学第三章 184(0)0,(0=yyyym）001110-(10).11,(10)myAmyB首先取初值 ，利用四阶龙格 库塔法求出 求出012.5212,11 10(10)1.0
26、8853 10m mmBmmyBA然后利用用弦截法求出下一个初值则（）求出3.7442731.82644412/3/2022371.打靶法打靶法 化工应用数学第三章 有限差分法有限差分法基本思想：基本思想：运用数值微分将导数用离散点上函数值表示，从而将边运用数值微分将导数用离散点上函数值表示，从而将边值问题的微分方程和边界条件转化为只含有限个未知数的差分方值问题的微分方程和边界条件转化为只含有限个未知数的差分方程组，并将此差分方程组的解作为该边值问题的数值解。程组，并将此差分方程组的解作为该边值问题的数值解。1.二阶常微分方程的第一边值问题二阶常微分方程的第一边值问题 其中其中q(x)(0)，
27、f(x)在在a，b)上连续，上连续，为常数。为常数。()()()(),(),()y xq x y xf xax by ay b 12/3/2022382 差分法差分法化工应用数学第三章 n设等距节点：设等距节点：xi=a+ih，i=0，1，2，n，n对其中内节点应用三点微分公式：对其中内节点应用三点微分公式：n(当当h充分小时，略去充分小时，略去O(h2)，并以，并以yi-1，yi，yi+1，代，代y(xi-1)，y(xi)，y(xi+1)，得计算，得计算yi的差分方程组的差分方程组).()()(2)()(2211hOhxyxyxyxyiiii)()()()()()(2)(2211hOxfxy
28、xqhxyxyxyiiiiii 1122iiiiiiyyyq yfh12/3/2022392 差分法差分法化工应用数学第三章 加上边界条件即得边值问题的差分方程组加上边界条件即得边值问题的差分方程组其中其中qi=q(xi)，f i=f(xi)。矩阵形式：矩阵形式：niiiiiyynifhyyqhy,1,2,1,)2(02121 1222121210122212.11211211211nnnnfhfhfhyyyyyqhqhqh12/3/2022402 差分法差分法化工应用数学第三章 注：注：可以证明三角方程组存在唯一解，且有误差估计，其中可以证明三角方程组存在唯一解，且有误差估计，其中 三角方程
29、组等价于：三角方程组等价于：n为对角占优的三对角方程组，可用追赶法求解为对角占优的三对角方程组，可用追赶法求解。.1,2,1,96)()(224 nihabMyxyii)(max)4(4xyMbxa 12222212122112222212.21121.12112nnnnnnfhfhfhfhyyyyqhqhqhqh412 差分法差分法化工应用数学第三章 例例 取取h=0.1，求边值问题：，求边值问题：的数值解。的数值解。解：本例的节点解：本例的节点 xi=0.1i，i=0，1，2，10，q(x)=1，r(x)=x，差分方程为差分方程为 1)1(,0)0(10,)()(yyxxxyxy 1,01,2,1,1.0)1.02(02121niiiiyynixyyy12/3/20224222110(2),iiiinyhyyh xyy化工应用数学第三章 1,01,2,1,01.001.2011niiiyyniiyyy 991.0008.0.002.0001.001.21101.21.101.21101.21221nnyyyy12/3/202243解此方程组。解此方程组。差分方程为差分方程为
