# finding frequency and domain of equation using ode45

2 views (last 30 days)
shahin hashemi on 9 Jun 2020
Commented: Star Strider on 17 Jul 2020
dear all
i use following code to find answer of the following equation :
u ̈+u+u^3=0
function dydt= vdp1(t,u)
dydt=[u(2);-u(1)-((u(1))^3)];
clc
clear all
for a=0.1:0.1:0.3
[t,y]=ode45(@vdp1,[0 60],[0 a]);
hold on
plot(t,y(:,1))
end
is there any way to find frequency and domain of this equation ? i know ode 45 gives nonuniform answer but can i use interpolation to finde the maximum of domain and The intersection with the x axis
In summary i want to find exact amount of red and green dot

Star Strider on 9 Jun 2020
Try this:
vdp1 = @(t,u) [u(2);-u(1)-((u(1))^3)];
tv = linspace(0, 60, 5000);
a=0.1:0.1:0.3;
zci2 = @(v) find(diff(sign(v)));
for k = 1:numel(a)
[t,y]=ode45(@vdp1,tv,[0 a(k)]);
ym(:,k) = y(:,1);
lmx(:,k) = islocalmax(y(:,1));
lmn(:,k) = islocalmin(y(:,1));
zx2(:,k) = find(diff(sign(y(:,1))));
end
figure
for k = 1:3
subplot(3,1,k)
plot(t,ym(:,k))
hold on
plot(t(lmx(:,k)),ym(lmx(:,k),k), '^r') % Plot Maxima
plot(t(lmn(:,k)),ym(lmn(:,k),k), 'vg') % Plot Minima
plot(t(zx2(:,k)),ym(zx2(:,k),k), 'dk') % Plot Zero-Crossings
hold off
grid
end
producing:

Star Strider on 9 Jun 2020
As always, my pleasure!
shahin hashemi on 17 Jul 2020
dear star srider sry to bother you again
again tanx for ur all help
if it is possible i have another question :
can u plz say y when i increase a from 0.3 to for example 0.4 i have the following error :
Unable to perform assignment because
the size of the left side is 20-by-1
and the size of the right side is
21-by-1.
Error in Untitled2 (line 11)
zx2(:,k) =
find(diff(sign(y(:,1))));
i run the following code :
clc
clear all
tv = linspace(0, 60, 5000);
a=0.1:0.1:0.4;
zci2 = @(v) find(diff(sign(v)));
for k = 1:numel(a)
[t,y]=ode45(@vdp1,tv,[0 a(k)]);
ym(:,k) = y(:,1);
lmx(:,k) = islocalmax(y(:,1));
lmn(:,k) = islocalmin(y(:,1));
zx2(:,k) = find(diff(sign(y(:,1))));
A=t(zx2(7,k))-t(zx2(5,k));
B=t(zx2(9,k))-t(zx2(7,k));
C=t(zx2(11,k))-t(zx2(9,k));
O(k)=(A+B+C)/3;
L=max(ym);
end
tanx alot again for all your attention
Star Strider on 17 Jul 2020
When ‘a’ is 0.4, there are more zero-crossings, so ‘zx2’ no longer has the same row size.
Creating ‘zx2’ as a cell array works, however much of that loop then has to be rewritten.
Try this:
tv = linspace(0, 60, 5000);
a=0.1:0.1:0.4;
zci2 = @(v) find(diff(sign(v)));
for k = 1:numel(a)
[t,y]=ode45(@vdp1,tv,[0 a(k)]);
ym(:,k) = y(:,1);
lmx(:,k) = islocalmax(y(:,1));
lmn(:,k) = islocalmin(y(:,1));
zx2{:,k} = find(diff(sign(y(:,1))));
zxk = zx2{:,k};
A=t(zxk(7))-t(zxk(5));
B=t(zxk(9))-t(zxk(7));
C=t(zxk(11))-t(zxk(9));
O(k)=(A+B+C)/3;
L=max(ym);
end
That appears to be robust to various values of ‘a’.
.

R2018a

### Community Treasure Hunt

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

Start Hunting!