Info

This question is closed. Reopen it to edit or answer.

Frequency sweep for DAE question

5 views (last 30 days)
haohaoxuexi1
haohaoxuexi1 on 30 Nov 2021
Closed: haohaoxuexi1 on 19 Jan 2022 at 2:36
%%
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

Answers (0)

Community Treasure Hunt

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

Start Hunting!