Info

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

Frequency sweep for DAE question

1 view (last 30 days)
haohaoxuexi1
haohaoxuexi1 on 30 Nov 2021
Closed: haohaoxuexi1 on 19 Jan 2022
%%
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)

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!