MATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组

House

0 p% T7 D* `, c# T% @% e6 P$ G
MATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组
function F=f(t,Y)
! w* o' C) k1 x  O! r" ~! o# t% 定义待求的微分方程组
$ S2 D* I9 ~0 \) t, M3 H4 Ex=Y(1);
y=Y(2);
f1=-x^2*(5-y);
f2=y*(4-3*x);
& h4 N* k( d) U& p* }F=[f1;f2];
; a) X/ L8 t# wend
4 @  P7 D3 H' m6 E! ~; r


7 C2 p9 e) @7 ^. c8 @
%% 定义计算的步长, 设置变量的初始值
1 Z' Q6 c6 u( u6 A) j3 a& n" o7 ]clear;clc;close all
, J  C% ?/ Z1 qDelta=0.001;     % 步长
t=0: Delta :8;     % 定义自变量 t
n=length(t);
Y(:,1)=[0.5;3];  % 变量 x y 的初始值

%% 自定义龙格库塔法, 求解微分方程组
for k=1: n-1
- u! z& q( r, Y    z1=f(t(k),Y(:,k));
* I0 I% G7 I  `5 Q' c( t: r  B    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
  s+ S  G$ [4 F. B    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
% m( h4 _" V/ m- }; K. S    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
4 r/ W- `) l1 f3 @, uend
& Z1 l3 p/ ~" Y- P& E' Px=Y(1,: );
y=Y(2,: );

- h2 @6 Z1 z! v( |& y. `, K; n%% 绘制 x y 的求解结果
$ t# P$ B3 p" j; H; M0 _6 X. T' Vfigure
- |* H- `! N8 n1 @. _$ _% Uset(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
subplot(1,2,1)
plot(t,x,'b')
xlabel('t')
/ P# T+ x# |+ hylabel('x')
5 H! @* r. J3 k2 X" O/ m* M8 {9 E
subplot(1,2,2)
plot(t,y,'r')
xlabel('t')
7 |0 ]6 [7 i: o$ @3 N9 qylabel('y')
( V5 O1 Q. h! E8 j1 P- }
! x1 z7 z7 T7 a) z/ J7 y+ D% X6 X( D0 w+ w

helendcany

MATLAB实现四阶龙格库塔法求解常微分方程组