?? guaidian1.m
字號:
function [tm,xm]=guaidian1(dt,k1,m,eta,xg,x1,xs);%第一類拐點處理;
%%%%精確拐點非迭代法;紐馬克法求解;
% gama=0.505;
% beta=0.255;
gama=0.5;
beta=0.25;
w=sqrt(k1/m);
cc=2*eta*w*m;
dxs=xs-x1(1);
%采用精確拐點非迭代法求出屈服拐點;
pa=cc*(gama/2/beta-1)*x1(3)-m*(xg(2)-xg(1))/dt;
pb=m/2/beta*x1(3)+gama*cc/beta*x1(2)-k1*dxs;
pc=m*x1(2)/beta-gama*cc*dxs/beta;
pd=-m/beta*dxs;
if (pd==0)
tm=0;
else
if pa==0 %%%解二次方程;
a=pc/pb;b=pd/pb;
%%%%%%x=solve('x^2+a*x+b=0');
det(1)=-1/2*a+1/2*(a^2-4*b)^(1/2);
det(2)=-1/2*a-1/2*(a^2-4*b)^(1/2);
tm=dt;
ndet=0;
for ni=1:2
if det(ni)<dt&det(ni)>0
tm=det(ni);
ndet=ndet+1;
end
end
if(ndet>1|tm>=dt|tm<=0)
cuowu1='計算機精度所限,采用迭代處理為妥';
%牛頓法求解拐點時間;
t0=dt/2;
t1=t0-(t0^2+a*t0+b)/(2*t0+a);
while abs((t1-t0)/t0)>1e-12&(t1<dt&t1>0)
t0=t1;
t1=t0-(t0^2+a*t0+b)/(2*t0+a);
end
tm=t1;
end
else
p=[pa pb pc pd];
tm=roots(p);
tm=tm(imag(tm)==0 &tm<dt&tm>0);
if(imag(tm)~=0)
tt=1
end
if(length(tm)~=1)
tm=min(tm);
end
end
end
if(tm>dt|tm<0)
cuowu='111111?????'
end
txg(1)=xg(1);
txg(2)=xg(1)+(xg(2)-xg(1))*tm/dt;
if(tm==0)
xm=x1;
else
[xm]=newmark_single(tm,k1,m,eta,txg,x1);
end
%%%不計誤差則應該xm(1)=xs;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -