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

House

+ Y5 b3 l. \! H6 j7 L, y
MATLAB源程序代码分享：MATLAB实现四阶龙格库塔法求解常微分方程组
function F=f(t,Y)
% 定义待求的微分方程组
x=Y(1);
/ J: `" c: I# [% w$ W( |4 a9 G( B2 R5 Ey=Y(2);
f1=-x^2*(5-y);
f2=y*(4-3*x);
F=[f1;f2];
4 z: }( t, X; T8 c  Z4 v, _3 Yend




%% 定义计算的步长, 设置变量的初始值
; R8 S# s) e2 O2 \9 i2 C$ mclear;clc;close all
Delta=0.001;     % 步长
4 q% d* t( j1 C  M* v8 It=0: Delta :8;     % 定义自变量 t
9 D; o+ u2 c; j  sn=length(t);
0 P& q4 C! B, d+ n+ d' a/ {Y(:,1)=[0.5;3];  % 变量 x y 的初始值
; D6 M8 |' X0 f( [6 r; n
%% 自定义龙格库塔法, 求解微分方程组
- [( A, _# a) Z; Z+ O' T0 pfor k=1: n-1
; W8 w% `0 {! Y0 u7 n    z1=f(t(k),Y(:,k));
    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
& R0 l4 q& f( z) o) Z, T: K    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
5 \& c* v8 P. ^( `/ v    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
1 o$ P& b# q4 x; @end
( C5 s! ~3 e/ h! Ax=Y(1,: );
y=Y(2,: );
: b; g, U  q+ a7 I! ^/ [3 ]1 z
%% 绘制 x y 的求解结果
0 ]3 m* {! B( ~# q, Lfigure
8 t! R1 q3 N: s) N- U8 E# S* `set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]);  % 设置 figure 窗口的位置和尺寸
6 ?8 E9 A5 m  o1 M$ u) _* t4 bsubplot(1,2,1)
plot(t,x,'b')
$ h3 D/ J- R! [xlabel('t')
* q  J: B) C: r! \% kylabel('x')
2 a+ s1 V! u3 w3 A" _4 ]
subplot(1,2,2)
' S0 Z8 Z9 j4 G& xplot(t,y,'r')
xlabel('t')
) O5 |, e1 V2 E' l0 N' q6 \ylabel('y')
3 F& V. w, G# N0 c) H
0 c) Z' a% d, V6 C9 A

helendcany

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