-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathduffing_HB.m
142 lines (105 loc) · 3.25 KB
/
duffing_HB.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
clc
clear
close all
% sampling time
Ts = 0.01;
% number of terms considered in the fourier series
n_terms = 31;
% known parameters
F = 2;
% ode generation tspan
tspan = 0:Ts:500;
% initial condition
y0 = zeros(2,1);
% use ODE45 to generate the data
[tval, yval] = ode45(@duffing, tspan, y0, odeset("RelTol",1e-6));
displacement = yval(:,1);
velocity = yval(:,2);
% length of the signal
L = length(tval);
% start of steady state
% truncate the first 2/5 of the signal
% determine the lag using cross correlation of the reconstructed signal and
% the original signal
% [crosscorr, lags] = xcorr(displacement, p3)
start = ceil((2/5)*length(displacement));
start = start+420;
% steady state displacement
displacement_ss = displacement(start:end);
% find the period of the signal if unknown
T = findPeriod(displacement_ss,Ts);
omega = 2*pi/T;
% length of the period in time steps
periodLength = T/Ts;
% One period from start index
OPFS = floor(start+periodLength);
% plots
figure(1)
plot(tval, displacement)
xlabel("$t$",'interpreter','latex')
ylabel('$x$','Interpreter','latex')
set(gca, 'fontsize',16)
figure(2)
plot(tval, velocity)
xlabel("$t$",'interpreter','latex')
ylabel('$\dot{x}$','Interpreter','latex')
set(gca, 'fontsize',16)
% Phase Portrait
figure(3)
plot(displacement(start:end),velocity(start:end))
xlabel("$x$",'interpreter','latex')
ylabel('$\dot{x}$','Interpreter','latex')
set(gca, 'fontsize',16)
% one cycle of displacement
figure(4)
plot(tval(start:OPFS),displacement(start:OPFS))
xlabel("$t$",'interpreter','latex')
ylabel('$x$','Interpreter','latex')
set(gca, 'fontsize',16)
% get the fourier series coefficients of the data
[A,B] = fourierCoeff(displacement(start:OPFS));
% one cycle
singlePeriod = displacement(start:OPFS);
% find the fourier series reconstruction of the cycle
fourierReconstructionDisp = zeros(1,length(singlePeriod));
t_fourier = (0:floor(periodLength))*Ts;
for i = 1:n_terms
fourierReconstructionDisp = fourierReconstructionDisp + A(i)*cos(2*pi*(i-1)*t_fourier/T) ...
+ B(i)*sin(2*pi*(i-1)*t_fourier/T);
end
figure(5)
plot(tval(start:OPFS), displacement(start:OPFS),'b-','linewidth',2)
hold on
plot(tval(start:OPFS), fourierReconstructionDisp,'r--','linewidth',3)
xlabel("$t$",'interpreter','latex')
ylabel('$x$','Interpreter','latex')
legend('original signal', 'fourier reconstruction')
set(gca, 'fontsize',16)
% HARMONIC BALANCE
p1 = zeros(length(tval),1);
for i = 1:n_terms
p1 = p1 + ((i-1)^2)*(A(i)*cos((i-1)*omega*tval)+B(i)*sin((i-1)*omega*tval));
end
p1 = -(omega^2)*p1;
p2 = zeros(length(tval),1);
for i = 1:n_terms
p2 = p2 + (i-1)*(-A(i)*sin((i-1)*omega*tval)+B(i)*cos((i-1)*omega*tval));
end
p2 = omega*p2;
p3 = zeros(length(tval),1);
for i = 1:n_terms
p3 = p3 + A(i)*cos((i-1)*omega*tval)+B(i)*sin((i-1)*omega*tval);
end
p4 = p3.^3;
p5 = F*cos(omega*tval);
G = [p1 p2 p3 p4];
D = G'*G;
r = (D\(G'))*p5;
figure(6)
plot(displacement(start:end),velocity(start:end),'b-','linewidth',3)
hold on
plot(p3,p2,'r--','linewidth',2)
xlabel("$x$",'interpreter','latex')
ylabel('$\dot{x}$','Interpreter','latex')
set(gca, 'fontsize',16)
legend('original signal','fourier reconstruction')