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

House

MATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组
& v" R! r2 C& M3 wfunction F=f(t,Y)
% 定义待求的微分方程组
x=Y(1);
9 h  e; a: y) ~  ]y=Y(2);
f1=-x^2*(5-y);
$ e' ]4 Y2 q1 ^0 D2 T* wf2=y*(4-3*x);
F=[f1;f2];
& u" o1 K- J+ I- [" C+ J" S1 a( Yend
* ^2 q6 p  Z' m+ n1 c
8 c  k3 S8 h3 i  W; N; A
0 W! o; R5 B/ t( @' l

1 o2 r9 s; E1 N# f/ ~, |* X6 W%% 定义计算的步长, 设置变量的初始值
clear;clc;close all
Delta=0.001;     % 步长
2 l& l  l0 m( a, `; g& At=0: Delta :8;     % 定义自变量 t
! r5 E$ m4 v' p7 e! [% q) O/ ]# }n=length(t);
Y(:,1)=[0.5;3];  % 变量 x y 的初始值

3 H( x0 F( r! L. f8 M%% 自定义龙格库塔法, 求解微分方程组
- ^. k* l7 v0 y6 H/ {+ S+ Sfor k=1: n-1
    z1=f(t(k),Y(:,k));
5 s0 h9 U! E9 s( S8 N7 D    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
9 X, I+ z, k( O  J( l    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
end
) I1 w& V1 o' {3 c- g) R3 }! Cx=Y(1,: );
y=Y(2,: );
$ ^8 Y: L( J5 j( ]6 W  @
' l5 z3 p' z. a% g  X5 t8 l4 r7 ?%% 绘制 x y 的求解结果
- L) v+ Z2 ?+ o$ A3 Rfigure
7 R! {- ]+ h& o4 P- Bset(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
subplot(1,2,1)
plot(t,x,'b')
) i0 @9 c) H. J+ n0 Pxlabel('t')
ylabel('x')

subplot(1,2,2)
/ B/ W1 b9 J1 y5 F9 eplot(t,y,'r')
2 S5 R' C" _3 k6 D, y5 s/ \% \xlabel('t')
2 Q0 X0 c1 L* C6 k0 G8 d, Jylabel('y')

, M) B9 n* i; P6 q8 e: Q

helendcany

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