finding frequency and domain of equation using ode45

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

Answer 1

[Star Strider]

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:

[shahin hashemi]

tanx alot star strider

can u plz explain that in which line you use interpolation ? cause as i know ode 45 gives u nonuniform answer and u cant finde exact answer ?

[Star Strider]

As always, my pleasure!

I do not use interpolation. I let ode45 do that instead.

I define:

tv = linspace(0, 60, 5000);

as ‘tspan’ and that produces integrated results of the same length for each call to ode45. If ‘tspan’ has more than 2 elements, ode45 (and the other solvers) evaluate and return the integrated results at those time points. Choosing a vector of 5000 evenly-spaced points means that the indices returned in ‘zx2’ are close to the actual x-values (with a sampling interval of 0.012), so no interpolation is required there, either.

[shahin hashemi]

got it tanx :X

[Star Strider]

As always, my pleasure!

[shahin hashemi]

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]

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’.

.