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

House

" R: R/ h& s' ~  p; h1 B& _  ~$ G
6 X% R% Y3 E: w% R5 t- U+ dMATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组
function F=f(t,Y)
. a! i4 g! u. Z* r$ `; L  y, x8 u( D% 定义待求的微分方程组
2 b. g$ s3 H$ i( w1 X" Z1 Bx=Y(1);
y=Y(2);
: p; ]; _, k" Y& v0 [5 L) b3 Ef1=-x^2*(5-y);
$ u- K* k0 x. `, b3 Kf2=y*(4-3*x);
, h: H( q3 X' Y% X! j1 HF=[f1;f2];
end
! h- Z/ r' g( @. ~
7 s* K0 a+ ~1 v* W3 z- |8 p


%% 定义计算的步长, 设置变量的初始值
" a8 V. N1 r, z2 yclear;clc;close all
Delta=0.001;     % 步长
t=0: Delta :8;     % 定义自变量 t
% B6 q- i3 l( F$ Q1 _0 x1 vn=length(t);
. f" [" o, r1 g' h2 S6 L) hY(:,1)=[0.5;3];  % 变量 x y 的初始值

%% 自定义龙格库塔法, 求解微分方程组
# C+ i7 W# n# p1 K% L( \- |for k=1: n-1
) }+ b3 J: K) p, l    z1=f(t(k),Y(:,k));
6 T/ ?4 O  J1 W5 s# Q( g    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
9 z+ t) G4 x; X    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
& y3 A- k- N8 _0 S    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
end
x=Y(1,: );
+ C1 O7 s: L. p& [. `, @# e. `y=Y(2,: );
4 p3 T6 N3 F3 D0 I& Z
, O8 C3 `( t1 x! O%% 绘制 x y 的求解结果
figure
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
9 L6 v# o7 j9 _: j. \subplot(1,2,1)
plot(t,x,'b')
4 M4 o$ r4 }  n7 }7 wxlabel('t')
' |, m" K# R' O" \% Tylabel('x')

subplot(1,2,2)
# M) z: P8 L4 ]- K& n5 bplot(t,y,'r')
xlabel('t')
$ B+ q# q- d) j+ ?7 u. Sylabel('y')
- L. l. A5 p* N9 d  ~6 X& _8 Y

helendcany

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