常微分方程初值问题数值解的实现和分析—四阶Rungekutta方法与预估校正算法毕业论文

合集下载

常微分方程初值问题的解法及应用

常微分方程初值问题的解法及应用常微分方程是数学中非常重要的一部分，它涉及了许多领域的模型建立和问题求解。

本文将介绍常微分方程初值问题的解法及其应用。

一、常微分方程初值问题的定义常微分方程初值问题是指给定一个常微分方程，以及它在某一点上的初始条件，求解该方程的解曲线。

通常，一个常微分方程初值问题可以表示为:y'(x) = f(x,y), y(x0) = y0,其中，y(x)是未知函数，f(x,y)是已知函数，y(x0) = y0是初始条件。

二、常微分方程初值问题的解法常微分方程初值问题的解法有多种，下面我们将介绍几种常用的方法。

1.欧拉法欧拉法是最简单的一种求解常微分方程初值问题的方法。

该方法基于初始条件，通过不断迭代计算得到近似解曲线。

具体步骤如下：步骤1：设定步长h，确定求解区间[x0, xn]，计算步数n。

步骤2：初始化，即确定初始点(x0, y0)。

步骤3：根据方程dy/dx = f(x,y)和初始点(x0, y0)，计算斜率k = f(x0, y0)。

步骤4：根据已知的斜率和步长h，计算下一个点的坐标(xi+1,yi+1)。

步骤5：重复步骤3和步骤4，直到达到步数n。

步骤6：得到近似解曲线。

2.改进的欧拉法（改进欧拉法）改进的欧拉法是对欧拉法的改进，其求解精度比欧拉法更高。

具体步骤如下：步骤1：设定步长h，确定求解区间[x0, xn]，计算步数n。

步骤2：初始化，即确定初始点(x0, y0)。

步骤3：根据方程dy/dx = f(x,y)和初始点(x0, y0)，计算斜率k1 =f(x0, y0)。

步骤4：根据已知的斜率k1和步长h/2，计算中间点的坐标(x0+h/2, y0+k1*h/2)。

步骤5：根据方程dy/dx = f(x,y)和中间点的坐标(x0+h/2, y0+k1*h/2)，计算斜率k2= f(x0+h/2, y0+k1*h/2)。

步骤6：根据已知的斜率k2和步长h，计算下一个点的坐标(xi+1,yi+1)。

四阶Runge-Kutta方法

实验题目3 四阶Runge-Kutta 方法摘要一阶常微分方程初值问题0(,)()dy f x y dxy x y ⎧=⎪⎨⎪=⎩ （6.1）的数值解法是近似计算中很重要的部分。

常微分方程初值问题的数值解法是求方程（6.1）的解在点列1(0,1,)n n n x x h n -=+=L 上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。

常微分方程初值问题的数值解法一般分为两大类：（1）单步法：这类方法在计算n y 时，只用到1n x +、n x 和n y ，即前一步的值。

因此，在有了初值以后就可以逐步往下计算。

典型方法如龙格–库塔()R K -方法。

（2）多步法：这类方法在计算1n y +时，除用到1n x +、n x 和n y 以外，还要用(1,2,,;0)n p y p k k -=>L ，即前面k 步的值。

典型方法如Adams 方法。

经典的R K -方法是一个四阶的方法，它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （6.2） R K -方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。

