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

House

/ e( g) G# ~0 `% J, f
' M% F8 ]! v- A8 QMATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组
5 t/ B, z* t1 y( r( R, e& G$ Cfunction F=f(t,Y)
% 定义待求的微分方程组
x=Y(1);
y=Y(2);
1 T+ X' @( X$ |# l( r' bf1=-x^2*(5-y);
8 ~" k1 {+ ?& Zf2=y*(4-3*x);
F=[f1;f2];
end
$ V8 l5 F( d. u& v4 ]  f
: o4 x9 g8 r/ H! P* C

' p4 _' x1 J- ]+ i. s# u% A
' V: U8 |0 W' \& \. d6 s3 t5 M%% 定义计算的步长, 设置变量的初始值
8 T! h9 q8 t- Z1 Aclear;clc;close all
: ?; _4 B, Y8 E0 zDelta=0.001;     % 步长
/ }5 I/ M$ t. H( S# l& z% j" m3 G1 K3 |t=0: Delta :8;     % 定义自变量 t
0 Y# ~2 J4 s/ i! t! w1 Z5 ~n=length(t);
# B' a. W4 D! {, ]$ \Y(:,1)=[0.5;3];  % 变量 x y 的初始值

%% 自定义龙格库塔法, 求解微分方程组
& k4 X; W- Y* p3 d0 C8 G7 Mfor k=1: n-1
    z1=f(t(k),Y(:,k));
    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
/ v* d- u2 M# l0 B( {    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
4 R% m/ ]) A1 ]: P+ y    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
  F% Y, t9 w6 |: a    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
1 O% @* c+ b4 L% @, |- Xend
- u* k$ I! D1 C$ wx=Y(1,: );
) R' f  R! x; Z  n9 @- Y0 e+ gy=Y(2,: );
" W+ _3 i8 ^4 ]; N
%% 绘制 x y 的求解结果
& ~" Q' \3 n3 _figure
& H! t( p3 `; R  sset(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
subplot(1,2,1)
3 _* ?3 F/ A% }; r- kplot(t,x,'b')
2 M) f+ ^8 P7 t$ j' i/ v; \xlabel('t')
ylabel('x')

subplot(1,2,2)
plot(t,y,'r')
xlabel('t')
ylabel('y')
, I: \' d" {% R
6 Z, M" M( r" K0 i

helendcany

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