第8讲 常微分方程数值解法

Ordinary Differential Equations

•一阶常微分方程的初值问题：

•节点：x 1<x 2< … <x n

•步长

为常数 ⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy

1--=i i x x h

•一 欧拉方法（折线法）

y i +1=y i +h f (x i ,y i ) (i =0,1, …, n-1) 优点：计算简单。 缺点：一阶精度。 •二 改进的欧拉方法

)

(2

1

)

,(),(11c p i p i i c i i i p y y y y x hf y y y x hf y y +=+=+=++

•改进的欧拉公式可改写为

它每一步计算f (x,y )两次，截断误差为O(h 3)

⎪⎪⎩⎪

⎪⎨⎧++==++=+),(),()(2121211

hk y h x f k y x f k k k h y y i i i i i i

y dt

dy

-=1

0,1)0(≤≤=t y t

e

y -=精确解：

function [t,y] = Heun(ode,tspan,h,y0) t = (tspan(1):h:tspan(end))'; n = length(t); y = y0*ones(n,1); for i=2:n

k1 = feval(ode,t(i-1),y(i-1)); k2 = feval(ode,t(i),y(i-1)+h*k1); y(i) = y(i-1)+h*(k1+k2)/2; end

•三 龙格—库塔法(Runge-Kutta)

欧拉公式可改写为

它每一步计算 f (x i ,y i ) 一次，截断误差为O(h 2)

⎩

⎨⎧=+=+),(111i

i i i y x f k hk y y

•标准四阶龙格—库塔公式

每一步计算 f (x, y ) 四次，截断误差为O(h 5)

⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨

⎧

++=++=++==++++=+),()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 i i i i i i i i i i 0 1/2 1/2 1/2 0 1/2 1

1

1/6

2/6

2/6

1/6

.

)()),(),(,()(;)()),(),(,()('

'

βα====b v x v x u x g x v a u x v x u x f x u ()

()

4321143211226

226c c c c h

v v k k k k h

u u i i i i ++++=++++=++)

,,(),,()

2

1

,21,2()

21

,21,2()

21

,21,2()21,21,2()

,,(),,(33433422322311211211hc v hk u h x g c hc v hk u h x f k hc v hk u h x g c hc v hk u h x f k hc v hk u h x g c hc v hk u h x f k v u x g c v u x f k i i i i i i i i i i i i i i i i i i i i i i i i +++=+++=+++=+++=+++=+++

===对于两个分量的一阶常微分方程组

用经典4阶 Runge-Kutta 法求解的格式为

,

1

1∑=++=n

j j j i i k b h y y .

,,3,2,),

,(),

,(1

1

1

1

1n j a c k a h y h c x f k y x f k j m jm j j m m jm i j i j i i ==++==∑∑-=-=n 级显式Runge-Kutta 方法的一般计算格式：

其中

.

1,,3,2),51623(12

211-=+-+=--+N i f f f h

y y i i i i i .

1,,4,3),9375955(24

3211-=-+-+=---+N i f f f f h

y y i i i i i i .

1,,3,2),5199(24

2111

-=+-++=--++N i f f f f h

y y i i i i i i Adams 外插公式（Adams-Bashforth 公式）是一类 k +1 步 k +1 阶显式方法 三步法（k =2）,

四步法（k =3）,

Adams 内插公式（Adams-Moulton 公式）是一类 k +1 步 k +2 阶隐式方法 三步法（k =2）,

Adams 预估-校正方法（Adams-Bashforth-Moulton 公式） 一般取四步外插法与三步内插法结合。

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define TRUE 1

main()

{

int nstep_pr, j, k;

float h, hh, k1, k2, k3, k4, t_old, t_limit, t_mid, t_new, t_pr, y, ya, yn; double fun();

printf( "\n Fourth-Order Runge-Kutta Scheme \n" );

while(TRUE){

printf( "Interval of t for printing ?\n" );

scanf( "%f", &t_pr );

printf( "Number of steps in one printing interval?\n" );

scanf( "%d", &nstep_pr );

printf( "Maximum t?\n" );

scanf( "%f", &t_limit );

y = 1.0; /* Setting the initial value of the solution */

h = t_pr/nstep_pr;

printf( "h=%g \n", h );

t_new = 0; /* Time is initialized. */

hh = h/2;

printf( "--------------------------------------\n" );

printf( " t y\n" );

printf( "--------------------------------------\n" );

printf( " %12.5f %15.6e \n", t_new, y );