finding frequency and domain of equation using ode45

9 次查看（过去 30 天）
shahin hashemi
shahin hashemi 2020-4-16
评论： Pedro Calorio 2020-8-13
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

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2020-4-16
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])
.
  8 个评论
显示 6更早的评论
Star Strider
Star Strider 2020-4-17
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
Pedro Calorio 2020-8-13
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?

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品


版本

R2017a

Community Treasure Hunt

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

Start Hunting!

Translated by