9299.net
大学生考试网 让学习变简单
当前位置:首页 >> >>

常微分方程的数值解.ppt_图文

常微分方程的数值解.ppt_图文

常微分方程的数值解法
电子科技大学

常微分方程的数值解
? ? ?

?

引言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程

6.1 引言

在高等数学中我们见过以下常微分方程:
? y ? ? f ( x , y ) a ? x ? b ( 2)? y ?? ? f ( x , y , y ?) a ? x ? b ? (1)? ? y(a ) ? y 0 , y ?(a ) ? ? ? y( a ) ? y 0 ? y ?? ? f ( x , y , y ?) a ? x ? b ( 3)? ? y(a ) ? y 0 , y( b ) ? y n

(1),(2)式称为初值问题,(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组:
? ? f 1 ( x , y1 , y 2 ) ? y1 ( 4 )? 2 ? f 2 ( x , y1 , y 2 ) ? y? y1 ( x 0 ) ? y10 y 2 ( x 0 ) ? y 20

本章主要研究问题(1)的数值解法,对(2)~(4)只 作简单介绍。

考虑一阶常微分方程初值问题

? y? ? f ( x, y) (1)? ? y( x 0 ) ? y 0

其中,y = y(x) 是未知函数,y(x0) = y0 是初值 条件,而f (x, y) 是给定的二元函数. 由常微分方程理论知,若f(x)在x?[a,b]连续且 f 满足对 y 的Lipschitz条件:

f ( x, y1 ) ? f ( x, y2 ) ? L y1 ? y2
(其中L为Lipschitz常数)则初值问题(1)存 在唯一的连续解。

求问题(1)的数值解,就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。
为了计算方便,可取 xn=x0+nh,(n=0,1,2,…), h称为步长。
?

常微分方程的数值解法有单步法和多步法之分:
?

单步法:在计算yn+1 时只用到前一点yn 的值 ;

?

多步法:计算yn+1 时不仅利用yn,还要利用yn-1, yn-2,……. 一般k步法要用到 yn, yn-1, yn-2,……., yn-k+1。

6.2 简单的数值方法
一、欧拉(Euler)方法

在x= x0 处,用差商代替导数:

y( x1 ) ? y ( x 0 ) y( x 1 ) ? y( x 0 ) y ?( x 0 ) ? ? x1 ? x 0 h
由 y?( x0 ) ? f ( x0 , y0 ), 得

y( x0 ) ? y0

y( x1 ) ? y0 ? hf ( x0 , y0 ) ? y1

同理,在x= xn 处,用差商代替导数: y ( x n ? 1 ) ? y ( x n ) y( x n ? 1 ) ? y( x n ) y ?( x n ) ? ? x n ?1 ? x n h 由 y?( x n ) ? f ( x n , yn ) 得

y( x n?1 ) ? y( x n ) ? hf ( x n , yn ) y( x n ) ? y n , y( x n?1 ) ? yn?1

若记

则上式可记为 y n?1 ? y n ? hf ( x n , y n ) 此即为求解初值问题的 Euler 方法,又称显式 Euler方法。

Euler方法的几何意义:(Euler折线法)
y
P2 P1 P0 Pn Pn+1

y=f(x) x

O

x0 x1

x2

xn

xn+1

例: 用Euler方法求解常微分方程初值问题
y ? 2 ? y ? ? 2 y (0 ? x ? 3 ) ? x ? ? ? y(0) ? 0.

并将数值解和该问题的解析解比较。 x 解析解: y( x ) ? 1 ? x2 解:Euler方法的具体格式:

y n?1

yn 2 ? y n ? h( ? 2 yn ) xn

取h=0.2, xn=nh,(n=0,1,2…,15), f(x,y)=y/x – 2y2 计算中取f(0,0)=1. 计算结果如下:
xn 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y(xn) 0 0.1923 0.3448 0.4412 0.4878 0.5000 0.4918 0.4730 yn 0 0.2000 0.3840 0.5170 0.5824 0.5924 0.5705 0.5354 yn-y(xn) 0 0.0077 0.0392 0.0758 0.0946 0.0924 0.0787 0.0624

xn 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

y(xn) 0.4494 0.4245 0.4000 0.3767 0.3550 0.3351 0.3167 0.3000

yn 0.4972 0.4605 0.4268 0.3966 0.3698 0.3459 0.3246 0.3057

yn-y(xn) 0.0478 0.0359 0.0268 0.0199 0.0147 0.0108 0.0079 0.0057

由表中数据可以看到,微分方程初值问题的数值解 和解析解的误差一般在小数点后第二位或第三位小 数上,这说明Euler方法的精度是比较差的。

h=0.2;y(1)=0.2;x=h:h:3; for n=1:14 xn=x(n);yn=y(n); y(n+1)=yn+h*(yn/xn-2*yn*yn); end x0=0:h:3;y0=x0./(1+x0.^2); plot(x0,y0,x,y,x,y,'o')
0 .8 0 .6

0 .4

- : 数值解 — : 准确解
0 0 .5 1 1 .5 2 2 .5 3

0 .2

0

若直接对y`=f(x,y)在[xn, xn+1]积分,
y( x n ? 1 ) ? y( x n ) ? ?
xn ? 1 xn

f ( x, y( x ))dx

利用数值积分中的左矩形公式:

?

xn ? 1

xn

f ( x, y( x ))dx ? h ? f ( x n , y( x n ))

设y(xn)= yn,则得

yn?1 ? yn ? hf ( xn , yn )
此即为Euler公式。 若用右矩形公式:

?

xn ? 1

xn

f ( x, y( x ))dx ? h ? f ( x n?1 , y( x n?1 ))

得 y n?1 ? y n ? hf ( x n?1 , y n?1 )
上式称后退的Euler方法,又称隐式Euler方法。 可用迭代法求解:
( 0) y 初值: n?1 ? yn ? hf ( xn , yn )
( k ?1) (k ) 迭代: yn?1 ? yn ? hf ( xn?1 , yn?1 ) k=0,1,……

( k ?1) (k ) 因 yn ? y ? h f ( x , y ?1 n?1 n?1 n?1 ) ? f ( x n?1 , y n?1 ) (k ) ? hL y n ?1 ? y n?1

故当hL<1时,迭代法收敛。

二、梯形方法
由 y( x n?1 ) ? y( x n ) ? ?x
xn ? 1
n

f ( x, y( x ))dx

利用梯形求积公式: xn ? 1 h ?xn f ( x, y( x ))dx ? 2 ? f ( x n , y( x n )) ? f ( x n?1 , y( x n?1 ))? 得 y n?1
h ? y n ? ? f ( x n , y n ) ? f ( x n ? 1 , y n ? 1 )? 2

上式称梯形方法,是一种隐式方法。

用迭代法求解: ( 0) 初值: yn?1 ? yn ? hf ( xn , yn )

迭代: y

( k ?1) n?1

h (k ) ? yn ? f ( x n , y n ) ? f ( x n?1 , y n?1 ) 2

?

?

k=0,1,……

h (k ) 因 y ( k ?1) ? y f ( x n?1 , y n ?1 ) ? f ( x n?1 , y n?1 ) n?1 n?1 ? 2 hL ( k ) ? y n?1 ? y n?1 2
故当hL/2<1时,迭代法收敛。

由以上分析可以看出,隐式方法的计算比显 式方法复杂,需要用迭代法求解非线性方程 才能得出计算结果。 可采用将显式 Euler 格式与梯形格式结合使用 的方法来避免求解非线性方程。 记 y n?1 ? y n ? hf ( x n , yn ) 再用梯形格式计算: --预测

y n?1

h ? y n ? ? f ( x n , y n ) ? f ( x n ? 1 , y n ? 1 )? 2
--校正

上面两式统称预测-校正法, 又称改进的Euler方法。

三、单步法的局部截断误差和精度 单步法的一般形式为 : (?与f 有关)

yn?1 ? yn ? h? ( x n , yn , yn?1 , h)
显式单步法形式为:

yn?1 ? yn ? h? ( xn , yn , h)
整体截断误差:从x0开始,考虑每一步产生的 误差,直到xn,则有误差 en ? y( xn ) ? yn 称为数值方法在节点xn处的整体截断误差。 但en不易分析和计算,故只考虑从xn到xn+1的局部 情况。

定义:设y(x)是初值问题(1)的精确解,则称

Tn?1 ? y( xn?1 ) ? y( xn ) ? h? ( xn , y( xn ), h)
为显式单步法在节点xn+1处的局部截断误差 。 若存在最大整数p使局部截断误差满足

Tn?1 ? y( x ? h) ? y( x) ? h? ( x, y, h) ? O(h

p?1

)

则称显式单步法具有p阶精度或称p阶方法。 注:将 Tn+1 表达式各项在 xn 处作 Taylor 展开,可 得具体表达式。

Euler方法的局部截断误差:

(设yn=y(xn))

Tn?1 ? y( xn?1 ) ? y( xn ) ? hf ( xn , y( xn ))

? y( xn ? h) ? y( xn ) ? hy?( xn )
2

h ? y( x n ) ? hy ?( x n ) ? y ??(? n ) ? y( x n ) ? hy ?( x n ) 2
h2 h2 ? y ??(? n ) ? y ??( x n ) ? O( h 3 ) ? n ? ( xn , xn?1 ) 2 2

其中 (h2 / 2) y??( xn ) 称局部截断误差主项。 故Tn+1= O(h2),p=1, 即Euler方法具1阶精度。

梯形方法的局部截断误差: (设yn=y(xn))
Tn?1 h ? y( x n?1 ) ? y( x n ) ? [ y ?( x n ) ? y ?( x n?1 )] 2 h2 h3 ? hy ?( x n ) ? y ??( x n ) ? y ??( x n ) 2 3!
h h2 ? [ y ?( x n ) ? y ?( x n ) ? hy ??( x n ) ? y ???( x n )] ? O( h 4 ) 2 2

h ?? y ???( x n ) ? O( h 4 ) 12

3

故Tn+1= O(h3),p=2, 梯形方法具2阶精度。

局部截断误差主项为:? (h / 12) y???( xn )
3

6.3 Runge-Kutta方法
一、Runge-Kutta方法的基本思想 由Taylor展式
y( x n?1 ) ? y( x n ? h) ? y n ?1 h2 h p ( p) ? y( x n ) ? hy ?( x n ) ? y ??( x n ) ? ? ? y ( xn ) 2 p!

Tn+1= O(hp+1),若提高p,可提高精度。 但因 y ? ? f ( x, y ) …… 高阶导数计算复杂,故可从另外角度考虑。

? ( x, y ) ? f y ? ( x , y ) ? f ( x, y ) y??( x) ? f x

分析Euler公式及改进的Euler公式:
? y n?1 ? y n ? hK 1 ? ? K 1 ? f ( xn , yn )
K1 K 2 ? ? y n ?1 ? y n ? h( 2 ? 2 ) ? ? K 1 ? f ( xn , yn ) ? K ? f ( x ? h, y ? hK ) n n 1 ? 2 ?

局部截断误差:O(h2)

局部截断误差:O(h3)

可用f(x,y)在某些点处值的线性组合得 yn+1,增加 计算f(x,y)的次数可提高阶数。

Runge-Kutta方法的基本思想:

设法计算 f(x,y) 在某些点上的函数值,然后对这 些函数值作线性组合,构造近似计算公式,再 把近似公式和解的泰勒展开式相比较,使前面 的若干项吻合,从而获得达到一定精度的数值 计算公式 。 r 设 ?
? y n?1 ? y n ? h? c i K i i ?1 ? ? ? K 1 ? f ( xn , yn ) ? i ?1 ? K i ? f ( x n ? ? i h, y n ? h? ? ij K j ), ( i ? 2,3,? , r ) ? j ?1 ?

ci , λi , μij 为待定常数。

上面第一个式子的右端在(xn,yn)作泰勒展开后,按h 的幂次作升序排列 :

y n?1

1 1 2 3 ? yn ? ? 1h ? ? 2 h ? ? 3 h ? ? 2! 3!

再与初值问题的精确解 y(x) 在点 x=xn处的泰勒展 开式 y( xn?1 ) ? y( xn ? h) ? y( xn ) ? hy?( xn )
h2 h p ( p) ? y ??( x n ) ? ? ? y ( x n ) ? O( h p ? 1 ) 2! p!

相比较,使其有尽可能多的项重合。

例如,要求

? 1 ? f n , ? 2 ? f n? , ? 3 ? f n??,?, ? p ? f n( p?1)
就得到p个方程,从而定出参数 ci ,?i,?ij ,再 代入K1,K2,…, Kr的表达式,就可得到计算 微分方程初值问题的数值计算公式 :
y n?1 ? y n ? h? c i K i
i ?1 r

上式称为r 级Runge-Kutta方法的计算公式。

若 Tn+1= O(hp+1),则称其为p 阶r 级R-K方法。
当r=1时,就是Euler方法。

要使Runge-Kutta公式具有更高的阶 p,就要增 加r 的值。下面我们只就r =2推导R-K方法。
二、二阶Runge-Kutta方法
? y n?1 ? y n ? h(c1 K 1 ? c 2 K 2 ) ? ? K 1 ? f ( xn , yn ) ? K ? f ( x ? ? h, y ? ? hK ) n 2 n 21 1 ? 2
其中 c1, c2, ?2, ?21 待定。

上式的局部截断误差为: Tn?1 ? y( x n?1 ) ? y n ? h[c1 f ( x n , y n ) ? c 2 f ( x n ? ? 2 h, y n ? ? 21 hf n )]

2 h 3 由 y( x n?1 ) ? y n ? hy ? ? ? ? y ? O ( h ) n n 2!

利用二元函数的Taylor展开,得
? y? n ? f ( xn , yn ) ? f n ? ? d ? y? f ( x n , y( x n )) ? f x? ( x n , y n ) ? f y? ( x n , y n ) ? f n n ? ? dx ?
又 f ( x n ? ? 2 h, y n ? ? 21 hf n )

2 ? ? ? f n ? f x ( x n , y n )? 2 h ? f y ( x n , y n )? 21 hf n ? O( h )

代入Tn+1的表达式,得

Tn?1 ? y( x n?1 ) ? y n ? h[c1 f ( x n , y n ) ? c 2 f ( x n ? ? 2 h, y n ? ? 21 hf n )]
h2 3 ? ? ? y n ? hy ? ? y ? O ( h ) n n 2! ? y n ? h[c1 f ( x n , y n ) ? c 2 f ( x n ? ? 2 h, y n ? ? 21 hf n )]

h ? hf n ? [ f x? ( x n , y n ) ? f y? ( x n y n ) f n ] ? h[c1 f n 2 ? c 2 ( f n ? ? 2 f x? ( x n , y n )h ? ? 21 f y? ( x n , y n ) f n h)] ? O(h 3 ) 1 ? (1 ? c1 ? c 2 ) f n h ? ( ? c 2 ? 2 ) f x? ( x n , y n )h 2 2 1 ? ( ? c 2 ? 21 ) f y? ( x n , y n ) f n h 2 ? O( h 3 ) 2

2

要使上式p=2阶,则需
? ?1 ? c1 ? c 2 ? 0 ? ?1 ? ? c2?2 ? 0 ?2 ?1 ? c 2 ? 21 ? 0 ? ?2



? ?c1 ? c 2 ? 1 ? 1 ? ?c 2 ? 2 ? 2 ? 1 ? c 2 ? 21 ? ? 2 ?

方程组解不唯一,可令c2=a ? 0 ,则 c1 = 1-a , ?2 = ?21 =1/(2a) 满足上述条件的公式都为2阶R-K公式。

如取a= ? ,则c1= c2= ?, ?2=?21 =1,即为改 进Euler公式。 若取a= 1,则c1= 0,c2= 1, ?2=?21 = ? ,得
? ? y n?1 ? y n ? hK 2 ? ? K 1 ? f ( xn , yn ) ? h h ? K 2 ? f ( xn ? , yn ? K 1 ) 2 2 ?

称中点公式,相当于数值积分的中矩形公式:

y n?1

h h ? y n ? hf ( x n ? , y n ? f ( x n , y n )) 2 2

例:蛇形曲线的初值问题
y ? 2 ? ? y ? ? 2 y (0 ? x ? 3 ) x ? ? ? y(0) ? 0.

令f(x,y)=y/x –2y2, 取f(0,0)=1, h=0.2, xn=nh , ( n = 1,2,…,15) 2阶龙格-库塔公式计算格式: k1=yn/xn – 2yn2, k2 = (yn+hk1)/(xn+h) – 2(yn+hk1)2

yn+1=yn + 0.5h[ k1 + k2]

x0=0;y0=0;h=.2;
x=.2:h:3; k1=1; k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2; y(1)=y0+.5*h*(k1+k2); for n=1:14 k1=y(n)/x(n)-2*y(n)^2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2; y(n+1)=y(n)+0.5*h*(k1+k2); end y1=x./(1+x.^2); subplot(221) plot(x,y1,x,y1,'b*') subplot(222) plot(x,y,x,y,'o') subplot(223) plot(x,y,x,y,'o',x,y1,x,y1,'b*‘)

三、三阶与四阶Runge-Kutta方法 当r=3时,R-K公式表示为
? y n?1 ? y n ? h(c1 K 1 ? c 2 K 2 ? c 3 K 3 ) ? ? K 1 ? f ( xn , yn ) ? ? K 2 ? f ( x n ? ? 2 h, y n ? ? 21 hK 1 ) ? ? K 3 ? f ( x n ? ? 3 h, y n ? ? 31 hK 1 ? ? 32 hK 2 )

其中 c1 , c2 , c3 , ?2 , ? 21 , ?3 , ?31 , ?32 为8个待定常数。 上式的局部截断误差为
Tn?1 ? y( x n?1 ) ? yn ? h[c1 K 1 ? c2 K 2 ? c3 K 3 ]

类似二阶的推导过程,将 K2, K3按二元函数展 开,使Tn+1= O(h4),得
?c1 ? c 2 ? c 3 ? 1 ? ?? 2 ? ? 21 ?? 3 ? ? 31 ? ? 32 ? ?c ? ? c ? ? 1 3 3 ? 2 2 2 ? ?c ? 2 ? c ? 2 ? 1 3 3 ? 2 2 3 ? 1 ?c 3 ? 2 ? 32 ? 6 ?

方程有8个未知数,解不唯一。 满足该条件的公式统称为三阶 R-K公式。 其中一个常用公式为:
h ? ? y n?1 ? y n ? 6 ( K 1 ? 4 K 2 ? K 3 ) ? ? K 1 ? f ( xn , yn ) ? ?K ? f ( x ? h , y ? h K ) n n 1 ? 2 2 2 ? K ? f ( x ? h, y ? hK ? 2hK ) n n 1 2 ? 3

当 r=4 时,利用相同的推导过程,经过较复杂的 计算,可以得出四阶R-K公式的成立条件。 下列经典公式是其中常用的一个:
h ? ? y n?1 ? y n ? 6 ( K 1 ? 2 K 2 ? 2 K 3 ? K 4 ) ? ? K 1 ? f ( xn , yn ) ? h h ? ? K 2 ? f ( xn ? , yn ? K 1 ) 2 2 ? h h ? ? K 3 ? f ( xn ? 2 , yn ? 2 K 2 ) ? ? ? K 4 ? f ( x n ? h, y n ? hK 3 )

四、一阶常微分方程组和高阶微分方程的 数值解法简介
一阶常微分方程组的数值解法: 下列包含多个一阶常微分方程的初值问题:
? ? f 1 ( x , y 1 , y 2 ,? , y n ) ? y1 ? ? ? y 2 ? f 2 ( x , y 1 , y 2 ,? , y n ) ? ? ??? ? n ? f n ( x , y 1 , y 2 ,? , y n ) ? y? y1 ( x 0 ) ? y10 y 2 ( x 0 ) ? y 20 y n ( x 0 ) ? y n0

称为一阶常微分方程组的初值问题。

引进向量记号:
? y1 ? ? ? ? ? y2 ? y?? ? ? ? ? ?y ? ? n?
? ? ? ? ? f ( x, y) ? ? ? ? ?

? f 1 ( x, y) ? ? ? f 2 ( x, y) ? ? ? ? ? f n ( x, y) ? ?

? y1 ( x 0 ) ? ? y1 0 ? ? ? ? ? ? ? y2 ( x0 ) ? ? y2 0 ? y0 ? ? ?? ? ? ?? ? ? ? ? ? y ( x )? ? y ? ? n 0 ? ? n0 ?

则上述一阶常微分方程组的初值问题化为矩阵 ? ? ? 形式: ? y ? ? f ( x, y)
?? ? ? y( x 0 ) ? y 0

它在形式上跟单个微分方程的初值问题形式完 全相同,只是函数变成了向量函数。故前面介 绍的一切数值方法都适用,只要把函数换成向 量函数即可。

四阶龙格-库塔公式
数值求解:

? y? ? f ( x, y) ? ? y( x 0 ) ? y 0

x ? x0

k1=f(xn,yn), k2=f(xn+0.5h,yn+0.5hk1) k3=f(xn+0.5h,yn+0.5hk2), k4=f(xn+h,yn+hk3)

yn+1= yn+ h[k1+2k2+2k3+k4]/6 (n = 0, 1, 2, · · · · · ,N)

狐狸和野兔问题——常微分方程组

在一个生物圈中,只有野兔和狐狸两种动物,记y1 为野兔数量, y2 为狐狸数量. 设有足够的食物供 野兔生存,而狐狸只以野兔为食物. 模型如下 ? dy1 自变量x ∈(0, 15), ? y ? 0 . 01 y y 1 1 2 ? dx 初值条件为: ? dy ? 2 ? ? y 2 ? 0.02 y1 y 2 y (0)=20, y (0)=20 1 2 ? dx
? y1 ? y?? ? ? y2 ?

? y1 ? 0.01y1 y 2 ? f ( x, y) ? ? ? ?? y 2 ? 0.02 y1 y 2 ?

求解常微分方程: function z=lot1(x,y) z(1,:)=y(1)-0.01*y(1).*y(2); (1)定义函数; z(2,:)=-y(2)+0.02*y(1).*y(2); (2)调用ode23 t0=0;t1=15; %设置区域 y0=[20 20]; %给定初值 [t,y]=ode23(‘lot1’,t0,t1,y0); %求解 plot(t,y) %绘图
400 300 200 100 0

0

5

10

15

一、高阶常微分方程的数值解法:
可化为一阶常微分方程组求解。 例如,二阶常微分方程初值问题
? y ?? ? f ( x, y, y ?) ? ? y ( x0 ) ? y 0 ? y ?( x ) ? y ? 0 0 ?

引进新的变量,令z=y`,即可将上述二阶方程 化为如下的一阶方程组的初值问题:

? z ? ? f ( x, y, z ), ? ? y ? ? z,

? z ( x0 ) ? y 0 y ( x0 ) ? y 0

例 求下列高阶微分方程的数值解:
y??? ? 3 y?? ? y?y ? 0
y(0) ? 0, y?(0) ? 1, y??(0) ? ?1
(0 ? x ? 2 )

解:显然 假设

y??? ? 3 y?? ? y?y
y1 ? y y2 ? y? y3 ? y??

? ? y2 则 ? y1

? y? ? y ? 2 3 ? ? ? 3 y3 ? y2 y1 ? y3 ? ? y1 (0) ? 0, y2 (0) ? 1, y3 (0) ? ?1

即二阶问题化为微分方程组的初值问题。

五、二阶常微分方程边值问题数值方法
考虑方程:
y?? ? f ( x, y, y?), a ? x ? b
( 5.1)

结合下述三种边界条件之一:
( 5 .2 ) y( a ) ? ? , y( b ) ? ? ; ( 5 .3 ) y?(a ) ? ? , y?(b) ? ? ; ( 5 .4 ) y?(a) ? ? 0 y(a) ? ?1 , y?(b) ? ? 0 y(b) ? ?1;

(5.4)式中 ? 0 ? 0, ? 0 ? 0,? 0 ? ? 0 ? 0 它们分别称为第一、第二、第三边界问题。 边界问题的解法:打靶法、有限差分法

打靶法
基本思路:将边值问题转化为初值问题考虑。或者 说适当选择初始值使初值问题的解满足边值条件。然 后用求解初值问题的任一种有效的数值方法求解。

以第一边界条件为例
考虑边值问题:? y ?? ? f ( x, y, y ?), a ? x ? b ? ? y( a ) ? ? , y ( b ) ? ? 取y0=? , 考虑初值问题

? y ?? ? f ( x, y, y ?), a ? x ? b ? ? ? y(a ) ? y 0 , y ?(a ) ? y 0

( 5 .5 )

? y0

待定,由数值解法求解(5.5)得到在 x ? xn ? b 处的 ? 解 y? n ? yn ( y0 , y0 ) ,

? ) ? g( y0 ?) yn ? yn ( y0 , y0
? ) ? ? , 则得所求数值解。当 若 g( y0
? ) ? ?时, g( y0

?)? ? ? 0 求解非线性方程 g( y0

( 5 .6 )

该方程可用二分法、正割法或Newton法等来求解。 ?, 若求得 y0
? ) ? ? ? ? ,这里 ? ? 0 为给定允许误差界,就 使 g( y0 停止迭代改进, 输出作为数值解。

对第二类边界问题,也可转化为考虑初值问题 (5.5), ? ? ? , 以 y0 为待定参数。 取 y0 对第三类边界问题,仍可转化为考虑初值问题 ? ? ?1 ? ?0 y0 ,以 y0 为待定参数。 (5.5),取 y0

有限差分法
b?a , xi ? a ? ih, i ? 0,1,?, n ? 1, 设在 将区间[a,b]进行等分:h ? n?1 x ? xi , i ? 0,1,?, n ? 1处的数值解为 yi 。用中心差分近似微分,即 yi ? 1 ? yi ? 1 y? ? , i ? 1,2,?, n i ? 2h ( 5 .7 ) ? yi ? 1 ? 2 yi ? yi ? 1 ? ? y ? , i ? 1,2,?, n ? i 2 h 二阶中心差商 则离散化成差分方程 yi ? 1 ? yi ? 1 2 yi ? 1 ? 2 yi ? yi ? 1 ? h f ( x i , yi , ), i ? 1,2,?, n (5.8) 2h

?) ?( ( y (a a) )? ?? ?,, yy (b b) )? ?? ? ? a ? ? 对应的边界条件也离散成 y?(y 0 y(a ) ? ?1 , y (b) ? ? 0 y(b) ? ?1; ( 5 .9 ) 第一边界问题: y0 ? ? , yn?1 ? ? (5.10) 第二边界问题: y1 ? y0 ? h? , yn?1 ? yn ? h? 第三边界问题: y1 ? (1 ? ? 0h) y0 ? ?1h, (5.11) (1 ? ? 0h) yn?1 ? yn ? ?1h 若 f ( x, y, y?) 是 y , y? 的线性函数时,f 可写成 f ( x, y, y?) ? p( x ) y?( x ) ? q ( x ) y( x ) ? r ( x ) 其中 p( x ), q( x ), r ( x )为已知函数,则由常微分方程的理论知,通过 变量替换总可以消去方程中的 y? 项,不妨设变换后的方程为 ? y??( x ) ? q( x ) y( x ) ? r ( x ) ? ? y( a ) ? ? , y( b ) ? ? ? y i ? 1 ? 2 y i ? y i ? 1 ? qi yi ? ri 2 则近似差分方程成离散差分方程为 ? ? h ? y0 ? ? , yn ? 1 ? ? ? q ? q ( x ), r ? r ( x ), i ? 1 , 2 , ? , n . 其中 i i i i

将以上方程合并同类项整理得方程组 y0 ? ? ? 2 2 y ? ( 2 ? q h ) y ? y ? r h i i i ?1 i ? i ?1 ? yn ? 1 ? ?

? yi ? 1 ? 2 yi ? yi ? 1 ? ? qi yi ? ri 2 ? h ? 0 ? ? (iy? 1,? 2,, ?, n)yn ? 1 ? ?

其中只要 qi ? 0 ,则方程组的系数矩阵为弱对角占优的三对角阵, 方程组为三对角线方程组,可以用追赶法求解。而且还有误差估 M 2 计: Ri ? y( xi ) ? yi ? h ( xi ? a )( b ? xi )
max y ( 4 ) ( x ) 。 其中 M ? x ?[ a , b ] 若 f ( x, y, y?) 不是 y , y? 的线性函数时,所得方程组是非线性 组,可以用解非线性方程组的方法求解。例如,可用简单迭代法、 Newton迭代法求解。
24


网站首页 | 网站地图 | 学霸百科 | 新词新语
All rights reserved Powered by 大学生考试网 9299.net
文档资料库内容来自网络,如有侵犯请联系客服。zhit325@qq.com