在下面的myfun2及主函式中，想要實作將ug處的方程形式變為一組地震動加速度激勵形式進行求解，
不知應該如何將地震動施加入方程中，還希望各位老師可以指點迷津，小弟在此謝過:D
撰寫的myfun2函式部分如下：
function dy = myfun2(t,y,x1,x2,g,H,alpha)
B = tan(alpha)*H;
R = sqrt(B^2+H^2);
p = sqrt(3*g/4/R);
wp = x1*p;
ap = x2*g*tan(alpha);
phi = asin(alpha*g/ap);
ug = ap*cos(wp*t).*(t<=2);
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -p^2*( sin(alpha*sign(y(1))-y(1)) + ug/g*cos(alpha*sign(y(1))-y(1)) );
end
主函式部分如下：
clear
clc
close all
g = 9.8;
H = 7.11;
alpha = 15/180*pi;
B = tan(alpha)*H;
R = sqrt(B^2+H^2);
p = sqrt(3*g/4/R);
T = 0:0.0001:10;
x1 = pi/p;
x2 = 0.53/tan(alpha);
aa = 0;
cc = 0;
theta = [];
T_new = T;
while length(T_new)
[t,y] = ode23s(@(t,y)myfun2(t,y,x1,x2,g,H,alpha),T_new,[aa,cc]);
index = find(sign(y(1:end-1,1)).*sign(y(2:end,1))==-1);
if isempty(index)
theta = [theta;y];
break;
else
theta = [theta;y(1:index(1)-1,:)];
T_new = T_new(index(1):end);
aa = y(index(1)+1,1);
cc = y(index(1),2)*(1-3/2*(sin(alpha))^2);
end
end
wp = x1*p;
ap = x2*g*tan(alpha);
phi = asin(alpha*g/ap);
ug = ap*cos(wp*T).*(T<=2);
subplot(3,1,1)
plot(T(1,1:length(ug))',ug/g)
axis([0,10,-0.5,0.5])
subplot(3,1,2)
plot(T,theta(:,1)/alpha)
axis([0,10,-1.5,1.5])
subplot(3,1,3)
plot(T,theta(:,2)/p)
axis([0,10,-0.5,0.5])
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/115288.html
標籤：其他開發語言
上一篇：7-4計算運算式求解答