osc1.m

文件是编辑好的m文件

function y=osc1(x0,tt);
ac=acos(x0(2));
 if x0(2)<1
     dd=sqrt(1-x0(2)^2);
     y=x0(1)*(1-exp(-x0(2)*x0(3).*tt).*sin(dd*x0(3).*tt+ac)/dd);
 elseif x0(2)==1
     y=x0(1)*(1-exp(-x0(3).*tt).*(1+x0(3).*tt));
 elseif x0(2)>1
     rr=sqrt(x0(2)^2-1);ee=1/x0(3)/(x0(2)-rr);ff=1/x0(3)/(x0(2)+rr);
     y=x0(1)*(1-ee*x0(3)/2/rr*exp(-tt/ee)+ff*x0(3)/2/rr*exp(-tt/ff));
 else y=x0(1)*(1-cos(x0(3)*tt));
 end;