欧拉法与龙格库塔法解常微分方程 设一阶微分方程及初值为：    y    ( y x  0 dy dx )   ( , f x y ) y 0 (1) 欧拉法 过点(x0,y0)以 y’(x0)=f(x0,y0)为斜率作切线，切线方程： y  y 0  ( , f x y 0 0 )( x  x 0 ) 欧拉法即是 f(x,y)在(x0,y0)处的一阶泰勒展开，即有： y n 1   y n  ( , f x y n n )( x n 1   x n ) 设 x0,x1,x2,…xn 的步长为 h，则欧拉法求解的公式可表示为： x n y    1  n 1    x n y n   h ( , f x y n ) h  n 欧拉法具有一阶精度，其局部阶段误差是关于步长的二阶无穷小量。 (2) 改进欧拉法 由微分中值定理，改进欧拉方程为： y n 1   y n   ( , f x y n n )  ( f x n 1  , y n 1   ) h 2  预报值 :   校正值 :  y n 1   y n  )  n ( f x n 1  , y n 1   ) ) h  n n , ( f x y h  2 n , ( f x y 由于等式两边都存有 yn+1 未知量，这种形式称为隐式形式。因为是近似计算， 可以由欧拉公式求得 yn+1 的一个近似值（预报值），然后将其带入公式中再进行 计算得到一个 yn+1 值（校正值）。   y n 1  y n 写成一个公式即为： y n 1   y n   h 2 , ( f x y n n )  ( f x n 1  , y n  ( , f x y n n )  h )  改进欧拉法具有二阶精度，其局部阶段误差是关于步长的三阶无穷小量。 (3) 龙格库塔法 根据微分中值定理，在(xn,xn+1)区间存在一个ζ满足： y  ( )   y n x n 1  1    n y x n   ( , x x n n 1  ) 称ζ处的斜率 f(ζ,y(ζ))为平均斜率 K，则有： K h     1n y y n 显然，取 xn 处的斜率 f(xn,yn)为 K 时即为欧拉法，取 xn 和 xn+1 处的斜率均 值为 K 时为改进欧拉法。若在[xn,xn+1]区间多取几个点的斜率，以其某种加权均 

值作为 K，可进一步提升 yn+1 精度。 在[xn,xn+1]区间取 m 个点，其中第 j 个点及对应的斜率 Kj 为： x j  x n   j h (0   j  1,1   j m K j  ( , f x y j j )  ( f x n   j , h y n  h (  jj K jj )) )  1  j jj 1       每个斜率 Kj 的加权值为λj，则平均斜率 K 及 yn+1 估值为： m  j 1  (  j K ) j  j  1 m   K         y 1  1  n j  y n  K h   y n  h m  j 1  ( K  j j ) 此即为龙格库塔算法，选取恰当的αj,βj 和λj 可使算法精度尽可能提升，算法 阶数由取点数 m 决定，阶数越高精度越高。 二阶龙格库塔公式：        y n K 1 K 2 y hK   2 1 n  , ) ( f x y  n h 2 ( f x   n n , y n  h 2 K 1 ) 三阶龙格库塔公式：          四阶龙格库塔公式：            K K K Matlab 代码  1  y n K 1 K K 2 3    K 1 n y  h  n 6 ) , ( f x y n h 2 , h y ( f x ( f x   , n n y n  4 K 2  K 3  n )  h K 1 2 ( h K   1  2 K )) 2  2 K 2  2 K 3  K 4  1  1  y n K 1     2 3 4 n y  K h  n 6 ) ( , f x y n h 2 h 2 , h y n ( f x ( f x ( f x    , , n n y y n n ) )   h 2 h 2 hK K 1 K 2 ) 3  n 

%% 初始条件 syms x y; f(x,y) = y*cos(x); % 条件 1：方程导数 y'=f(x,y) x0 = 0; y0 = 1; % 条件 2：初值 y(0)=1 % 原方程为 y=e^(sinx) h = 0.2; x = x0:h:10; y = exp(sin(x)); len = length(x); % 步长 % x 取值范围 % 真实的 y 值（待求） %% 二阶龙格库塔 y2 = zeros(size(x)); % 初始化 y y2(1) = y0; for ii = 2:len K1 = f(x(ii-1),y2(ii-1)); K2 = f(x(ii-1)+h/2,y2(ii-1)+h*K1/2); y2(ii) = y2(ii-1) + h*K2; K1 = f(x(ii-1),y4(ii-1)); K2 = f(x(ii-1)+h/2,y4(ii-1)+h*K1/2); K3 = f(x(ii-1)+h/2,y4(ii-1)+h*K2/2); K4 = f(x(ii-1)+h,y4(ii-1)+h*K3); y4(ii) = y3(ii-1) + h*(K1+2*K2+2*K3+K4)/6; end %% 绘图 figure %% 三阶龙格库塔 y3 = zeros(size(x)); y3(1) = y0; for ii = 2:len end end %% 三阶龙格库塔 y4 = zeros(size(x)); y4(1) = y0; for ii = 2:len K1 = f(x(ii-1),y3(ii-1)); K2 = f(x(ii-1)+h/2,y3(ii-1)+h*K1/2); K3 = f(x(ii-1)+h,y3(ii-1)+h*(K2*2-K1)); y3(ii) = y3(ii-1) + h*(K1+4*K2+K3)/6; 

plot(x,y,x,y2,x,y3,x,y4) title('\fontname{微软雅黑}\fontsize{10}龙格库塔求解值'); xlabel('\fontname{微软雅黑}\fontsize{10}x'); ylabel('\fontname{微软雅黑}\fontsize{10}y'); legend('原值','2 阶值','3 阶值','4 阶值') figure plot(x,y2-y,x,y3-y,x,y4-y) title('\fontname{微软雅黑}\fontsize{10}龙格库塔误差'); xlabel('\fontname{微软雅黑}\fontsize{10}x') ylabel('\fontname{微软雅黑}\fontsize{10}'); legend('2 阶误差','3 阶误差','4 阶误差')