zhdg1989 发表于 2013-11-4 11:07

使用fsolve出现问题,请大家帮我看看,谢谢!

用fsolve解下边的非线性方程组,调节不同的参数得到的解相同,哪里出错了呢?
function shiyanjidi
clc
syms rho A I E L C1 M K2 l0 C2 m2 l v
rho=7200;
A=4.8;
I=2.5498;
E=5*10^10;
L=40;
C1=0.02;
M=2000;
K2=2392;
l0=6.72;
C2=200;
m2=2390;
l=1;
v=1/2*sqrt(E*I/rho/A)*pi/L;
syms beta f omega omega1 phi gamma1 gamma2 gamma3 mu1 mu2 alpha
beta=A*l^2/8/I
f=2*M*9.8*L^3/(l*E*I*pi^4)
phi=2*M/(rho*A*L)
omega1=sqrt(E*I*pi^4/(rho*A*L^4))
gamma1=C1/(rho*A*omega1)
gamma2=C2/(omega1*m2)
gamma3=K2*rho*A*L^4/(m2*E*I*pi^4)
mu1=2*m2/(rho*A*L)
mu2=1
alpha=l0/l
omega=sqrt(rho*A/E/I)*L*v/pi
p=fsolve('myfun',,optimset('Display','off'))
function q=myfun(p)
syms R1 R2;
=ellipke(sqrt(p(4)^2/(alpha^2+(p(4))^2)));
R1=2*gamma3*p(4)+8*gamma3*alpha^2*K/pi/p(4)/sqrt(alpha^2+(p(4))^2)-8*gamma3*sqrt(alpha^2+(p(4))^2)*E/pi/p(4);
R2=-2*gamma2*omega*p(4);
q(1)=p(1)+beta*(p(1))^3+3/2*beta*(p(1))*(p(2))^2+3/2*beta*(p(1))*(p(3))^2+omega^2*phi*p(2)-2*f/pi;
q(2)=3*beta*(p(1))^2*p(2)+3/4*beta*(p(2))^3+3/4*beta*p(2)*(p(3))^2+(1-4*omega^2-2*omega^2*phi)*p(2)+2*gamma1*omega*p(3)+mu1*((R1)*cos(p(5))+(R2)*sin(p(5)))+4*f/3/pi;
q(3)=3*beta*(p(1))^2*p(3)+3/4*beta*(p(3))^3+3/4*beta*(p(2))^2*p(3)+(1-4*omega^2-2*omega^2*phi)*p(3)-2*gamma1*omega*p(2)+mu1*((R2)*cos(p(5))-(R1)*sin(p(5)));
q(4)=-4*omega^2*p(4)*cos(p(5))+(R1)*cos(p(5))+(R2)*sin(p(5))+4*omega^2*p(2);
q(5)=4*omega^2*p(4)*sin(p(5))+(R2)*cos(p(5))-(R1)*sin(p(5))+4*omega^2*p(3);
end
end

ChaChing 发表于 2013-11-6 14:31

1.求助完整格式:出错代码和出错提示
2.建议楼主稍微学下函数使用
3.建议数值与符号别混用

ChaChing 发表于 2013-11-7 11:55

function shiyanjidi
%syms rho A I E L C1 M K2 l0 C2 m2 l v
%syms beta f omega omega1 phi gamma1 gamma2 gamma3 mu1 mu2 alpha
p=fsolve(@myfun,,optimset('Display','off'))
end

function q=myfun(p)
%syms R1 R2;
rho=7200; A=4.8; I=2.5498; E=5*10^10; L=40; C1=0.02; M=2000; K2=2392; l0=6.72; C2=200; m2=2390; l=1;
v=1/2*sqrt(E*I/rho/A)*pi/L; beta=A*l^2/8/I; f=2*M*9.8*L^3/(l*E*I*pi^4);
phi=2*M/(rho*A*L); omega1=sqrt(E*I*pi^4/(rho*A*L^4));
gamma1=C1/(rho*A*omega1); gamma2=C2/(omega1*m2); gamma3=K2*rho*A*L^4/(m2*E*I*pi^4);
mu1=2*m2/(rho*A*L); mu2=1; alpha=l0/l; omega=sqrt(rho*A/E/I)*L*v/pi;
=ellipke(sqrt(p(4)^2/(alpha^2+(p(4))^2)));
R1=2*gamma3*p(4)+8*gamma3*alpha^2*K/pi/p(4)/sqrt(alpha^2+(p(4))^2)-8*gamma3*sqrt(alpha^2+(p(4))^2)*E/pi/p(4);
R2=-2*gamma2*omega*p(4);
q(1)=p(1)+beta*(p(1))^3+3/2*beta*(p(1))*(p(2))^2+3/2*beta*(p(1))*(p(3))^2+omega^2*phi*p(2)-2*f/pi;
q(2)=3*beta*(p(1))^2*p(2)+3/4*beta*(p(2))^3+3/4*beta*p(2)*(p(3))^2+(1-4*omega^2-2*omega^2*phi)*p(2)+2*gamma1*omega*p(3)+mu1*((R1)*cos(p(5))+(R2)*sin(p(5)))+4*f/3/pi;
q(3)=3*beta*(p(1))^2*p(3)+3/4*beta*(p(3))^3+3/4*beta*(p(2))^2*p(3)+(1-4*omega^2-2*omega^2*phi)*p(3)-2*gamma1*omega*p(2)+mu1*((R2)*cos(p(5))-(R1)*sin(p(5)));
q(4)=-4*omega^2*p(4)*cos(p(5))+(R1)*cos(p(5))+(R2)*sin(p(5))+4*omega^2*p(2);
q(5)=4*omega^2*p(4)*sin(p(5))+(R2)*cos(p(5))-(R1)*sin(p(5))+4*omega^2*p(3);
endp =

    0.0002   -0.0676    0.0467    0.0974   -2.5290
页: [1]
查看完整版本: 使用fsolve出现问题,请大家帮我看看,谢谢!