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

House

! h/ \/ a0 B/ r
MATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组
function F=f(t,Y)
& ]& R+ {8 k8 U5 k% 定义待求的微分方程组
x=Y(1);
8 N% X0 C) |* |* t8 Wy=Y(2);
f1=-x^2*(5-y);
f2=y*(4-3*x);
: Y8 s" r7 N0 aF=[f1;f2];
, _: N6 U: m6 v; A* s- nend
: y/ H" u# ^2 d5 |, H, \8 V
: m: k0 B3 z5 B7 p

) m* G+ l) u$ _
%% 定义计算的步长, 设置变量的初始值
clear;clc;close all
Delta=0.001;     % 步长
t=0: Delta :8;     % 定义自变量 t
n=length(t);
Y(:,1)=[0.5;3];  % 变量 x y 的初始值
3 r8 Q8 v- H* X$ A9 w
%% 自定义龙格库塔法, 求解微分方程组
( L( e+ X4 I) C; e' W/ a0 K9 Wfor k=1: n-1
  m: @9 F* {6 h2 ?. N1 z    z1=f(t(k),Y(:,k));
+ O; r/ D+ Q2 _$ j5 ?# Y3 q2 h    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
& N: D% x1 |! g) }7 i9 j6 Q    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
9 f; ]3 X4 ^" f: H    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
) e; ?9 M  z  Pend
x=Y(1,: );
1 t5 h" [* J, X# z& @y=Y(2,: );
8 I& {! J! j8 `
$ k1 b, s* c5 b: `, f% }%% 绘制 x y 的求解结果
$ a& A0 c* V! @  m% M2 K9 v1 P2 mfigure
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
subplot(1,2,1)
plot(t,x,'b')
xlabel('t')
  D' k6 y5 T" e  I+ p* Uylabel('x')

subplot(1,2,2)
plot(t,y,'r')
xlabel('t')
0 U( S0 p" P5 tylabel('y')

helendcany

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