高阶常微分方程的数值求解.doc

高阶常微分方程的数值求解 谷照升 （长春工程学院理学院，长春，130012） 摘 要 对经典初始条件的高阶常微分方程，给出其数值求解方法。该方法比Runge-Kutta法具有更好的适应性、易用性、计算速度和可控制的更高精度。 关键词 常微分方程；数值解；算法 中图分类号: O241.81　　　文献标识码: A 1 引言 求解复杂的1阶常微分方程，通常只能采用数值解法。数值解法一般又以Runge-Kutta法为主。对高阶常微分方程，则通常是将其转化为1阶常微分方程组，再用Runge-Kutta法求解[1, 2]。但这种通用性方法，在精度、软件计算的适应度方面，常常不够理想，甚至得不到结果。针对不同的广普性方程类型，可以建立更具针对性的计算方法。例如，针对三类边值条件和特定形式的方程，已经有有相应的差分法和有限元法。每一种更具针对性的方法，都有其更高的精度和更为健壮的算法，当然也存在其必然的局限性。 本文对如下形式的2阶常微分方程 x∈[a, b] (1) 给出了完整、方便、高效的单步法数值求解算法，并将其推广到任意高阶问题的求解。 2 算法 思路：将[a, b]分割为a=x0 x1…xn=b，步长。从i=1开始，利用数值积分公式，通过[xi-1, xi]上的数值积分，先求，再求出，最后求y(xi)。依次取i=2, 3, … ,n，得到各点、、y(xi)的近似值。 根据数值积分的算法不同，主要有两种不同的计算公式。 2.1、矩形数值积分算法 矩形算法形式简单，精度偏低。基于不同的积分形式，又可以得到3种不同的计算公式，其精度无显著差别。此处仅给出与2.2梯形算法最为接近的一种。 利用 在点得到 (2) (3) 这里 。 将(2) 、(3)代入(1), 在点得到 从而 (4) 算法： Step 1: i=1; 将(1)中 、代入 (4), 输出; 将代入 (2), 输出; 取, 输出。 Step 2: while in { i= i+1; 将、代入 (4), 输出; 将代入 (2), 输出; 取, 输出; } 2.2、梯形数值积分算法 根据(1)式，首先算出 (5) 对, 取 则点满足 (6) (7) 将(6) 、(7)代入(1), 得到 从而 (8) 算法： Step 1: 按(5)式求出。 Step 2: i=0; while in { i= i+1; 将、、代入 (8), 输出; 将、代入 (6), 输出; 将、、代入 , 输出; } 3 误差分析 记h=(实际计算中，步长通常是取定的，不随i变化)，在矩形算法中，根据(3)式，误差的阶为O(h3)。在梯形算法中，根据(7)式，误差的阶仍为O(h3)。虽然两种方法误差的阶是相同的，但由于在数值微积分中，没有可以确定的绝对误差与相对误差，所以仅凭误差的阶还不能完全说明实际误差水平。本文算法对h、、、是敏感的，像所有单侧初始条件的微分方程数值算法一样，其误差随求解区间的增长、步数的增加，会产生传递和累进。所以，控制误差的基本原则，是尽可能选择初始点和步进方向，让敏感因素的绝对值尽可能从小到大变化（参看后文）。通过多个计算实例的统计比较发现，对同水平的h，梯形算法的相对误差通常远远小于矩形算法相对误差的一半。由于梯形算法和矩形算法的运算速度没有显著差别，建议尽可能采用梯形算法。 4 算例 通过一个比较复杂的方程 先取，由于函数2阶导比较稳定，所以误差太小。 又取，得到方程 在 [0, 100] 内根据其在x=0的初始条件进行梯形算法数值求解，并与精确解进行比较，结果如下（图1、2、3）。 图1 的数值解、精确解、绝对误差、相对误差 图2 的数值解、精确解、绝对误差、相对误差 图3 的数值解、精确解、绝对误差、相对误差 可以看出，即使在x=100点，y(100)的近似值487187.30529与准确值487187.675007误差也仅为0.369715682，其1阶、2阶导数值的绝对误差也都在0.6以内，相对误差则在10-6之内。 对笔者而言，测试结果有些意外，因为相对误差太小，计算时间不到5秒。取步长Δx=0.005测试，、、的绝对误差分别为3.902、14.73和12.003。 与其相对照的是，使用最权威的专业数学软件Maple进行同样的计算，它只能在区间[0, 46.839703] 内求解，当让它求解y(100)时，就报错了（图4）。它算得y(46) = -60977.64，=147763.772594。用本文的方法得到y(46