用四阶龙格库塔法求解矩阵微分方程

global R M U
syms t
R=[
00247,0,0,0,0,0;
0,00247,0,0,0,0;
0,0,00247,0,0,0;
0,0,0,00193,-00193,0;
0,0,0,0,00193,-00193;
0,0,0,1,1,1
];
M=[
0,0,0,-15727/10000sin(5/12pi+80pit)pi,15727/10000cos(1/4pi+80pit)pi,15727/10000sin(1/12pi+80pit)pi;
0,0,0,15727/10000sin(1/12pi+80pit)pi, -15727/10000sin(5/12pi+80pit)pi, 15727/10000cos(1/4pi+80pit)pi;
0,0,0, 15727/10000cos(1/4pi+80pit)pi,15727/10000sin(1/12pi+80pit)pi, -15727/10000sin(5/12pi+80pit)pi;
-15727/10000sin(5/12pi+80pit),15727/10000sin(1/12pi+80pit)pi, 15727/10000cos(1/4pi+80pit)pi,0,0,0;
15727/10000cos(1/4pi+80pit)pi, -15727/10000sin(5/12pi+80pit)pi,15727/10000sin(1/12pi+80pit)pi,0,0,0;
15727/10000sin(1/12pi+80pit)pi, 15727/10000cos(1/4pi+80pit)pi, -15727/10000sin(5/12pi+80pit)pi,0,0,0
];
U=[
380cos(100pit);
380cos(100pit-2pi/3);
380cos(100pit+2pi/3);
0;
0;
0
];
DM=diff(M,t);
%%%要调节的参数在这里
%%注意，你的M有点奇异，计算很快发散掉了，你检察一下相关的参数吧。
%%det(subs(M,t,0))
I0=[0;0;0;0;0;0];
tstart=0;
tend=01;
dt=01;
%%%end
tout=tstart:dt:tend;
n=length(tout);
M_t_dt=subs(M,t,tstart);
U_t_dt=subs(U,t,tstart);
DM_t_dt=subs(DM,t,tstart);
II=I0;
for i=1:n-1
tt=tout(i);

M_t=M_t_dt;
U_t=U_t_dt;
DM_t=DM_t_dt;

M_t_dt_2=subs(M,t,tt+dt/2);
U_t_dt_2=subs(U,t,tt+dt/2);
DM_t_dt_2=subs(DM,t,tt+dt/2);

M_t_dt=subs(M,t,tt+dt);
U_t_dt=subs(U,t,tt+dt);
DM_t_dt=subs(DM,t,tt+dt);

k1=dtM_t \(U_t -(R+DM_t )(II(:,end) ));
k2=dtM_t_dt_2\(U_t_dt_2-(R+DM_t_dt_2)(II(:,end)+05k1));
k3=dtM_t_dt_2\(U_t_dt_2-(R+DM_t_dt_2)(II(:,end)+05k2));
k4=dtM_t_dt \(U_t_dt -(R+DM_t_dt )(II(:,end)+k3 ));

I_t=II(:,end)+(k1+2k2+2k3+k4)/6;
II=[II,I_t];
end
plot(tout',II')

假设求解初值问题：
y'=y-2x/y
(0<x<1)
y(0)=1
设步长h=02,从x=0直到x=1用四阶龙格库塔法：
private
sub
form_click()
dim
x
as
single,
y
as
single
dim
k1
as
single,
k2
as
single,
k3
as
single,
k4
as
single
dim
y1
as
single,
y2
as
single,
y3
as
single,
y4
as
single
x
=
0
y
=
1
h
=
02
for
i
=
0
to
1
step
02
k1
=
f(x,
y)
y1
=
y
+
h
/
2

k1
x
=
x
+
h
/
2
k2
=
f(x,
y1)
y2
=
y
+
h
/
2

k2
k3
=
f(x,
y2)
y3
=
y
+
h

k3
x
=
x
+
h
/
2
k4
=
f(x,
y3)
y
=
y
+
h
/
6

(k1
+
2

k2
+
2

k3
+
k4)
print
"x=";
x;
"
y=";
y
next
end
sub
private
function
f(a
as
single,
b
as
single)
as
single
f
=
b
-
2

a
/
b
end
function

>