24小时热门版块排行榜

返回列表

有累人1

金虫 (初入文坛)

应助: 0 (幼儿园)
金币: 744
帖子: 5
在线: 16.9小时
虫号: 3385453
注册: 2014-08-27
性别: GG
专业: 机构学与机器人

[求助] MATLA B求解微分方程

这个程序得不到想要的结果，运行可能导致电脑死机，求大神帮助，TXT在附件中。

clc;clear all;

V=textread('D:\Program Files (x86)\MATLAB71\work\1.txt');
i=1;
c=1;
re=1;
U=[V(i,1) V(i,2) V(i,3);V(i,4) V(i,5) V(i,6);V(i,7) V(i,8) V(i,9)];
r=(U+U')/2;
w=(U-U')/2;
Dr=(0.0184*exp(-0.7148*c*re))*(sqrt(2*sum(diag(r*r))));
A=[1/3 0 0;0 1/3 0;0 0 1/3];

syms a11 a12 a13 a21 a22 a23 a31 a32 a33;
f1=eval((w(1,1)*a11+w(1,2)*a21+w(1,3)*a31-w(1,1)*a11-w(2,1)*a12-w(3,1)*a13)/(-2)+(r(1,1)*a11+r(1,2)*a21+r(1,3)*a31+r(1,1)*a11+r(2,1)*a12+r(3,1)*a13-2*r(1,1)*(6*a11/7-3/35)-2*r(1,2)*3*a12/7-2*r(1,3)*3*a13/7-2*r(2,1)*(2*a12+a21)/7-2*r(2,2)*((a11+a22)/7-1/35)-2*r(2,3)*a23/7-2*r(3,1)*(2*a13+a31)/7-2*r(3,2)*a32/7-2*r(3,3)*((a11+a33)/7-1/35))/2+Dr*(1-3*a11));
f2=eval((w(1,1)*a12+w(1,2)*a22+w(1,3)*a32-w(1,2)*a11-w(2,2)*a12-w(3,2)*a13)/(-2)+(r(1,1)*a12+r(1,2)*a22+r(1,3)*a32+r(1,2)*a11+r(2,2)*a12+r(3,2)*a13-2*r(1,1)*(a12+2*a21)/7-2*r(1,2)*((a11+a22)/7-1/35)-2*r(1,3)*a23/7-2*r(2,1)*((a11+a22)/7-1/35)-2*r(2,2)*3*a12/7-2*r(2,3)*a13/7-2*r(3,1)*a23/7-2*r(3,2)*a13/7-2*r(3,3)*a12/7)/2-3*Dr*a12);
f3=eval((w(1,1)*a13+w(1,2)*a23+w(1,3)*a33-w(1,3)*a11-w(2,3)*a12-w(3,3)*a13)/(-2)+(r(1,1)*a13+r(1,2)*a23+r(1,3)*a33+r(1,3)*a11+r(2,3)*a12+r(3,3)*a13-2*r(1,1)*(a13+2*a31)/7-2*r(1,2)*a32/7-2*r(1,3)*((a11+a33)/7-1/35)-2*r(2,1)*a32/7-2*r(2,2)*a13/7-2*r(2,3)*a12/7-2*r(3,1)*((a11+ a33)/7-1/35)-2*r(3,2)*a12/7-2*r(3,3)*3*a13/7)/2-3*Dr*a13);
f4=eval((w(2,1)*a11+w(2,2)*a21+w(2,3)*a31-w(1,1)*a21-w(2,1)*a22-w(3,1)*a23)/(-2)+(r(2,1)*a11+r(2,2)*a21+r(2,3)*a31+r(1,1)*a21+r(2,1)*a22+r(3,1)*a23-2*r(1,1)*3*a21/7-2*r(1,2)*((a11+a22)/7-1/35)-2*r(1,3)*a23/7-2*r(2,1)*((a11+a22)/7-1/35)-2*r(2,2)*(a21+2*a12)/7-2*r(2,3)*a13/7-2*r(3,1)*a23/7-2*r(3,2)*a13/7-2*r(3,3)*a21/7)/2-3*Dr*a21);
f5=eval((w(2,1)*a12+w(2,2)*a22+w(2,3)*a32-w(1,2)*a21-w(2,2)*a22-w(3,2)*a23)/(-2)+(r(2,1)*a12+r(2,2)*a22+r(2,3)*a32+r(1,2)*a21+r(2,2)*a22+r(3,2)*a23-2*r(1,1)*((a11+a22)/7-1/35)-2*r(1,2)*(2*a21+a12)/7-2*r(1,3)*a13/7-2*r(2,1)*3*a21/7-2*r(2,2)*(6*a22/7-3/35)-2*r(2,3)*3*a23/7-2*r(3,1)*a31/7-2*r(3,2)*(2*a23+a32)/7-2*r(3,3)*((a22+a33)/7-1/35))/2+Dr*(1-3*a22));
f6=eval((w(2,1)*a13+w(2,2)*a23+w(2,3)*a33-w(1,3)*a21-w(2,3)*a22-w(3,3)*a23)/(-2)+(r(2,1)*a13+r(2,2)*a23+r(2,3)*a33+r(1,3)*a21+r(2,3)*a22+r(3,3)*a23-2*r(1,1)*a23/7-2*r(1,2)*a31/7-2*r(1,3)*a21/7-2*r(2,1)*a31/7-2*r(2,2)*(a23+2*a32)/7-2*r(2,3)*((a22+a33)/7-1/35)-2*r(3,1)*a21/7-2*r(3,2)*((a22+a33)/7-1/35)-2*r(3,3)*3*a23/7)/2-3*Dr*a22);
f7=eval((w(3,1)*a11+w(3,2)*a21+w(3,3)*a31-w(1,1)*a31-w(2,1)*a32-w(3,1)*a33)/(-2)+(r(3,1)*a11+r(3,2)*a21+r(3,3)*a31+r(1,1)*a31+r(2,1)*a32+r(3,1)*a33-2*r(1,1)*3*a31/7-2*r(1,2)*a32/7-2*r(1,3)*((a11+a33)/7-1/35)-2*r(2,1)*a32/7-2*r(2,2)*a31/7-2*r(2,3)*a12/7-2*r(3,1)*((a11+a33)/7-1/35)-2*r(3,2)*a12/7-2*r(3,3)*(a31+2*a13)/7)/2-3*Dr*a31);
f8=eval((w(3,1)*a12+w(3,2)*a22+w(3,3)*a32-w(1,2)*a31-w(2,2)*a32-w(3,2)*a33)/(-2)+(r(3,1)*a12+r(3,2)*a22+r(3,3)*a32+r(1,2)*a31+r(2,2)*a32+r(3,2)*a33-2*r(1,1)*a32/7-2*r(1,2)*a31/7-2*r(1,3)*a21/7-2*r(2,1)*a31/7-2*r(2,2)*3*a32/7-2*r(2,3)*((a22+a33)/7-1/35)-2*r(3,1)*a21/7-2*r(3,2)*((a22+a33)/7-1/35)-2*r(3,3)*(a32+2*a23)/7)/2-3*Dr*a32);
f9=eval((w(3,1)*a13+w(3,2)*a23+w(3,3)*a33-w(1,3)*a31-w(2,3)*a32-w(3,3)*a33)/(-2)+(r(3,1)*a13+r(3,2)*a23+r(3,3)*a33+r(1,3)*a31+r(2,3)*a32+r(3,3)*a33-2*r(1,1)*((a11+a33)/7-1/35)-2*r(1,2)*a12/7-2*r(1,3)*(2*a31+a13)/7-2*r(2,1)*a21/7-2*r(2,2)*((a22+a33)/7-1/35)-2*r(2,3)*(2*a32+a23)/7-2*r(3,1)*3*a31/7-2*r(3,2)*3*a32/7-2*r(3,3)*(6*a33/7-3/35))/2+Dr*(1-3*a33));
[a11,a12,a13] = dsolve(strcat('Da11 =',char(f1)),strcat('Da12 =',char(f2)),strcat('Da13 =',char(f3)),'a11(0) = num2str(A(1,1))','a12(0) = num2str(A(1,2))','a13(0) = num2str(A(1,3))','t')
[a21,a22,a23] = dsolve(strcat('Da21 =',char(f4)),strcat('Da22 =',char(f5)),strcat('Da23 =',char(f6)),'a21(0) = num2str(A(2,1))','a22(0)= num2str(A(2,2))','a23(0) = num2str(A(2,3))','t')
[a31,a32,a33] = dsolve(strcat('Da31 =',char(f7)),strcat('Da32 =',char(f8)),strcat('Da33 =',char(f9)),'a31(0)=num2str(A(3,1))','a32(0)=num2str(A(3,2))','a33(0)=num2str(A(3,3))','t')

[a11,a12,a13,a21,a22,a23,a31,a32,a33] = dsolve(strcat('Da11 =',char(f1)),strcat('Da12 =',char(f2)),strcat('Da13 =',char(f3)),strcat('Da21 =',char(f4)),strcat('Da22 =',char(f5)),strcat('Da23 =',char(f6)),strcat('Da31 =',char(f7)),strcat('Da32 =',char(f8)),strcat('Da33 =',char(f9)),'a11(0)=num2str(A(1,1))','a12(0)=num2str(A(1,2))','a13(0)=num2str(A(1,3))','a21(0)=num2str(A(2,1))','a22(0)=num2str(A(2,2))','a23(0)=num2str(A(2,3))','a31(0)=num2str(A(3,1))','a32(0)=num2str(A(3,2))','a33(0)=num2str(A(3,3))','t')

t=0.01325;
A=[a11 a12 a13;a21 a22 a23;a31 a32 a33]

回复此楼