用四阶龙格库塔法求解矩阵微分方程要求电流就是求解矩阵微分方程-查字典问答网

来自胡德安的问题

　　用四阶龙格库塔法求解矩阵微分方程要求电流就是求解矩阵微分方程：（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的矩阵.已

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

　　要求电流就是求解矩阵微分方程：（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,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*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];

　　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个电流与时间的关系

1回答

2020-10-21 17:33

我要回答

请先登录

万军

　　globalRMU

　　symst

　　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/10000*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);

　　];

　　DM=diff(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);

　　U_t_dt=subs(U,t,tstart);

　　DM_t_dt=subs(DM,t,tstart);

　　II=I0;

　　fori=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=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')

2020-10-21 17:37:41

最新问答

推荐文章

猜你喜欢

附近的人在看

推荐阅读

拓展阅读