前言利用四阶龙格-库塔方法求解微分方程的初值问题程序设计流程龙格-库塔法流程图问题1（1）TestRK4('ode1', 1, [0 -1], 5, inline('-x-1'))TestRK4('ode1', 1, [0 -1], 10, inline('-x-1'))TestRK4('ode1', 1, [0 -1], 20, inline('-x-1'))（2）TestRK4('ode2', 1, [0 1], 5, inline('1./(x+1)'))TestRK4('ode2', 1, [0 1], 10, inline('1./(x+1)'))TestRK4('ode2', 1, [0 1], 20, inline('1./(x+1)'))问题2（1）TestRK4('ode3', 3, [1 0], 5, inline('x.^2.*(exp(x)-x)'))TestRK4('ode3', 3, [1 0], 10, inline('x.^2.*(exp(x)-x)'))TestRK4('ode3', 3, [1 0], 20, inline('x.^2.*(exp(x)-x)'))（2）TestRK4('ode4', 3, [1 -2], 5, inline('2*x./(1-2*x)'))TestRK4('ode4', 3, [1 -2], 10, inline('2*x./(1-2*x)'))TestRK4('ode4', 3, [1 -2], 20, inline('2*x./(1-2*x)'))问题3（1）TestRK4('ode5', 1, [0 1/3], 5, inline('x.^2+1/3*exp(-20*x)'))TestRK4('ode5', 1, [0 1/3], 10, inline('x.^2+1/3*exp(-20*x)'))TestRK4('ode5', 1, [0 1/3], 20, inline('x.^2+1/3*exp(-20*x)'))（2）TestRK4('ode6', 1, [0 1], 5, inline('exp(-20*x)+sin(x)'))TestRK4('ode6', 1, [0 1], 10, inline('exp(-20*x)+sin(x)'))TestRK4('ode6', 1, [0 1], 20, inline('exp(-20*x)+sin(x)'))（3）TestRK4('ode7', 1, [0 0], 5, inline('exp(x).*sin(x)'))TestRK4('ode7', 1, [0 0], 10, inline('exp(x).*sin(x)'))TestRK4('ode7', 1, [0 0], 20, inline('exp(x).*sin(x)'))实验所用函数function [x,y] = RK4ODE(fun, xEnd, ini, h)% RK4ODE 用四阶Runge-Kutta法解初值问题dy/dx = f(x,y),y(x0) = y0，在x处y的值%% Synopsis: [x,y] = RK4ODE(fun, xEnd)% [x,y] = RK4ODE(fun, xEnd, ini)% [x,y] = RK4ODE(fun, xEnd, ini, h)%% Input: fun = (string) 初值问题的函数% xEnd = 使用Euler法的截止点% ini = (optional)初始条件[x0 y0]，默认为[0 0]% h = (optional)步长，默认为0.05%% Output: y = 初值问题在x处y的近似值if nargin < 3ini = [0 0]; %若未给初始条件，将初始条件设为[0 0]endif nargin < 4h = 0.05; %若未给出步长，将步长设为0.05endini = ini(:); %将初始条件转为列向量，便于判断是否正确[m,n] = size(ini);if m ~= 2 | n~= 1error('初始值必须是一个含两个元素的向量[x0 y0]');endx0 = ini(1); %初始化x0y0 = ini(2); %初始化y0x = (x0:h:xEnd)'; %构建x向量y = y0*ones(length(x), 1); %初始化y向量for j=2:length(x)k1 = h * feval(fun, x(j-1), y(j-1)); %三阶Runge-Kutta法的递推公式：y(n+1) = y(n) + (k1 + 2*k2 + 2*k3 + k4) / 6k2 = h * feval(fun, x(j-1)+h/2, y(j-1)+k1/2); % k1 = h * f( x(n), y(n) )k3 = h * feval(fun, x(j-1)+h/2, y(j-1)+k2/2); % k2 = h * f( x(n)+h/2, y(n)+k1/2 )k4 = h * feval(fun, x(j-1)+h, y(j-1)+k3); % k3 = h * f( x(n)+h/2, y(n)+k2/2 )y(j) = y(j-1) + (k1+2*k2+2*k3+k4)/6; % k4 = h * f( x(n)+h, y(n)+k3 )endfunction TestRK4(fun, xEnd, ini, N, result)h = (xEnd - ini(1))/N;[x,y] = RK4ODE(fun, xEnd, ini, h);y0 = feval(result, x);plot(x,y,'ro',x,y0,'b-');function dydx = ode1(x,y)dydx = x+y;function dydx = ode2(x,y)dydx = -y.^2;function dydx = ode3(x,y)dydx = 2*y./x + x.^2.*exp(x);function dydx = ode4(x,y)dydx = (y.^2 + y)./x;function dydx = ode5(x,y)dydx = -20*(y-x.^2) + 2*x;function dydx = ode6(x,y)dydx = -20*y + 20*sin(x) + cos(x);function dydx = ode7(x,y)dydx = -20*(y - exp(x).*sin(x)) + exp(x).*(sin(x) + cos(x));思考题1、实验一中(1)数值解和解析解相同，(2)数值解和解析解稍有不同，因为四阶Runge-Kutta方法是以小段的线性算法来近似获得微分方程的数值解，(1)的准确解是1阶的，(2)的准确解是无限阶的，因此对于(1)数值解和解析解相同。

龙格-库塔法,求解常微分方程

隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用Matlab编程进行计算求解，为了体现计算结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为比较.通过比较三种方法的计算精度，发现四阶经典龙格-库塔方法的误差最小，预估较正法其次，欧拉方法误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格库塔法求解常微分方程初值问题的计算精度的影响.总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度也最高.关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；MatlabRunge Kutta Method For Solving Ordinary Differential EquationsABSTRACTProblem solving ordinary differential equations are often needed in science andtechnology. the problem in the simplest form is the initial value problem of first order equations in this chapter ,which will be discussed. Although there are various analytical methods for solving ordinary differential equations, the analytical method can only be used to solve some special types of equations.differential equations can be summed up the actual problems whichThis paper discusses the initial value problem of Runge Kutta Barclays by solving a differential equation, using the four order Runge Kutta method with high accuracy.for instance through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correction method for instance by calculation through the calculation precision. The comparison of three kinds of methods, found that the error of four order Runge Kutta method of minimum, prediction correction method secondly, Euler method error is relatively large. Finally, by selecting different step, study the affect the calculation accuracy of different step of Runge Kutta method to solve initial value problems of ordinary differential equations.In short, this paper comprehensively analyzes the application of Runge Kutta method for solving ordinary differential equations, compared with the numerical solution of other, higher accuracy Runge Kutta method, good convergence, the Runge Kutta method of order four of the highest efficiency and its precision is the highest.Key words: Four order Runge Kutta method; Euler method; prediction correction method; first order ordinary differential equations; Matlab目录1 问题的提出 (1)1.1 问题背景............................................... . (1)1.2 问题的具体内容 (1)2 问题假设 (2)3 符号系统 (2)4 问题的分析 (3)4.1 欧拉格式 (3)4.2 预估较正法 (3)4.3 四阶隆格库塔法的格式 (4)5 模型的建立与求解 (4)5.1 隆格库塔法的基本原理 (4)5.1.1 Taylor级数 (4)5.1.2 隆格库塔法的基本思想 (4)5.1.3 四阶的隆格库塔法 (5)5.2 其他求解常微分方程边值问题算法的简介 (6)5.3 模型求解 (8)5.3.1 运用MATLAB软件对模型求解结果及析 (8)6 模型的评价 (16)7 课程设计的总结与体会 (16)参考文献 (17)附录 (18)一、问题的提出1.1 问题背景：科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩ （1）其中a ，b 为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重讨论一类高精度的单步法——隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较.1.2 问题的具体内容实例一：在区间[0,1]上采用经典的四阶隆格库塔方法求解微分方程1(0)1dy y x dx y ⎧=-++⎪⎨⎪=⎩，其精确解为x y x e -=+，步长为0.5,然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.实例二：在区间[0,1]上用经典的四阶龙格库塔方法求解初值问题 (0)1x dy xe y dx y -⎧=-⎪⎨⎪=⎩，其精确解为21(2)2xx e -+，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响.二、问题假设2.1 假设数值方法本身的计算是准确的.2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解.2.3 假设步长的增加不会导致舍入误差的严重积累.三、符号系统3.1 符号说明符号含义h选取的步长*K平均斜率p精度的阶数∆前后两次计算结果的偏差y第n个节点的实验值n()y x第n个节点的精确值nδ实验值与精确值的绝对误差四、问题的分析问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优越性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，分析在选取不同的步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及经典的四阶隆格库塔法的计算公式.4.1欧拉格式（1）显式欧拉格式1(,)n n n n y y hf x y +=+ （2）局部截断误差22211()()()()22n n n h h y x y y y x o h ξ++''''-=≈= （3）（2）隐式欧拉格式111(,)n n n n y y hf x y +++=+ （4）局部截断误差2211()()()2n n n h y x y y x o h ++''-≈-= (5) 4.2 预估校正法预估 1(,)n n n n y y hf x y +=+ (6)校正 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ (7) 统一格式 1[(,)(,(,))]2n n n n n n n n h y y f x y f x h y hf x y +=++++ (8) 平均化格式 11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧⎪=+⎪=+⎨⎪⎪=+⎩ (9)4.3 四阶龙格库塔方法的格式（经典格式）112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩(10)五、模型的建立与模型求解5.1 隆格库塔法的基本原理隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联系.5.1.1 Taylor 级数设初值问题 '00(,)()y f x y y x y ⎧=⎨=⎩有解，按泰勒展开，有2'''1()()()()....2n n n n h y x y x hy x y x +=+++；（11）其中()y x 的各阶导数依据所给方程可以用函数f 来表达，下面引进函数序列(,)j f x y 来描叙求导过程，即(0)(1)'(0)''(1)(1)(1)'''(2),f f y f f y f f x x f f y f f x y ∂∂=≡=+≡∂∂∂∂=+≡∂∂ （12）(2)(2)()(1)j j j j f f y f f x y ---∂∂=+≡∂∂ （13）根据上式，结果导出下面 Taylor 格式2'''()1...2!!pp n n nn n h h y y hy y y p +=++++ （14）其中一阶Taylor 格式为： '1n n n y y hy +=+提高Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶Taylor 格式的局部截断误差为：11(1)1(1)!p p n n h y y y p ζ+++-+=+ （15）因此它具有p 阶精度.5.1.2 隆格库塔方法的基本思想隆格库塔法实质就是间接地使用Taylor 级数法的一种方法,考察差商1()()n n y x y x h+- 根据微分中值定理，这有'1()()()n n n y x y x y x h h θ+-=+ （16）利用所给方程 '(,)y f x y =1()()(,())n n n n y x y x hf x h y x h θθ+=+++ （17）设平均斜率*(,())n n K f x h y x h θθ=++，由此可见，只要对平均斜率一种算法，便相应地可以导出一种计算格式.再考察改进的Euler 格式，它可以改写成平均化的形式：1121211()2(,)(,)n n n n n n h y y K K K f x y K f x y hK ++⎧=++⎪⎪=⎨⎪=+⎪⎩（18）这个处理过程启示我们，如果设法在1(,)n n x x +内多预测几个点的斜率值，然后将它们加权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基本思想.5.1.3 四阶的隆格库塔法为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如下：112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩（19）下面为其具体的算法流程图：5.2 其他求解常微分边值问题算法的简介5.2.1欧拉数值算法（显式）微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：(,)()dyf x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩(20) 其中a ，b 为常数.开始输入x 0,y 0,h,N x 1=x 0+hk 1=f(x 0,y 0),k 2=f(x 0+h/2,y 0+hk 1/2)k 3=f(x 0+h/2,y 0+hk 2/2),k 4=f(x 1,y 0+hk 3)y 1=y 0+h(k 1+2k 2+2k 3+k 4)/6n=1输出x 1,y 1n=N? 结束n=n+1x 0=x 1,y 0=y 1否是图5.1 龙格-库塔法流程图因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法.所有算法中的f 就是代表上式中(,)f x y ，而y f 表示(,)f x y y ∂∂，x f 表示(,)f x y x∂∂. 简单欧拉法是一种单步递推算法.简单欧拉法的公式如下所示：1(,)n n n n y y hf x y +=+ (21)简单欧拉法的算法过程介绍如下:给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+5.2.2欧拉数值算法（隐式）隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：111(,)n n n n y y hf x y +++=+ (22)隐式欧拉法是一阶精度的方法，比它精度高的公式是：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++ (23)隐式欧拉的算法过程介绍如下.给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+得出1k y +.5.2.3 欧拉预估-校正法改进欧拉法是一种二阶显式求解法，其计算公式如下所示：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++11(,)[(,)(,)]2n n n n n n n n t y hf x y h y y f x y f x t ++=+⎧⎪⎨=++⎪⎩ (24)四阶龙格-库塔法有多种形式，除了改进的欧拉法外还有中点法.中点法计算公式为：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++ (25)5.3 模型求解5.3.1运用MATLAB 软件对模型求解结果及分析用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题，并且比较其精度（具体的MATLAB 源程序见附录）以下进行实例分析：实例一. 1(0)1dyy x dx y ⎧=-++⎪⎨⎪=⎩由题可知精确解为：x y x e -=+，当x=0时，y(x)=0.在这里取步长h 为0.1，通过MATLAB 程序的计算，相应的结果如下：表5-1 步长为0.1时各方法的预测值与精确值的比较(精确到6位小数)初值 Euler 法相对误差预估校正法相对误差经典四阶库相对误差精确值 0 -- -- -- -- -- -- 1.00000 0.1 1.00910 0.00424 1.00500 0.00016 1.00484 0.00000 1.00484 0.2 1.02646 0.00759 1.01903 0.00029 1.01873 0.00000 1.01873 0.3 1.05134 0.01011 1.04122 0.00038 1.04082 0.00000 1.04082 0.4 1.08304 0.01189 1.07080 0.00045 1.07032 0.00000 1.07032 0.5 1.12095 0.01303 1.10708 0.00049 1.10653 0.00000 1.10653 0.6 1.16451 0.01366 1.14940 0.00052 1.14881 0.00000 1.14881 0.7 1.21319 0.01388 1.19721 0.00052 1.19659 0.00000 1.19659 0.8 1.26655 0.01378 1.24998 0.00052 1.24933 0.00000 1.24933 0.9 1.32414 0.01344 1.30723 0.00050 1.30657 0.00000 1.30657 1.01.38558 0.01294 1.36854 0.00048 1.36788 0.000001.36788步长为0.1时的精确值与预测值的比较精确值欧拉法改进欧拉格式四阶龙格库塔轴Y00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4X轴图5.2 步长为0.1时各方法的预测值与精确值的比较原函数图像轴YX轴图5.3 步长为0.1时原函数图像在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-2 h=0.05时三个方法与精确值的真值表步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0 -- -- -- -- -- --1.00000 0.05 1.00250 0.00911 1.00125 0.01035 1.00123 0.01037 1.00484 0.10 1.00738 0.01711 1.00488 0.01954 1.00484 0.01958 1.01873 0.15 1.01451 0.02405 1.01076 0.02765 1.01071 0.02770 1.04082 0.20 1.02378 0.03001 1.01880 0.03473 1.01873 0.03479 1.07032 0.25 1.03509 0.03507 1.02889 0.04085 1.02880 0.04093 1.10653 0.30 1.04834 0.03930 1.04092 0.04610 1.04082 0.04619 1.14881 0.35 1.06342 0.04277 1.05480 0.05053 1.05469 0.05063 1.19659 0.40 1.08025 0.04555 1.07044 0.05422 1.07032 0.05432 1.24933 0.45 1.09874 0.04772 1.08775 0.05724 1.08763 0.05734 1.30657 0.50 1.11880 0.04933 1.10666 0.05964 1.10653 0.05975 1.36788 0.55 1.14036 0.05045 1.12709 0.06150 1.12695 0.06161 1.00484 0.60 1.16334 0.05113 1.14895 0.06286 1.14881 0.06298 1.01873 0.65 1.18768 0.05143 1.17219 0.06379 1.17205 0.06391 1.04082 0.70 1.21329 0.05139 1.19674 0.06433 1.19659 0.06445 1.07032 0.75 1.24013 0.05106 1.22252 0.06453 1.22237 0.06465 1.10653 0.80 1.26812 0.05048 1.24949 0.06443 1.24933 0.06455 1.14881 0.85 1.29722 0.04969 1.27757 0.06408 1.27742 0.06419 1.19659 0.90 1.32735 0.04871 1.30673 0.06349 1.30657 0.06361 1.249330.95 1.35849 0.04759 1.33690 0.06272 1.33674 0.06283 1.306571.00 1.39056 0.04634 1.36804 0.06178 1.36788 0.06189 1.36788欧拉格式改进欧拉格式四阶龙格库塔精确值图5.4 步长为0.05时各方法的预测值与精确值的比较11.051.11.151.21.251.31.351.4X 轴Y 轴原函数图像图5.5 步长为0.05时原函数图像实例2 (0)1xdy xe ydx y -⎧=-⎪⎨⎪=⎩由题可知精确解为：21(2)2x x e -+当x=0时，y(x)=0.在这里取步长h为0.1，通过MATLAB程序的计算，相应的结果如下：步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.1 0.9000 0.0318 0.9096 0.0214 0.9094 0.0216 0.9295 0.2 0.8192 0.0611 0.8359 0.0420 0.8356 0.0424 0.8726 0.3 0.7544 0.0876 0.7761 0.0614 0.7757 0.0619 0.8268 0.4 0.7027 0.1109 0.7277 0.0793 0.7272 0.0799 0.7903 0.5 0.6617 0.1310 0.6886 0.0956 0.6881 0.0963 0.7615 0.6 0.6294 0.1480 0.6572 0.1103 0.6567 0.1110 0.7387 0.7 0.6040 0.1621 0.6320 0.1234 0.6315 0.1241 0.7209 0.8 0.5841 0.1737 0.6116 0.1349 0.6111 0.1355 0.70690.9 0.5686 0.1830 0.5951 0.1450 0.5946 0.1456 0.69591.0 0.5563 0.1905 0.5815 0.1538 0.5811 0.1544 0.6872原函数图像轴Y00.10.20.30.40.50.60.70.80.91X轴图5.6 步长为0.1时原函数图像各方法的预测值与精确值的比较欧拉格式改进的格式四阶龙格库塔精确值轴Y0.10.20.30.40.50.60.70.80.91X轴图5.7 步长为0.1时各方法的预测值与精确值的比较在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-4 步长为0.1时各方法的预测值与精确值的比较(精确到5位小数)步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.050.95000 0.01342 0.95245 0.01088 0.95242 0.01090 0.96292 0.100.90490 0.02650 0.90947 0.02158 0.90942 0.02163 0.92953 0.150.86428 0.03916 0.87067 0.03206 0.87061 0.03213 0.89951 0.200.82774 0.05137 0.83567 0.04228 0.83559 0.04237 0.87256 0.250.79491 0.06307 0.80414 0.05219 0.80405 0.05230 0.84842 0.300.76545 0.07423 0.77576 0.06176 0.77566 0.06188 0.82682 0.350.73904 0.08482 0.75023 0.07096 0.75013 0.07110 0.80754 0.40 0.71541 0.09482 0.72730 0.07977 0.72719 0.07991 0.79035 0.45 0.69427 0.10422 0.70672 0.08817 0.70660 0.08832 0.77505 0.50 0.67539 0.11302 0.68825 0.09615 0.68813 0.09630 0.76146 0.55 0.65855 0.12123 0.67168 0.10370 0.67156 0.10386 0.74940 0.60 0.64352 0.12886 0.65683 0.11084 0.65671 0.11100 0.73871 0.65 0.63012 0.13593 0.64351 0.11756 0.64340 0.11773 0.72925 0.70 0.61819 0.14245 0.63157 0.12388 0.63145 0.12405 0.72087 0.75 0.60754 0.14847 0.62085 0.12981 0.62073 0.12998 0.71347 0.80 0.59805 0.15400 0.61121 0.13537 0.61110 0.13553 0.70691 0.85 0.58957 0.15908 0.60254 0.14058 0.60243 0.14073 0.70109 0.90 0.58198 0.16374 0.59470 0.14545 0.59460 0.14560 0.695930.95 0.57517 0.16802 0.58761 0.15001 0.58751 0.15016 0.691321.00 0.56904 0.17194 0.58117 0.15429 0.58107 0.15443 0.687190.10.20.30.40.50.60.70.80.910.550.60.650.70.750.80.850.90.951X 轴Y 轴各方法下的预测值与精确值的比较欧拉格式改进欧拉格式四阶龙格库塔精确值图5.8 步长为0.05时各方法的预测值与精确值的比较六、模型的评价本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库塔法在求解初值问题的优越性.从上面的实例比较可知，在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小.这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对y的估计上，而精度则有所欠缺，以上各方法的选择，都取决于具体的情况.七、课程设计的总结与体会本文着重采用隆格库塔法运用MATLAB编程来求解微分方程，相比于欧拉法以及预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量不大、算法组织容易.其次，每一次的课程设计总是让我学到了更多的知识，不论是C++、SPASS还是MATLAB软件，这些都让我学到了如何解决实际问题的好工具，通过这些工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会：（1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从一点一点做起.以往训练得少只是把握的不牢靠.所以做起来感到有点吃力.所以，无论什么学科，一定要打好基础.（2）程序设计要靠多练，多见识，那样形成一种编程思维，我想对我是有很大好处的.尤其像我这种平时学得不扎实的人.（3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力，但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉.参考文献[1] 李庆杨，王能超，易大义编．数值分析（第四版）[M]：华中科技大学出版社，2006.[2] 姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：高等教育出版社，2005[3] 刘琼荪，数学实验[M]，高等教育出版社，2004[4] 王建伟，MATLAB7.X程序设计[M]，中国水利水电出版社，2007[5] 王高雄，周之铭等编.常微分方程（第三版） [M]:高等教育出版社，2006[6]何坚勇编著. 运筹学基础（第二版）[M]. 北京：清华大学出版社，2008.附录附录一：显示欧拉法matlab程序%欧拉法clear allclcx=[];y=[];y1=[];h=0.1;x=0:h:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold ony1(1)=1;for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);迭代n次的后退的欧拉格式matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%迭代n次的后退的欧拉格式clear allclcN=input('请输入迭代次数:');x=[];y=[];y1=[];h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); T=y1(j);for k=1:NT=y1(j-1)+h*f(x(j),T);endy1(j)=T;endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y); fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)附录二：预估校正法matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%预估校正法%欧拉法clear allclcx=[];y1(1)=1;y2=[];y2(1)=1;h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:ny1(j)=y2(j-1)+h*f(x(j-1),y2(j-1));y2(j)=y2(j-1)+(h/2).*[f(x(j-1),y2(j-1))+f(x(j),y1(j))]; endY=[x;y2];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y2,'r-');figure(2)DT=y-y2;plot(x,DT)附录三：四阶龙格库塔法matlab程序%四阶龙格库塔法clear allclcn=length(x);y=[];y1=[];y1(1)=1;for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:nK1=f(x(j-1),y1(j-1));K2=f(x(j-1)+h/2,y1(j-1)+h/2*K1);K3=f(x(j-1)+h/2,y1(j-1)+h/2*K2);K4=f(x(j-1)+h,y1(j-1)+h*K3);y1(j)=y1(j-1)+(h/6)*(K1+2*K2+2*K3+K4); endY=[x;y1];fid=fopen('data1.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)。

常微分方程数值解法-欧拉法、改进欧拉法与四阶龙格库塔法常微分方程数值解法

y( xn1)
y( xn
Байду номын сангаас
h)
y(xn )
hy'( xn )
h2 2!
y''( )
进一步: 令
h2 y( xn ) hy'( xn ) 2! y''( xn )
常微分方 yn1 y( xn1 ) , yn y( xn )
程数值解
法-欧拉法 yn1 yn hf ( xn , yn ) h2
、改进欧 y( xn1 ) yn1
2
max y''( x)
a xb
拉法和四
三、Euler方法
已知初值问题的一般形式为:
dy
dx
f (x, y)
a xb
(1)
y( x0 ) y0
常微分方用差商近似导数程数值解问题转化为
yn1 yn dy
h
dx
法-欧拉法 yn1 yn hf ( xn , yn )
法-欧 y(拉0) 法1
、改进欧
拉法和四
四、几何意义
由 x0 , y0 出发取解曲线 y yx 的切线（存在！），则斜率
dy
f x0, y0
dx x y
,
0
0
常微分方由于 f x0, y0 及 x0, y0 已知，必有切线方程。
由点斜式写出切程线方数程：值解
法、-改欧进拉欧法 ddxy y y0 x x0
常微分方程数值解能用解析方法求出精确解的微分方程为数不多,
而且有的方程即使有解析解,也可能由于解的表达
法-欧拉法式非常复杂而不易计算,因此有必要研究微分方程

matlab经典的4级4阶runge kutta法

MATLAB是一种用于算法开发、数据分析、可视化和数值计算的高级技术计算语言和交互式环境。

作为一个强大的工具，MATLAB提供了许多数值计算方法，其中4级4阶Runge-Kutta方法就是其中之一。

1. Runge-Kutta方法简介Runge-Kutta方法是求解常微分方程（ODE）的数值方法之一。

在MATLAB中，用户可以使用内置的ode45函数来调用4级4阶Runge-Kutta方法。

具体来说，4级4阶Runge-Kutta方法是一种单步迭代方法，通过在每个步骤中计算斜率来逐步逼近解析解。

它的优点是数值稳定性好，适用于多种类型的微分方程。

2. Runge-Kutta方法的公式4级4阶Runge-Kutta方法的一般形式如下：$$k_1 = hf(t_n, y_n)$$$$k_2 = hf(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_1)$$$$k_3 = hf(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_2)$$$$k_4 = hf(t_n + h, y_n + k_3)$$$$y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$其中，$t_n$是当前的独立变量值，$y_n$是当前的解向量，h是步长，$f(t_n, y_n)$是给定点$(t_n, y_n)$处的斜率。

通过不断迭代上述公式，可以逐步求解微分方程的数值解。

3. MATLAB中的4级4阶Runge-Kutta方法的应用在MATLAB中，用户可以使用ode45函数调用4级4阶Runge-Kutta方法来求解常微分方程。

使用ode45函数的基本语法如下：```matlab[t, y] = ode45(odefun, tspan, y0)```其中，odefun是用户定义的ODE函数句柄，tspan指定了求解的时间范围，y0是初始条件。

常微分方程初值问题数值解的实现和分析—四阶Rungekutta方法与预估校正算法毕业论文

《数值分析》课程设计常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法摘要求解常微分方程的初值问题，Euler方法，改进的Euler方法及梯形方法精度比较低，所以本文构造高精度单步的四级Runge-kutta方法及高精度的多步预估—校正算法及其Matlab编程来实现对常微分方程初值问题的求解，使在求解常微分方程时，对以前积分方法的收敛速度及精度都有了很高的提高。

关键词：Runge-kutta方法，Adams方法，预估—校正算法，Matlab目录1.前言 (1)2. 几个简单的数值积分法 (2)2.1Runge-kutta方法 (2)2.1.1 Runge-kutta方法的应用 (5)2.2预估—校正算法 (7)2.2.1 Adams数值积分方法简介及预估—校正算法 (7)2.2.2 预估—校正算法的应用 (12)3. 结果分析 (16)总结 (17)参考文献 (18)英文原文和中文翻译 (19)1英文原文 (19)2中文翻译 (20)1.前言常微分方程的初值问题是微分方程定解问题的一个典型代表，以下面的例子介绍常微分方程初值问题数值解的基本思想和原理。

例1.1 一重量垂直作用于弹簧所引起的震荡，当运动阻力与速度的平方成正比时，可借助如下二阶常微分方程描述若令和，则上述二阶常微分方程可化成等价的一阶常微分方程组类似于例1.1，对于m阶常微分方程其中。

若定义可得如下等价的一阶常微分方程组我们知道多数常微分方程主要靠数值解法。

所谓数值解法，就是寻求解在一系列离散节点上的近似值。

相邻两个节点之间的间距称为步长[1]。

2. 几个简单的数值积分法2.1 Runge-kutta方法Runge在1985年提出了一种基于Euler折线法的新的数值方法，此后这种新的数值方法又经过其同胞K.Heun和Kutta的努力[2]，发展完善成为后世所称的Runge-kutta 方法。

四阶龙格库塔法(Runge-Kutta)求解微分方程

四阶龙格库塔法（Runge-Kutta ）求解微分方程张晓颖（天津大学材料学院学号：1012208027）1 引言计算传热学中通常需要求解常微分方程。

这类问题的简单形式如下：{),(')(00y x f y y x y == （1）虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中的多数微分方程需要采用数值解法求解。

初值问题（1）的数值解法有个基本特点，它们采取“步进式”，即求解过程顺着节点排序一步一步向前推进。

这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。

2 龙格库塔法（Runge-Kutta ）介绍假设对于初值问题（1）有解 y = y (x ) ，用 Taylor 展开有：......)(!3)(!2)()()(321+'''+''+'+=+n n n n n x y h x y h x y h x y x y （2）龙格库塔法（Runge-Kutta ）实质上是间接的使用 Taylor 级数法的一种方法。

