题目
用四阶龙格库塔法求解矩阵微分方程
要求电流就是求解矩阵微分方程:(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;
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个电流与时间的关系
要求电流就是求解矩阵微分方程:(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;
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个电流与时间的关系
提问时间:2020-12-14
答案
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/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;
0;
0
];
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;
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=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')
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/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;
0;
0
];
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;
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=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')
举一反三
已知函数f(x)=x,g(x)=alnx,a∈R.若曲线y=f(x)与曲线y=g(x)相交,且在交点处有相同的切线,求a的值和该切线方程.
我想写一篇关于奥巴马的演讲的文章,写哪一篇好呢?为什么好
奥巴马演讲不用看稿子.为什么中国领导演讲要看?
想找英语初三上学期的首字母填空练习……
英语翻译
最新试题
- 1一道英文题:You can ’t work all the time.You need a h_____for your free time.(填补空格单词)
- 2一项工程,甲乙合作8天完成,如果甲先独做6天,然后乙再做9天完成任务,乙需独做多少天完成任务?
- 32011年12月20日(星期二)怎么翻译成英语
- 4某项工程,甲乙两队合做20天完成,甲队单独做30天完成.现在两队合做15天后,余下的由乙队完成,还要多少?
- 5土字旁加一个劳字读什么
- 6已知函数f(x)=1/2cos^2x+√3/2sinxcosx+1,x属于R,(1)求函数f(x)的最小正周期
- 7生活中常见的固体,液体,气体热胀冷缩的例子 并分析对人们生活造成的影 响
- 8英语翻译
- 9缩醛磷脂的功能是什么?
- 10vb问题利用随机函数产生80个10~99之间的随机整数,
热门考点