[求助]基于MATLAB的高速储能飞轮强度设计
 

本人在做毕业设计，要用到MATLAB，但现在有段程序，总是得不到理想的结果，求高手帮忙修改下，谢谢！
syms r dr v vx E Et Er A t F Beta1 lam result_1 x delt P1 ii
dr=0.001;
ri=0.112;
ro=0.145;
v=0.3;
vx=0.28;
E=202*10^9;
Et=38.6*10^9;
Er=8.27*10^9;
A=3.14*dr^2;
t=2*dr;
F=3; %Beta1=112/115;
lam=sqrt(Et/Er);
for ii=0.115:0.001:0.145
Beta1=0.112/ii;
eq=x*(ii-x)*A-(1-v)/E*F*t*ii-F/Et*t*(ii-x)*((1+Beta1^(2*lam))/(1-Beta1^(2*lam))*lam+vx);
result_1=solve(eq,x);
s=vpa(result_1,6);
delt=s(1); %s=double(result_1)
%ff=subs(result_1,{ii},{115/1000:1/1000:145/1000})
%
P1=delt/((1-v)/E*ii+(ii-delt)/Et*((1+ Beta1^(2*lam))/(1-Beta1^(2*lam))*lam+vx));
%
% %内环 P=0 Q=P1%

huan_1=(0*Beta1^(lam+1)-P1)/(1-Beta1^(2*lam))*lam*(ii/ro)^(lam-1)+(0-P1*Beta1^(lam-1))/(1-Beta1^(2*lam))*lam*Beta1^(lam+1)*(ii/ro)^(-lam-1);
huan_1a=double(huan_1)
%
jing_1=(0*Beta1^(lam+1)-P1)/(1-Beta1^(2*lam))*(ii/ro)^(lam-1)+(0-P1*Beta1^(lam-1))/(1-Beta1^(2*lam))*Beta1^(lam+1)*(ii/ro)^(-lam-1);
jing_1a=double(jing_1)
%外环 P=P1 Q=0
huan_2=(P1*Beta1^(lam+1)-0)/(1-Beta1^(2*lam))*lam*(ii/ro)^(lam-1)+(P1-0*Beta1^(lam-1))/(1-Beta1^(2*lam))*lam*Beta1^(lam+1)*(ii/ro)^(-lam-1);
huan_2a=double(huan_2)
jing_2=(P1*Beta1^(lam+1)-0)/(1-Beta1^(2*lam))*(ii/ro)^(lam-1)+(P1-0*Beta1^(lam-1))/(1-Beta1^(2*lam))* Beta1^(lam+1)*(ii/ro)^(-lam-1);
jing_2a=double(jing_2)

%jing=jing_1

%huan=huan_li+huan_n;
%jing=jing_li+jing_n;


plot(ii,huan_1a, '.k')
plot(ii,jing_1a, 'ob')
plot(ii,huan_2a, '^r')
plot(ii,jing_2a, '*c')
hold on
grid on
xlabel('r')
ylabel('yingli')
legend('huan_la','jing_1a','huan_2a','jing_2a')
end