矩阵龙格库塔
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
要求电流就是求解矩阵微分方程：（R+pM(t)）*I(t)+M(t)*pI(t)-U(t)=0,
其中p是求导，
R是6*6常数矩阵，
M(t)是6*6的时变矩阵，
U(t)是6*1的时变矩阵，
求I(t)，也是6*1的矩阵。
已知条件：
M=[0,0,0, -15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi, 15 727/10000*sin(1/12*pi+80*pi*t)*pi;0,0,0,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/100 00*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi;0,0,0,15727/10000*cos(1/4 *pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*p i;-15727/10000*sin(5/12*pi+80*pi*t), 15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/100 00*cos(1/4*pi+80*pi*t)*pi, 0,
0, 0;
15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi, 15727/1 0000*sin(1/12*pi+80*pi*t)*pi, 0,
0, 0;
15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/ 10000*sin(5/12*pi+80*pi*t)*pi, 0,
0, 0];
U=[380*cos(100*pi*t);
380*cos(100*pi*t-2*pi/3);
380*cos(100*pi*t+2*pi/3);
0;
0;
0];
R=[0.0247,0,0,0,0,0;
0,0.0247,0,0,0,0;
0,0,0.0247,0,0,0;
0,0,0,0.0193,-0.0193,0;
0,0,0,0,0.0193,-0.0193;
0,0,0,1,1,1];
电流I的初值是：I(0)=[0;0;0;0;0;0];
t是时间
预期最后的结果是6个电流与时间的关系
谢谢
分享到：
举报| 2010-06-22 16:47提问者采纳
global R M U
syms t
R=[
0.0247,0,0,0,0,0;
0,0.0247,0,0,0,0;
0,0,0.0247,0,0,0;
0,0,0,0.0193,-0.0193,0;
0,0,0,0,0.0193,-0.0193;
0,0,0,1,1,1
];
M=[
0,0,0,-15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi,15727/10 000*sin(1/12*pi+80*pi*t)*pi;
0,0,0,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi;
0,0,0, 15727/10000*cos(1/4*pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi;
-15727/10000*sin(5/12*pi+80*pi*t),15727/10000*sin(1/12*pi+80*pi*t)*pi,
15727/10000*cos(1/4*pi+80*pi*t)*pi,0,0,0;
15727/10000*cos(1/4*pi+80*pi*t)*pi,
-15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi,0,0,0;
15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi,0,0,0
];
U=[
380*cos(100*pi*t);
380*cos(100*pi*t-2*pi/3);
380*cos(100*pi*t+2*pi/3);
0;
0;
];
DM=diff(M,t); %就是M对t求导
%%%要调节的参数在这里
%%注意，你的M有点奇异，计算很快发散掉了，你检察一下相关的参数吧。
%%det(subs(M,t,0))
I0=[0;0;0;0;0;0];
tstart=0;
tend=0.1;
dt=0.1;
%%%end
tout=tstart:dt:tend;
n=length(tout);
M_t_dt=subs(M,t,tstart); %P=subs(P,'t',x)就是把P表达式中所有't',都用具体的x值代替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); %（R+pM(t)）*I(t)+M(t)*pI(t)-U(t)=0
k1=dt*M_t \(U_t -(R+DM_t )*(II(:,end) ));
k2=dt*M_t_dt_2\(U_t_dt_2-(R+DM_t_dt_2)*(II(:,end)+0.5*k1));
k3=dt*M_t_dt_2\(U_t_dt_2-(R+DM_t_dt_2)*(II(:,end)+0.5*k2));
k4=dt*M_t_dt \(U_t_dt -(R+DM_t_dt )*(II(:,end)+k3 ));
I_t=II(:,end)+(k1+2*k2+2*k3+k4)/6; II=[II,I_t];
end
plot(tout',II')。