|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我是matlab的新手,前几天要编个程序需要用matlab先算出初值。不过结果始终不满意,我先简单介绍下问题,再给出自己程序,希望高手们帮我看看。谢谢了!
位移u(数组),然后是中间量a,b,c (也都是数组)是u 的函数。根据这些算出应力q1,q2,q3..再是根据应力算出力和力偶n,m并得到加速度。最后是用a,b,c,m,n得到下一时刻位移u.
近似后,u = (1/4 )*a*t^2,作为方程。
这是梁的位移,从u[21]到u[41]是挠度,应该对称。可以用来检验结果是否正确。
我是新手,可能是最基本的地方错了,具体的某个量应该没什么问题,如果不理解可以不用考虑。
程序:
function work
format long
u0 = zeros(42,1);
u = zeros(42,1);
a = zeros(21,1);
b = zeros(21,1) ;
c = zeros(21,1);
q1 = zeros(21,5);
p2 = zeros(21,5);
q3 = zeros(21,1);
n = zeros(21,1);
m = zeros(21,1);
u = fsolve(@myfun,u0)
function f = myfun(u)
u(1) = 0; u(21) = 0; u(22) = 0; u(42) = 0;
a(1) = u(2)/0.01;
a(21) = -u(20)/0.01;
for i=2:20
a(i)=(u(i+1)-u(i-1))/(2*0.01);
end
b(1) = u(23)/0.01;
b(21) = -u(41)/0.01;
for i=2:20
b(i)=(u(i+22)-u(i+20))/(2*0.01);
end
c(1) = (u(24)-2*u(23))/(0.01)^2;
c(21) = (u(40)-2*u(41))/(0.01)^2;
for i=2:20
c(i) = (u(i+22)-2*u(i+21)+u(i+20))/(0.01)^2;
end
for i=1:21
for j=1:5
if abs(80e9*(a(i)-j*(4e-4)*c(i)+0.5*b(i)^2))<=0.3e9
q1(i,j)=80e9*(a(i)-j*(4e-4)*c(i)+0.5*b(i)^2);
elseif 80e9*(a(i)-j*(4e-4)*c(i)+0.5*b(i)^2)>0.3e9
q1(i,j)=0.3e9;
else
q1(i,j)=-0.3e9;
end
end
end
for i=1:21
for j=1:5
if abs(80e9*(a(i)+j*(4e-4)*c(i)+0.5*b(i)^2))<=0.3e9
p2(i,j)=80e9*(a(i)+j*(4e-4)*c(i)+0.5*b(i)^2);
elseif 80e9*(a(i)+j*(4e-4)*c(i)+0.5*b(i)^2)>0.3e9
p2(i,j)=0.3e9;
else
p2(i,j)=-0.3e9;
end
end
end
for i=1:21
if abs(80e9*(a(i)+0.5*b(i)^2))<=0.3e9
q3(i)=80e9*(a(i)+0.5*b(i)^2);
elseif 80e9*(a(i)+0.5*b(i)^2)>0.3e9
q3(i)=0.3e9;
else
q3(i)=-0.3e9;
end
end
for i=1:21
n(i) = (q1(i,5) + p2(i,5))*0.02*(4e-4)/2+q3(i)*0.02*(4e-4);
for j=1:4
n(i) = n(i) + (q1(i,j)+p2(i,j))*0.02*(4e-4);
end
end
for i = 1:21
m(i) = -(q1(i,5)-p2(i,5))*0.02*5*(4e-4)^2/2;
for j=1:4
m(i) = m(i) -(q1(i,j)-p2(i,j))*j*0.02*(4e-4)^2;
end
end
for i=2:20
f(i) = (1/4)*(((1+a(i+1))*n(i+1) - (1+a(i-1))*n(i-1) + m(i+1)*c(i+1) - m(i-1)*c(i-1))/(0.01*2*0.02*(4e-3)*2700) +b(i)*96500/((4e-3)*2700))*(5e-4)^2 -u(i);
end
for i=23:41
f(i) = (1/4)*(-(m(i-20) - 2*m(i-21) + m(i-22)/(0.02*(4e-3)*2700*(0.01)^2)) + 96500/((4e-3)*2700))*(5e-4)^2 - u(i);
end |
|