From f6641a280e34f9c35c8c644a425ef1e8f36ed8cf Mon Sep 17 00:00:00 2001
From: "Nolte, David D" <nolte@purdue.edu>
Date: Wed, 17 Nov 2021 10:03:24 -0500
Subject: [PATCH] Add files via upload

Driven van der Pol
---
 drivenvanderpol.m | 28 +++++++++++++++++++++
 dvdptest.m        | 63 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 91 insertions(+)
 create mode 100644 drivenvanderpol.m
 create mode 100644 dvdptest.m

diff --git a/drivenvanderpol.m b/drivenvanderpol.m
new file mode 100644
index 0000000..31a62f6
--- /dev/null
+++ b/drivenvanderpol.m
@@ -0,0 +1,28 @@
+% function [y t] = drivenvanderpol(beta,mu,f0,fdriv,b,y0)
+% called by dvdptest.m
+
+function [y t] = drivenvanderpol(beta,mu,f0,fdriv,b,y0)
+
+w02 = (2*pi*f0)^2;
+w = 2*pi*fdriv;
+
+tspandum = [1 100];
+[tdum,ydum] = ode45(@f5,tspandum,y0);
+[sz dum] = size(tdum);
+
+y02 = [ydum(sz,1) ydum(sz,2) ydum(sz,3)];
+tspan = [1 200];        % 400
+[t,y] = ode45(@f5,tspan,y02);
+
+    function yd = f5(t,y)
+
+        yp(1) = y(2);
+        yp(2) = b*y(3) + 2*mu*y(2)*(1-beta*y(1)^2)-w02*y(1);
+        yp(3) = w*cos(w*t);
+        yd = [yp(1);yp(2);yp(3)];
+
+    end     % end f5
+
+
+end % end vandpol
+
diff --git a/dvdptest.m b/dvdptest.m
new file mode 100644
index 0000000..09bcbcb
--- /dev/null
+++ b/dvdptest.m
@@ -0,0 +1,63 @@
+% dvdptest.m
+
+clear
+
+y1 = 1;y2 = 1;y3 = 0;
+for floop = 1:50
+
+    fre = 0.4 + 0.025*(floop-1);
+
+    beta = 1;   %1
+    mu = 2.5;   %2.5
+    f20 = 1;    %1
+    b = 10;     %10
+
+    y0 = [y1 y2 y3];
+
+
+    [y t] = drivenvanderpol(beta,mu,f20,fre,b,y0);
+
+
+    [st dum] = size(y);
+
+    y1 = y(st,1);
+    y2 = y(st,2);
+    y3 = y(st,3);
+
+    figure(1)
+    plot(t,y(:,1))
+    pause(0.2)
+
+
+    tshort = t(400:st,1);
+    yshort = y(400:st,1);
+    [st dum] = size(tshort);
+    dt = (max(tshort)-min(tshort))/st;
+    df = 1/st/dt;
+    fmax = 1/dt;
+    freq = df:df:fmax;
+
+    yf = yshort(:,1);
+    F = fft(yf);
+    Pow = F.*conj(F);
+
+    mid = round(st/4);
+    figure(2)
+    loglog(freq(1:mid),Pow(1:mid))
+    title('Power Spectrum')
+    xlabel('Frequency (Hz)')
+
+
+    [Pmax I] = max(Pow(10:mid));
+    f0 = freq(I+1)
+    T0 = 1/f0;
+
+    fdep(floop) = f0;
+    fvar(floop) = fre;
+
+
+end
+
+figure(3)
+plot(fvar,fdep,'o','MarkerFaceColor','b')
+