柠檬小巫 发表于 2008-4-7 21:09

请教这个关于微分方程组的程序出错在哪

方程A(r)w=0
根据DET（A（r))=0,求得{ri} i=1:12
然后代入原方程，求得{wi=(w1i,w2i,w3i,w4i,w5i,w6i)}i=1:12
又u1(x)=c1*w1i*exp(ri*x)
u2(x)=c2*w2i*exp(ri*x)
p1(x)=c3*w3i*exp(ri*x)
p2(x)=c4*w4i*exp(ri*x)
v1(x)=c5*w5i*exp(ri*x)
v2(x)=c6*w6i*exp(ri*x)
边界条件为：
u1'(0)=0
u1'(L)=0
u2'(0)=0
u2'(L)=0
p1'(0)=0
p1'(L)=0
p2'(0)=0
p2'(L)=0
v1'(0)-p1(0)=0
v1'(L)-p1(L)=0
v2'(0)-p2(0)=0
v2'(L)-p2(L)=0


将w,r代入上述方程组,可得Mc=0
求DET（M）

具体程序如下：


m=2*pi*58.339
L=3.5
A1=3*10^-2
J1=9*10^-6
E1=4.539*10^10
b1=2600
G1=1.945*10^10
A2=1.64*10^-3
J2=5.41*10^-6
E2=2.1*10^11
b2=7850
G2=8.08*10^10
es=0.1
ec=0
EA=8.59*10^8
K=2.858*10^8
d=0.22
K1=1.2
K2=2.49
k=K/d
p=EA/d
b3=b1+b2
r=sym('r')


A=[E1*A1*r^2-k+b1*A1*m^2,k,0,k*es,k*ec*r/2,k*ec*r/2;
k,E2*A2*r^2-k+b2*A2*m^2,0,-k*es,-k*ec*r/2,-k*ec*r/2;
0,0,-r,0,r^2+b1*K1*m^2/G1-K1*p/(G1*A1),K1*p/(G1*A1);
0,0,0,-r,K2*p/(G2*A2),r^2+b2*K2*m^2/G2-K2*p/(G2*A2);
k*ec/2,-k*ec/2,E1*J1*r^2+b1*J1*m^2-G1*A1/K1,-k*ec*es/2,-k*ec^2*r/2+G1*A1*r/K1,-k*ec^2*r/6;
k*es+k*ec/2,-k*es-k*ec/2,0,E2*J2*r^2+b2*J2*m^2-G2*A2/K2-k*es^2,-k*ec*es*r/2-k*ec^2*r/6,-k*ec^2*r/3+G2*A2*r/K2-k*ec*ec*r/2]
B=det(A)
C=sym2poly(B)
E=roots(C)
Matrix=[]
for i=1:12
Temp=subs(A,r,E(i))
H=null(Temp)
Matrix=
end
w1=Matrix(1,:)
w2=Matrix(2,:)
w3=Matrix(3,:)
w4=Matrix(4,:)
w5=Matrix(5,:)
w6=Matrix(6,:)
syms c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 x
c='
for i=1:12
y1(i)=c(i)*w1(i)*exp(E(i)*x)
y2(i)=c(i)*w2(i)*exp(E(i)*x)
y3(i)=c(i)*w3(i)*exp(E(i)*x)
y4(i)=c(i)*w4(i)*exp(E(i)*x)
y5(i)=c(i)*w5(i)*exp(E(i)*x)
y6(i)=c(i)*w6(i)*exp(E(i)*x)
u1=sum(y1)
u2=sum(y2)
q1=sum(y3)
q2=sum(y4)
v1=sum(y5)
v2=sum(y6)
end
d1u1=vpa(diff(u1,x,1),5)
d1u2=vpa(diff(u2,x,1),5)
d1q1=vpa(diff(q1,x,1),5)
d1q2=vpa(diff(q2,x,1),5)
d1v1=vpa(diff(v1,x,1),5)
d1v2=vpa(diff(v2,x,1),5)
d1u10=vpa(subs(d1u1,x,0),5)
d1u1L=vpa(subs(d1u1,x,L),5)
d1u20=vpa(subs(d1u2,x,0),5)
d1u2L=vpa(subs(d1u2,x,L),5)
d1q10=vpa(subs(d1q1,x,0),5)
d1q1L=vpa(subs(d1q1,x,L),5)
d1q20=vpa(subs(d1q2,x,0),5)
d1q2L=vpa(subs(d1q2,x,L),5)
d1v10=vpa(subs(d1v1,x,0),5)
d1v1L=vpa(subs(d1v1,x,L),5)
d1v20=vpa(subs(d1v2,x,0),5)
d1v2L=vpa(subs(d1v2,x,L),5)
q10=vpa(subs(q1,x,0),5)
q1L=vpa(subs(q1,x,L),5)
q20=vpa(subs(q2,x,0),5)
q2L=vpa(subs(q2,x,L),5)
T10=d1v10-q10
T1L=d1v1L-q1L
T20=d1v20-q20
T2L=d1v2L-q2L
SV=

