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 and if it is imposible i really appreciate if someone can help me finde the frequance and domain of this equation

Answer 1

[Star Strider]

Try this:

vdp1 = @(t,u) [u(2); -u(1)-((u(1))^3)];
tspan = linspace(0, 60, 240);

a=0.1:0.1:0.3;
for k = 1:numel(a)
    [t,y{k}]=ode45(@vdp1,tspan,[0 a(k)]);
end
ym = cell2mat(y);

figure
plot(t,ym(:,1:2:size(ym,2)))
grid
xlabel('Time')
lgdstrc = sprintfc('a = %.1f',a);
title('Time Domain Plot')
legend(lgdstrc)

Ts = mean(diff(tspan));                             % Sampling Interval
Fs = 1/Ts;                                          % Sampling Frequency
Fn = Fs/2;                                          % Nyquist Frequency
L = numel(tspan);                                   % Signal Length
FTvdp1 = fft(ym(:,1:2:size(ym,2)))/L;               % Fourier Trasnsform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                 % Frequency Vector
Iv = 1:numel(Fv);                                   % Index Vector

figure
plot(Fv, abs(FTvdp1(Iv,:)))
grid
xlabel('Frequency')
title('Frequency Domain Plot')
legend(lgdstrc)
xlim([0 1.2])

.

[shahin hashemi]

tanx alot Star Strider

i think its work

but if it is possible i have some question :

1- y do u pick 240 point in tspan

2-what is Nyquist Nyquist Frequency

3-and last:

i really appreciate if u could help me with FV formulation

[shahin hashemi]

i mean :

Fv = linspace(0, 1, fix(L/2)+1)*Fn;

[Star Strider]

1. That was arbitrary. It creates 4 points for every time unit (from 0 to 60). The ‘tspan’ vector can be almost anything you want it to be.

2. The Nyquist frequency is one-half of the sampling frequency. It is the highest frequency that can be uniquely resolved in a sampled signal.

3. The ‘Fv’ vector is the vector of frequencies for the one-sided fft plot (and other possible analyses that could be done). It creates a vector that goes from 0 to the Nyquist frequency, and is one element longer than half the length of the fft, to be certain that it spans 0 Hz to the Nyquist frequency. The fix call is necessary in case ‘L’ is odd so it will work for any ‘L’ in the event you want to experiment with the code. (The ‘Iv’ vector is related to that, and is used to be certain that the vector lengths match.)

.

[shahin hashemi]

tanx alot for ur answer

is there an way to have more smooth curve ?

[Star Strider]

As always, my pleasure!

To get better frequency resolution, the easiest way is to increase the length of the fft.

Experiment with this:

Ts = mean(diff(tspan));                             % Sampling Interval
Fs = 1/Ts;                                          % Sampling Frequency
Fn = Fs/2;                                          % Nyquist Frequency
L = numel(tspan);                                   % Signal Length
N = nextpow2(L);                                    % Increasing ‘fft’ Length Increases Frequency Resolution
FTvdp1 = fft(ym(:,1:2:size(ym,2)),2^N)/L;           % Fourier Trasnsform
Fv = linspace(0, 1, fix(2^N/2)+1)*Fn;               % Frequency Vector
Iv = 1:numel(Fv);                                   % Index Vector

The rest of the code is unchanged.

.

[shahin hashemi]

tanx a again u are my hero :D

i search the net and is it called the Frequency response for hardening plot ?

[Star Strider]

As always, my pleasure!

I have only heard of it as a frequency response plot. I do not know what ‘hardening’ is in this context.

[Pedro Calorio]

I have one question.

In your function used to solve the ode, how do you consider the second order nature of the equation?

As written, d^2u/dt^2 + u + u^3 = 0. You defined as:

vdp1 = @(t,u) [u(2); -u(1)-((u(1))^3)];

tspan = linspace(0, 60, 240);

How does it works? I'm kind of lost here

And the funcion requires a u vector as input. But where do you insert these values? It is when you define the a variable to use as initial value?