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;