M1=[]
for i=1:12
B1=eye(12)
Tem1=subs(SV(1),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B1(i,:))
M1=
end
M2=[]
for i=1:12
B2=eye(12)
Tem2=subs(SV(2),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B2(i,:))
M2=
end

M3=[]
for i=1:12
B3=eye(12)
Tem3=subs(SV(3),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B3(i,:))
M3=
end
M4=[]
for i=1:12
B4=eye(12)
Tem4=subs(SV(4),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B4(i,:))
M4=
end
M5=[]
for i=1:12
B5=eye(12)
Tem5=subs(SV(5),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B1(i,:))
M5=
end
M6=[]
for i=1:12
B6=eye(12)
Tem6=subs(SV(6),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B6(i,:))
M6=
end
M7=[]
for i=1:12
B7=eye(12)
Tem7=subs(SV(7),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B7(i,:))
M7=
end
M8=[]
for i=1:12
B8=eye(12)
Tem8=subs(SV(8),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B8(i,:))
M8=
end
M9=[]
for i=1:12
B9=eye(12)
Tem9=subs(SV(9),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B9(i,:))
M9=
end
M10=[]
for i=1:12
B10=eye(12)
Tem10=subs(SV(10),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B10(i,:))
M10=
end
M11=[]
for i=1:12
B11=eye(12)
Tem11=subs(SV(11),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B11(i,:))
M11=
end
M12=[]
for i=1:12
B12=eye(12)
Tem12=subs(SV(12),{c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12},B12(i,:))
M12=
end

M=[M1
M2
M3
M4
M5
M6
M7
M8
M9
M10
M11
M12]
P=det(M)



偶计算的结果和书上的差很多，不知道问题出在哪?请高手指点～拜托了～～



[ 本帖最后由 eight 于 2008-4-7 21:12 编辑 ]

eight 发表于 2008-4-7 21:12

原帖由 柠檬小巫 于 2008-4-7 21:09 发表 http://www.chinavib.com/forum/images/common/back.gif
方程A(r)w=0
根据DET（A（r))=0,求得{ri} i=1:12
然后代入原方程，求得{wi=(w1i,w2i,w3i,w4i,w5i,w6i)}i=1:12
又u1(x)=c1*w1i*exp(ri*x)
u2(x)=c2*w2i*exp(ri*x)
p1(x)=c3*w3i*exp(ri*x)
p2(x)=c4*w4i*exp(ri* ... 新手若然要得到及时的、有用的帮助，请逐个置顶帖仔细阅读一次

柠檬小巫 发表于 2008-4-7 21:32

多谢提醒～～

谢谢了！
页: [1]
查看完整版本: 请教这个关于微分方程组的程序出错在哪