对于差商hx y x y n n )()(1-+，根据微分中值定理，存在 0 < θ < 1 ，使得：)()()(1h x y hx y x y n n θ+'=-+ （3）于是对于 y = y (x ) ，可得到：))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ （4）设))(,(*h x y h x f K n n θθ++=，为区间 [x n , x n +1 ] 上的平均斜率。

四阶龙格库塔格式中的*K 由下式计算得到：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n nn n n n n （5）四阶龙格库塔法（Runge-Kutta ）的每一步需要四次计算函数值f ，其截断误差更低，计算的精度更高。

Runge-Kutta 四阶

6 1.200000 1.778615
7 1.400000 1.516896
8 1.600000 1.245982
9 1.800000 1.062639
10 2.000000 0.999959
请选择函数
1.Y'(x)=y*sin(PIE*x) Y(0)=1
2.Y'(x)=x+y Y(0)=1.442
3.结束
1
请输入区间[a,b]的端点值
a=0
b=0
输入区间[a,b]端点值有误，请重新输入
a=0
b=2
请输入区间剖分点数n
n=10
n Xn Yn
1 0.200000 1.062680
2 0.400000 1.246016
3 0.600000 1.516916
4 0.800000 1.778619
5 1.000000 1.890103
goto main;
}
else;
}
printf("*****谢谢使用*****\n");
}
请ห้องสมุดไป่ตู้择函数
1.Y'(x)=y*sin(PIE*x) Y(0)=1
2.Y'(x)=x+y Y(0)=1.442
3.结束
4
选择错误，请重新选择
请选择函数
1.Y'(x)=y*sin(PIE*x) Y(0)=1
2.Y'(x)=x+y Y(0)=1.442
2.Y'(x)=x+y Y(0)=1.442
3.结束
3
*****谢谢使用*****
请按任意键继续. . .
for(i=0;i<MAX;i++){

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

若定义可得如下等价的一阶常微分方程组我们知道多数常微分方程主要靠数值解法。

所谓数值解法，就是寻求解在一系列离散节点上的近似值。

相邻两个节点之间的间距称为步长[1]。

Runge重要发现的灵感主要来源于把Euler方法应用于初值问题(2.1.1) 和f不显含y时的初值问题(2.1.2) 之间的类比。

Runge观察到，当把Euler方法应用到(2.1.2)型初值问题时，得到差分方程事实上，这是在计算积分问题(2.1.3) 时采用了左矩形法则即用高为，宽为h的矩形代替了(2.1.3)式中的积分值。

类似的，当把Euler方法应用到初值问题(2.1.1)时，得到差分方程这是在计算积分问题(2.1.4) 时采用了左矩形法则显然(2.1.3)和(2.1.4)式右侧数值积分的精度决定着数值解的精度。

此时，Runge发现，这个矩形公式的逼近程度并不高。

如果在计算积分问题(2.1.4)时采用中点法则或梯形法则，则数值解的精度会更高。

采用中点法则和梯形法则分别代替(2.1.4)式中的积分会得到(2.1.5)(2.1.6) (2.1.5)和(2.1.6)是两个隐式的方程，它们在早期发现Runge-kutta方法过程中扮演了一个重要角色。

Runge的出发点是：分别用一个Euler步代替(2.1.5)中未知的和(2.1.6)式右侧未知的值，这样Runge得到如下的中点格式：和梯形格式：这两个方程具有很好的几何解释。

它们是有折线构成的，这些折线假定微分方程所确定的斜率在前面的点上已经计算出来了。

与他的后继者一样，Runge用Taylor展开和系数对比说明上述两个方法的阶为2. Runge意识到，用中点法和梯形法得出的两个方法的阶数还不够高，所以，他设想，如果使用比中点法和梯形法精度更高的Simpson法则得到的方法阶数会提高，如果用M和T分别表示用中点法和梯形法算得的数值积分。

Simpson法则可以写成众所周知的形式S=M+(T-M)/3。

Taylor展开表明，如果f依赖于y，则这个表达形式只是2阶的，接着Runge发现，通过重复使用Euler步骤对梯形法则做适当的调整，会使格式=M+(-M)/3成为3阶方法，他还把他的方法及其Taylor 展开式拓展到微分方程组[3]。

五年后，Heun在其1900年的文章中评论Runge的方法时说，Runge获得的上述方法是归纳性的而且是令人费解的，他主使用更具一般性的Gauss求积公式其中于是可以把一般的Gauss格式扩充为(2.1.7)把(2.1.7)式的右端进行二元Taylor展开并与的Taylor展开式的对应的系数比较，适当选取参数使方法具有尽可能高的精度。

上述算法公式中的系数希望如此确定，使得(2.1.7)的Taylor展开式中所有h幂次不超过p的那些项与中的相应项的系数相等[2]。

假定p=4并取s=3，用类似的推导，可建立下述常用的显型四级Runge-kutta方法：截断误差为，当右端函数f不依赖于y时，上述公式可简化为2.1.1 Runge-kutta方法的应用用四级Runge-kutta方法求初值问题=y2 e-x，y(1)=1在区间[1,2]上的数值解(取h=0.1) (1) 先求出常微分方程的精确解，再用四阶Runge-kutta方法进行求解。

Y1=dsolve(‘Dy=y^2*exp(-x)’,’y(1)=1’,’x’)运行后的精确解为：Y1=1/(1/exp(x) - 1/exp(1) + 1)现给出常用的四阶Runge-kutta方法求解常微分方程的Matlab主程序如下：function [X,Y]=Rungek(funfcn,x0,b,y0,h)x=x0;y=y0;n=fix((b-x0)/h); %求步数i=1;X=zeros(n,1);Y=zeros(n,1);X(i)=x0;Y(i)=y0; %赋初值for i=2:nk1=feval(funfcn,x,y);k2=feval(funfcn,x+h/2,y+h*k1/2);k3=feval(funfcn,x+h/2,y+h*k2/2);k4=feval(funfcn,x+h,y+h*k3);y=y+h/6*(k1+2*k2+2*k3+k4);Y(i)=y;x=x+h;X(i)=x;end %按照龙格-库塔方法进行求解X,Y1=1./(1./exp(X) – 1./exp(1) + 1), %精确解plot(X,Y,'mh',X,Y1,'bo') %绘图gridlegend('用四阶龙哥库塔方法计算的数值解','精确解y(x)'), %图形说明wcha=abs(Y-Y1), %截断误差自定义函数：function z=funfcn(x,y)z=y.^2.*exp(-x);在matlab工作窗口输入如下程序：x0=1,b=2;y0=1;h=0.1;[X,Y]=Rungek(funfcn,x0,b,y0,h)运行后得：X =1.00001.10001.20001.30001.40001.50001.60001.70001.80001.9000Y1 =1.00001.03631.07141.10541.1380 1.1692 1.1990 1.2273 1.2540 1.2793 wcha = 1.0e-007 * 0 0.0111 0.0290 0.0518 0.0777 0.1054 0.1338 0.1620 0.18960.215911.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9211.051.11.151.21.251.31.35用四阶龙哥库塔方法计算的数值解精确解y=(x)2.2 预估—校正算法2.2.1 Adams 数值积分方法简介及预估—校正算法构造高阶精度格式的一个重要途径是，将方程改写成如下的积分形式(2.2.1.1) 然后以适当的差值多项式近似替代上式中的被积函数，这时后面将要介绍的Adams差值方法的出发点和思想。

选取k+1个插值结点：。

假定y(x)在上述结点上的近似值已被求出，记为，并令，它被视为在x=上的近似值。

利用上述结点和近似值作k次插值多项式用它近似(2.2.1.1)中的被积函数，得到近似公式(2.2.1.2) 将由此公式确定的当做的近似值。

在上述插值中，被插值点x位于包含所有插值结点的最小区间[]的外部，故称(2)为外插公式。

容易看出，当k=0时，(2.2.1.2)就是向前Euler公式。

为了得到(2.2.1.2)的具体计算公式，令x=，并将在区间写成Newton向后插值公式的形式其中代表j阶向前差分算子。

引进记号则有(2.2.1.3)将表达式(2.2.1.3)代入(2.2.1.2)得到(2.2.1.4) 其中系数(2.2.1.5) 它们与k和n无关。

例如，由于所以外插方法(2.2.1.4)和(2.2.1.5)的截断误差为(2.2.1.6)已知差分可用函数值的线性组合表示，故又可将外插公式(2.2.1.4)改写成(2.2.1.7) 其中系数已造表，应用时可查。

Adams外插公式表K 公式0 ,123仍然从积分关系式(2.2.1.1)出发，将其中的被积函数改用以为插值结点的k次Lagrange插值多项式近似。

这里，由于被插值点x位于包含所有插值结点的最小区间[]的部，故称相应的积分公式为插公式[1]。

用类似的方法可得到Adams插公式表如下Adams插公式表K 公式123在求解常微分方程初值问题的数值积分法中，还有一类高阶精度方法是多步方法。

一般来说，隐式多步方法（如Adams插法）比相应的显式多步方法精度高，并有较好的稳定性质。

然而，隐式方法计算的每一步需要求解非线性方程，所以它的计算比显式多步方法要复杂得多。

为了克服隐式多步方法计算方面的困难，人们经常隐式多步法和显式多步法结合起来使用，即预估—校正的算法。

下面介绍预估和校正的具体做法和计算公式：隐式k步方法可写成如下形式：(2.2.1.8) 其中，可认为是已知的，只是关于的非线性方程。

通常采用迭代法求解，计算公式如下：(2.2.1.9)容易证明，当h适当小并且恰当地选取初值时，上述迭代过程是收敛的。

显而易见，隐式方法(2.2.1.8)的每一步的计算量由(2.2.1.9)的迭代次数决定，因此选取好的初始值是十分重要的。

一种非常自然的做法就是取为用相应的k 步外插法计算的值，即令(2.2.1.10) 由(2.2.1.10)和(2.2.1.9)构成的算法被称为预估—校正算法，其中(2.2.1.10)为预估公式(记为P)，(2.2.1.9)为校正公式(记为C)。

由于这里的预估值比较精确，因此校正中所需的迭代次数一般不会很多。

如果出现迭代收敛慢的情况，应缩小步长在进行计算。

预估—校正算法的计算步骤：首先用预估公式(2.2.1.10)得到，然后计算，紧接着使用校正公式(2.2.1.9)求出，这样完成了一步校正。