CFD2020-第10讲-常微分方程数值解法
t x
y
z
x
y
z
FDM
空间离散
U Q(U)
FVM
时间推进 t
DG
推荐方法： Runge-Kutta法
2阶
U (1) U n tQ(U n )
U n1 1/ 2U n 1/ 2U (1) 1/ 2tQ(U (1) )
3 阶 （TVD型）
U (1) U n tQ(U n ) U (2) 3/ 4U n 1/ 4U (1) 1/ 4tQ(U ) (1) U n1 1/ 3U n 2 / 3U (2) 2 / 3tQ(U ) (2)
仅迭代一次(n=1),将
y(0) m 1
写入第二式
ym1
ym
h 2
f (xm , ym ) f (xm1, ym hf (xm , ym ))
改进的Euler方法
Copyright by Li Xinliang
7
P阶Taylor级数法
设初值问题的解具有p+1阶连续导数，利用Taylor公式
Euler公式
K1 f (xm , ym )
ym1 ym hf (xm , ym )
N 2
ym1 ym h b1K1 b2K2
K1 f (xm , ym )
K2 f (xm c2h, ym ha21K1)
• 将 K2 在(xm , ym ) 处展开
K2
f
h(c2
fx
a21K1 f y )
y'''(xm ) O(h4 )
y' f
利 用
y'' f x ff y
y''' fxx 2 ffxy f 2 f yy ( fx ffy ) f y
二阶：O(h3 )
b1 b2 1 b2c2 1/ 2
c2 a21
无穷多组解
b2a21 1/ 2
Copyright by Li Xinliang
ym1 ym h(xm , ym , h)
d dx
f
(x,
y(x))
fx
f fy
d2 dx2
f (x, y(x))
f xx 2 ff xy
f 2 f yy
fy( fx
ff y )
……
(*)式导数的计算比较复杂 离散情形（非解析形式）
（很少使用）
Copyright by Li Xinliang
向前Euler方法
ym1 ym hf (xm1, ym1)
向后Euler方法
ym1
ym
h 2
f (xm , ym ) f (xm1, ym1)
梯形法
迭求代解
y(0) m1
y(n) m1
ym ym
hf (xm , ym )
h 2
f (xm , ym )
f (xm1,
y ) (n1) m1
K1 f (xm , ym )
K2 f (xm c2h, ym ha21K1)
f
h(c2
fx
a21K1 f y )
h2 2
(c22
f xx
2c2a21K1 f xy
a221K12
fyy)
h3 6
(c23 f xxx 3c22a21K1 f xxy 3c2a221K12
f
xyy
y(xm1)
p j0
hj j!
y
(
j)
(
xm
)
(
h p
p1
1)!
y(
p1)
(m
)
y' f (x, y)
(x, y(x),h)
p h j1 j1 j!
d j1 dx j1
f (x, y(x))
(*)
y ( xm 1 )
y ( xm
)
h(xm ,
y(xm ),
h)
(
h p1 p 1)!
y( p1) (m )
f yyy)
m
O(h4 )
ym1 ym b1hK1 b2hK2 b3hK3 O(h4 )
ym b1hf b3h[ f
b2h[ h(c3 fx
f
h(c2 fx (a31K1 a32
a21K1 K2) fy
f )
y)
h2
2
h2 2
(c22
(c32 f xx
f xx 2c2a21K1 f xy 2c3 (a31K1 a32K2
• ym1 ym h(b1 b2 ) f h2b2 (c2 fx a21 ffy )
h3 2
b2 (c22
f xx
2c2a21
ff
xy
a221
f
2
f yy)
O(h4 )
• 将 y(xm1) 在 xm 处展开
y ( xm1 )
y(xm ) hy'(xm )
h2 2
y''(xm )
h3 3!
Rm1 y(xm1) y(xm ) h(xm , y(xm ), h)
从准确值y( xm )计算1步后得 到的值与准确值y(xm1) 之差
整体截断误差:
m1 y(xm1) ym1
P阶方法:
从初值出发，计算m+1步后 得到的值与准确值之差
局部截断误差 O(h p1) 显式方法:
右端项不包含 ym1 隐式方法:
xm1 xm
f (x, y(x))dx
• 利用Taylor公式
y ( xm1 )
y(xm ) hy'(xm )
h2 2
y''(m )
Copyright by Li Xinliang
5
单步法
计算ym1 的值只用到前一步的值 ym
ym1 ym h(xm , ym , h)
增量函数
局部截断误差:
h3 2
b2 (c22
f xx
2c2a21
ff
xy
a221
f
2
f yy)
O(h4 )
Copyright by Li Xinliang
10
Runge-Kutta方法
N阶
N
ym1 ym h bi Ki i 1
N 2 ym1 ym h b1K1 b2K2
K1 f (xm , ym ) K2 f (xm c2h, ym ha21K1)
4阶
U (1) U n 1/ 2tQ(U n ) U (2) U n 1/ 2tQ(U ) (1) U (3) U n tQ(U (2) ) U n1 1/ 3(U n U (1) 2U (2) U ) (3) 1/ 6tQ(U ) (3)
更高阶 ……
Copyright by Li Xinliang
右端项含有 ym1
• 一阶向后差商代替一阶导数
y(xm1) h
y(xm )
y' ( xm1 )
• 两边从xm 到xm1 积分,积分采用 右矩形公式
y(xm1) y(xm )
xm1 xm
f (x, y(x))dx
Copyright by Li Xinliang
6
改进的Euler方法
ym1 ym hf (xm , ym )
y
y' f (x, y)
y2
y1
y(x0 ) y0
初值条件
y0
因为f(x,y(x))恰好是待求精确解y(x)的斜率y’(x)
y(x)是在x-y平面上经过点(x0,y0)的曲线,该曲线
在点(x0,y0)处的切线的斜率为f(x0,y0).
y y(x)
x0
x1 x2
x
Euler方法
y y0 f (x0 , y0 )( x x0 )
x x1 x h
y1 y0 hf (x0 , y0 ) y2 y1 hf (x1, y1)
ym1 ym hf (xm , ym )
• 一阶向前差商代替一阶导数
y(xm1) h
y(xm )
y'(xm )
• 两边从 xm到xm1 积分,积分采用
左矩形公式
y(xm1) y(xm )
11
Runge-Kutta方法
N阶
N 2
ym1 ym h b1K1 b2K2
K1 f (xm , ym ) K2 f (xm c2h, ym ha21K1)
2阶格式满足的条件：
N
ym1 ym h bi Ki i 1
b1 b2 1 b2c2 1/ 2 b2a21 1/ 2
y Ce f (x)dx
y(x) c'(x) g(x)e f (x)dx
• 常微分方程的数值计算
P(x, y) Q(x, y) y' 0
U P, x
U Q y
d U (x, y(x)) 0 dx
U (x, y(x)) C
计算机的发展、绝大多数常微分方程难以求解析解
U f1(U) f2 (U) f3 (U) g1(U) g2 (U) g3 (U)
f (x, y) 1
(y x)
i 1
ci aij , i 2,, N j 1
新的ym1与原来的 x 在计 算过程中有相同的增量
N
bi 1
i 1
Copyright by Li Xinliang
13
N
Runge-Kutta方法
ym1 ym h bi Ki
i 1
N 3 将 K2 K3 在 (xm, ym ) 处展开
K1 f (xm , ym )
i 1
Ki f (xm cih, ym h aij K j ) j 1
第10章_常微分方程数值解
迭代法的收敛条件：
计
算 方
|
y y (k 1) i 1
(k ) i 1
|
h| 2
f
( xi1,
y(k) i 1
)
f
( xi1,
y (k 1) i 1
)
|
法 课
h 2
L|
y(k) i 1
y (k 1) i 1
|
件
其中L为Lipschitz常数。
当
0 h L 1
2
时迭代过程收敛。即迭代过程中误差没有扩大。
称为欧拉中点公式.容易看出,中点公式计算yi+1时,不仅需要yi的 值，还需要yi-1的值.
8
结束
凡是从已知(或已算出)的y0,y1,…,yi能直接从公式算出yi+1的公 式称显式，否则称隐式公式(8.4),(8.6)都是显式，(8.5)是隐式
在计算yi+1时，只需要yi的值，则称公式为单步法；若除yi之外 还需要以前的yi-1等多个值，则称多步法公式.(8.4)、(8.5)是单 计 步法，(8.6)是多步法，确切说是二步法.
0.2 xi yi
计 算
yi
1
yi
0.05
yi
2xi yi
y(0) i 1
2 xi 1 y (0)
i 1
方
y0 1
法 课
件 结果（略）
15
结束
§8.3 龙格-库塔(Runge-Kutta)法
欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出
高精度的显式一步法：龙格-库塔法,简称R-K法. 计 算 8.3.1 二阶R-K法
8.2.1 欧拉方法的导出
把区间［a,b］分为n个小区间，取步长h=(b-a)/ n ，节点 xi=x0+ih,i=0,1,2,…，n，其中x0=a，又设y (x)为上述问题的解.
常微分方程数值解法
用分段的折线逼近函数,此为 “折线法”而非“切线法”， 除第一个点是曲线上的切线，
其它都不是。
2、Euler方法的误差估计
1)局部截断误差。 在一步中产生的误差而非累积误差：
~
T x y y
n1
n1
n1
其中
~
y
是当
y
n
y(
x
)
n
（精确解！）时
n1
由Euler法求出的值，即y 无误差！ n
y y x , y x y 则得： h f
f ,
n1
n2
n
n
n1
n1
同样与Euler法结合，形成迭代算法，对n 0，1，2，
y y x , y 0 hf
n1
n
n
n
y y x , y x y
k 1
推出总体误差与步长的关系。
由微分方程解的存在唯一性,自然假定 （ f x，y）
充分光滑，或满足 Lipschitz条件：
f
x,ny源自xn
f
x
,
n
y
n
L
yxn
y n
第 n 步 的 总 体 截 断 误 差 记 为
en y
xn
y n
则 对 n 1 步：
e x y x y y y T y y ~ ~
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx
hf
(xn ,
yn )
常微分方程的数值解法
常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时na b h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1．对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
常微分方程的数值解法
常微分方程的数值解法…………江南大学信计1203柯恒一、前言对于很多微分方程，我们很难求出解析解，这时我们需要采取数值手段求解。
在数值分析这门课中，老师讲到dy dt =f (t ,y )的数值解法，我们采用了欧拉格式，后退欧拉格式，梯形公式，改进欧拉格式及4阶龙格-库塔方法求解。
而老师没讲关于微分方程组和高阶微分方程的解法。
现在我来粗虐的说一下微分方程组及高阶微分方程的数值解，这里主要是借助matlab 中的相关函数进行计算并仿真出图像。
二、相关问题初值问题问题1，微分方程组问题描述：给定一个微分方程，并且告诉我们初始状态。
我们便可求出整个过程的解。
给定下面方程（L 表示省略号）：11122112112(,,,,)(,,,,)(,,,,)n n n n dy f x y y y dx dy f x y y y dx dy f x y y y dx ⎧=⎪⎪⎪=⎪⎨⎪⎪⎪=⎪⎩而且初始状态y ′i (x 0)已知记成y(0).问题解决：四阶龙格库塔方法：K i1=f i (x n ,y 1n ,y 2n ,……,y nn );K i2=f i (x n +h 2,y 1n +h 2∗k 11,y 2n +h 2∗K 21,……,y nn +h 2∗K n1) K i3=f i (x n +h ,y 1n +h ∗K 12,y 2n +h ∗K 22,……,y nn +h ∗K n2) K i4=f i (x n +h 2,y 1n +h 2∗K 13y 2n +h 2∗K 32,……,y nn +h 2∗K n3),y i,n+1=y i,n+h6∗(K i1+2∗K i2+2∗k i3+2∗K i4)通过龙格库塔方法，我们可以计算后面点的值，在matlab中调用ode45来实现四阶龙格库塔方法的调用。
应用举例：有一个同步地球卫星，现在加速到4km/s进行变轨，试分析变轨后卫星运动的轨迹。
已知地球的质量M=5.97e24,引力常量的值为6.672e-11问题建模：我们已经知道行星的运动是在一个平面上的。
数值分析常微分方程数值解
许多实际问题的数学模型是微分方程或微分方程的定解问题。
如物体运动、电路振荡、化学反映及生物群体的变化等。
常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。
若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。
因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。
常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。
用解析方法只能求出线性常系数等特殊类型的方程的解。
对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。
因此研究微分方程的数值解法是非常必要的。
只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。
由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。
在下面的讨论中，我们总假定方程满足以上两个条件。
所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。
),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由到1+n x 的步长。
今后如无特别说明，我们总假定步长为常量。
建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
常微分方程数值解法5262115页PPT文档
r 表示食饵独立生存时的增长率；
d 表示捕食者独立生存时的死亡率；
a 表示捕食者的存在对食饵增长的影响系数，反映捕
食者对食饵的捕获能力；
b 表示食饵的存在对捕食者增长的促进系数，反映食
饵对捕食者的喂养能力
150 100
令 y 1 y ,y 2 y ',y 3 y '', ,y n y ( n 1 )
可以将以上高阶微分方程化为如下一阶常微分方程组
y1 ' y2 y2 ' y3 yn ' an(x)y1
a1(x)yn f (x)
例：P120,1(a),Bessel方程
常微分方程的数值解
一般地，凡表示未知函数，未知函数的导 数与自变量之间的关系的方程叫做微分方 程．未知函数是一元函数的，叫常微分方 程;未知函数是多元函数的，叫做偏微分方 程．
如
y ' x y'x2y2 y''y'xy
Matlab实现 [t,x]=ode45(f,ts,x0,options,p1,p2,......)
50 0 0
30 20 10
0 0
10
20
50
30
20
10
0
30
0
10
8
6
4
2
100
0
50
100
150
50
100
高阶常微分方程的解法
高阶常微分方程
y ( n ) a 1 ( x ) y ( n 1 ) a ( n 1 ) ( x ) y ' a n ( x ) y f( x )
常微分方程的数值解法
第六章 常微分方程的数值解法 §6.0 引言§6.1 算法构造的主要途径§6.2 Runge-Kutta Method 算法§6.3 线性多步法§6.4 线性多步法的一般形式§6.5 算法的稳定性、收敛性§6.0 引 言1． 主要考虑如下的一阶常微分方程初值问题的求解:()()00,dy f x y dx y x y ⎧=⎪⎪⎨=⎪⎪⎩ 微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。
2. 例如微分方程:xy '-2y=4x ；初始条件: y(1)=-3。
于是可得一阶常微分方程的初始问题24(1)3y y x y ⎧'=+⎪⎨⎪=-⎩。
显然函数y(x)=x 2-4x 满足以上条件,因而是该初始问题的微分方程的解。
3. 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。
因此，只能依赖于数值方法去获得微分方程的数值解。
4． 微分方程的数值解：设微分方程问题的解y(x)的存在区间是[a,b ]，初始点x 0=a ，将[a,b ]进行划分得一系列节点x 0 , x 1 ,...,x n ，其中a= x 0< x 1<…< x n =b 。
y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点x k 的近似值y(x k )，即y≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
如图所示:§6.1 算法构造的主要途径x 0 x 1 x 2 ...1 欧拉公式1.1 构造的思想：利用差商代替一阶导数，即010()()x x y x y x dy dx h =-≈，则 1000()()(,)y x y x f x y h -≈。
常微分方程数值解法
第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。
另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。
因此，有必要探讨常微分方程初值问题的数值解法。
本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。
显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。
定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。
一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。
以上方法称 为欧 拉法或欧拉折线法。
常微分方程数值解法
i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。
我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。
通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。
本章讨论常微分方程最常用的近似数值解法--差分方法。
先来讨论最简单的Euler 法。
为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。
在微积分课程中我们熟知，微商（即导数）是差商的极限。
反过来，差商就是微商的近似。
在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
下面我们用数值积分法重新导出 Euler 法以及其它几种方法。
为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。