找回密码
 注册
Simdroid-非首页
查看: 76|回复: 0

解微分方程组

[复制链接]
发表于 2012-4-29 04:57:05 | 显示全部楼层 |阅读模式 来自 美国
本帖最后由 wintermei 于 2012-4-29 06:19 编辑

我要解微分方程的是
d V/d t=AV + VA^{T}+ D,这里 A, D 是已知6*6的矩阵,

A=[-0.5*K 0 0 0 0 G1;0 -0.5*K 0 0 -G1 0;0 0 -0.5*K 0 0 -G2;0 0 0 -0.5*K -G2 0;0 G1 0 -G2 -0.5*RM 0;-G1 0 -G2 0 0 -0.5*RM]
D=[0.5*K 0 0 0 0 0;0 0.5*K 0 0 0 0;0 0 0.5*K 0 0 0;0 0 0 0.5*K 0 0;0 0 0 0 0 0;0 0 0 0 0 0.5*RM*(2*N+1)];
这里 K,G1,G2,RM,N 是已知的.



我的方法是让
V= [V(1) V(2) V(3) V(4) V(5) V(6);
V(7) V(8) V(9) V(10) V(11) V(12);
V(13) V(14) V(15) V(16) V(17) V(18);
V(19) V(20) V(21) V(22) V(23) V(24);
V(25) V(26) V(27) V(28) V(29) V(30);
V(31) V(32) V(33) V(34) V(35) V(36)]
然后将V代入 d V/d t=AV + VA^{T}+D 得到36 个方程如下:

dV(1)/dt =K/2 - (K*V(1))/2 + V(6)*conj(G1) - (V(1)*conj(K))/2 + G1*V(31);
dV(2)/dt =G1*V(32) - V(5)*conj(G1) - (V(2)*conj(K))/2 - (K*V(2))/2;
dV(3)/dt =G1*V(33) - V(6)*conj(G2) - (V(3)*conj(K))/2 - (K*V(3))/2;
dV(4)/dt =G1*V(34) - V(5)*conj(G2) - (V(4)*conj(K))/2 - (K*V(4))/2;
dV(5)/dt =V(2)*conj(G1) - (K*V(5))/2 - V(4)*conj(G2) - (V(5)*conj(RM))/2 + G1*V(35);
dV(6)/dt =G1*V(36) - V(1)*conj(G1) - V(3)*conj(G2) - (V(6)*conj(RM))/2 - (K*V(6))/2;
dV(7)/dt = V(12)*conj(G1) - (K*V(7))/2 - (V(7)*conj(K))/2 - G1*V(25);
dV(8)/dt = K/2 - (K*V(8))/2 - V(11)*conj(G1) - (V(8)*conj(K))/2 - G1*V(26);
dV(9)/dt =- (K*V(9))/2 - V(12)*conj(G2) - (V(9)*conj(K))/2 - G1*V(27);
dV(10)/dt = - (K*V(10))/2 - V(11)*conj(G2) - (V(10)*conj(K))/2 - G1*V(28);
dV(11)/dt =V(8)*conj(G1) - (K*V(11))/2 - V(10)*conj(G2) - (V(11)*conj(RM))/2 - G1*V(29);
dV(12)/dt =- (K*V(12))/2 - V(7)*conj(G1) - V(9)*conj(G2) - (V(12)*conj(RM))/2 - G1*V(30);
dV(13)/dt =V(18)*conj(G1) - (K*V(13))/2 - (V(13)*conj(K))/2 - G2*V(31);
dV(14)/dt =- (K*V(14))/2 - V(17)*conj(G1) - (V(14)*conj(K))/2 - G2*V(32);
dV(15)/dt =K/2 - (K*V(15))/2 - V(18)*conj(G2) - (V(15)*conj(K))/2 - G2*V(33);
dV(16)/dt =- (K*V(16))/2 - V(17)*conj(G2) - (V(16)*conj(K))/2 - G2*V(34);
dV(17)/dt = V(14)*conj(G1) - (K*V(17))/2 - V(16)*conj(G2) - (V(17)*conj(RM))/2 - G2*V(35);
dV(18)/dt=- (K*V(18))/2 - V(13)*conj(G1) - V(15)*conj(G2) - (V(18)*conj(RM))/2 - G2*V(36);
dV(19)/dt=V(24)*conj(G1) - (K*V(19))/2 - (V(19)*conj(K))/2 - G2*V(25);
dV(20)/dt= - (K*V(20))/2 - V(23)*conj(G1) - (V(20)*conj(K))/2 - G2*V(26);
dV(21)/dt=- (K*V(21))/2 - V(24)*conj(G2) - (V(21)*conj(K))/2 - G2*V(27);
dV(22)/dt=K/2 - (K*V(22))/2 - V(23)*conj(G2) - (V(22)*conj(K))/2 - G2*V(28);
dV(23)/dt=V(20)*conj(G1) - (K*V(23))/2 - V(22)*conj(G2) - (V(23)*conj(RM))/2 - G2*V(29);
dV(24)/dt =- (K*V(24))/2 - V(19)*conj(G1) - V(21)*conj(G2) - (V(24)*conj(RM))/2 - G2*V(30);
dV(25)/dt=V(30)*conj(G1) - (RM*V(25))/2 - (V(25)*conj(K))/2 + G1*V(7) - G2*V(19);
dV(26)/dt= G1*V(8) - V(29)*conj(G1) - (V(26)*conj(K))/2 - (RM*V(26))/2 - G2*V(20);
dV(27)/dt=G1*V(9) - V(30)*conj(G2) - (V(27)*conj(K))/2 - (RM*V(27))/2 - G2*V(21);
dV(28)/dt=G1*V(10) - V(29)*conj(G2) - (V(28)*conj(K))/2 - (RM*V(28))/2 - G2*V(22);
dV(29)/dt= V(26)*conj(G1) - (RM*V(29))/2 - V(28)*conj(G2) - (V(29)*conj(RM))/2 + G1*V(11) - G2*V(23);
dV(30)/dt=G1*V(12) - V(25)*conj(G1) - V(27)*conj(G2) - (V(30)*conj(RM))/2 - (RM*V(30))/2 - G2*V(24);
dV(31)/dt=V(36)*conj(G1) - (RM*V(31))/2 - (V(31)*conj(K))/2 - G1*V(1) - G2*V(13);
dV(32)/dt= - (RM*V(32))/2 - V(35)*conj(G1) - (V(32)*conj(K))/2 - G1*V(2) - G2*V(14);
dV(33)/dt=- (RM*V(33))/2 - V(36)*conj(G2) - (V(33)*conj(K))/2 - G1*V(3) - G2*V(15);
dV(34)/dt=- (RM*V(34))/2 - V(35)*conj(G2) - (V(34)*conj(K))/2 - G1*V(4) - G2*V(16);
dV(35)/dt =V(32)*conj(G1) - (RM*V(35))/2 - V(34)*conj(G2) - (V(35)*conj(RM))/2 - G1*V(5) - G2*V(17);
dV(36)/dt=(RM*(2*N + 1))/2 - V(31)*conj(G1) - V(33)*conj(G2) - (V(36)*conj(RM))/2 - (RM*V (36))/2 - G1*V(6) - G2*V(18);
然后用ODE45命令解这36个方程.

我的方法看起来很笨,谁有没有简单的方法?谢谢!


您需要登录后才可以回帖 登录 | 注册

本版积分规则

Simapps系列直播

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-10-2 02:53 , Processed in 0.026321 second(s), 10 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表