Info

この質問は閉じられています。 編集または回答するには再度開いてください。

Frequency sweep for DAE question

1 ビュー (過去 30 日間)
haohaoxuexi1
haohaoxuexi1 2021 年 11 月 30 日
終了済み: haohaoxuexi1 2022 年 1 月 19 日
%%
close all;
clear all;
clc;
format long;
k1=25000;
knl_1=300;
knl_3=2000000;
m=0.1;
c=1;
w0=20;
kAmp = 0.5;
M=[0,0,0; 0,1,0; 0,0,1];
opt = odeset('mass', M, 'RelTol',1.e-14);
a = w0-5;
b = w0+5;
np = 200;
f = linspace(a,b,np);
omega_rng = f;
y0 = [0; 0; 0]; % initial condition
for mm=1:length(omega_rng)
format long;
omega=omega_rng(mm);
w = omega;
period=2*pi/omega;
step=period/100;
endperiod=1000*period;
tspan=step: step: endperiod;
[t,usol] = ode23t(@(t,y) Test123(t, y, m, c, k1, knl_1, knl_3, kAmp, w), tspan, y0, opt);
y0=usol(end,:);
i=endperiod/period;
NUM1=1024;
step=2*pi/omega/NUM1;
tspan=[i*period: step: (i+1)*period];
Y1(mm)=norm(usol(:,2))/sqrt(1024);
[t,usol] = ode23t(@(t,y) Test123(t, y, m, c, k1, knl_1, knl_3, kAmp, w), tspan, y0, opt);
Y1(mm)=norm(usol(:,2))/sqrt(1024);
disp([' omega:', num2str(omega)]);
end
%%
%%
figure(2);
plot(f,Y1);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
%%
%%
function dy = Test123(t, y, m, c, k1, knl_1, knl_3, kAmp, w)
format long;
dy(1)=k1*(y(2)-y(1))-(knl_1*y(1)+knl_3*y(1)^3);
dy(2)=y(3);
dy(3)=(-c*y(3)-k1*(y(2)-y(1))+kAmp*cos(w*t))/m;
dy = dy';
end
%%
I want to do a frequency sweep for a DAE probelm. The code can work when I am not doing frequency sweep but calculating individual omega value. I know the reason why the code won't run is the new initial value didn't meet the requirement that the left side is 0 (it will be very close, but it may never reach).
dy(1)=k1*(y(2)-y(1))-(knl_1*y(1)+knl_3*y(1)^3);
Can anyone help me to solve the problem

回答 (0 件)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by