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

#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 );
n 级显式Runge-Kutta 方法的一般计算格式：
yi 1 yi h b j k j ,
其中
k1 f ( x i , y i ),
j 1
n
k j f ( x i c j h, y i h a jm k m ),
m 1
j 1
c j a jm , j 2, 3, , n.
• 三
龙格—库塔法(Runge-Kutta)
欧拉公式可改写为
yi 1 yi hk1 k f ( x , y ) i i 1
它每一步计算 f (xi,yi) 一次，截断误差为O(h2)
• 标准四阶龙格—库塔公式
yi 1 h yi (k1 2k 2 2k3 k 4 ) 6 k1 f ( xi , yi ) h h k 2 f ( xi , yi k1 ) 2 2 h h k3 f ( xi , yi k 2 )Байду номын сангаас2 2 k 4 f ( xi h, yi hk3 )
y p y i hf ( x i , y i ) yc y i hf ( x i 1 , y p ) 1 y i 1 ( y p yc ) 2
• 改进的欧拉公式可改写为
h y y ( k k ) i 1 i 1 2 2 k1 f ( xi , yi ) k f ( x h, y hk ) i i 1 2
0 1/2 1/2 1/2 0 1/2
1
0
1/6
0
2/6
1
2/6 1/6
每一步计算 f (x, y) 四次，截断误差为O(h5)
对于两个分量的一阶常微分方程组
u' ( x ) f ( x , u( x ), v ( x )), u(a ) ; v ' ( x ) g( x , u( x ), v ( x )), v (b) .
它每一步计算f(x,y)两次，截断误差为O(h3)
dy y dt
y(0) 1, 0 t 1
t y e 精确解：
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
• 一阶常微分方程的初值问题：
dy f ( x, y) dx y ( x 0 ) y0
• 节点：x1<x2< … <xn • 步长 h xi xi 1 为常数
• 一
欧拉方法（折线法） yi+1=yi+h f(xi,yi) (i =0,1, …, n-1) 优点：计算简单。 缺点：一阶精度。 • 二 改进的欧拉方法
用经典4阶 Runge-Kutta 法求解的格式为
h ui 1 ui k1 2k 2 2k 3 k 4 6 h v i 1 v i c1 2c 2 2c3 c4 6
k1 f ( x i , ui , v i ) c1 g ( x i , ui , v i ) h 1 1 , ui hk1 , v i hc1 ) 2 2 2 h 1 1 c2 g( x i , ui hk1 , v i hc1 ) 2 2 2 h 1 1 k 3 f ( x i , ui hk 2 , v i hc2 ) 2 2 2 h 1 1 c3 g( x i , ui hk 2 , v i hc2 ) 2 2 2 k4 f ( x i h, ui hk 3 , v i hc3 ) k2 f ( xi c4 g( x i h, ui hk 3 , v i hc3 )
Adams 内插公式（Adams-Moulton 公式）是一类 k+1 步 k+2 阶隐式方法 三步法（k=2）,
yi 1
h yi (9 f i 1 19 f i 5 f i 1 f i 2 ), i 2, 3, , N 1. 24
Adams 预估-校正方法（Adams-Bashforth-Moulton 公式） 一般取四步外插法与三步内插法结合。
j 1
Adams 外插公式（Adams-Bashforth 公式）是一类 k+1 步 k+1 阶显式方法 三步法（k=2）,
m 1
h yi 1 yi ( 23 f i 16 f i 1 5 f i 2 ), i 2, 3, , N 1. 12 四步法（k=3）, h yi 1 yi (55 f i 59 f i 1 37 f i 2 9 f i 3 ), i 3, 4, , N 1. 24