How can I plot a graph with a time dependent ode45?

4 views (last 30 days)
clc
clear all
A = 0.05; %Excitation Amplitude
l_r = 2; %Wave length of the road
v = 45; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
k_l = 26400; %Linear stiffness
m = 483; %Mass
d = -0.1; %Stretching condition
l = 0.5; %Length of the spring
k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
f = @(t,x) [ x(2); ...
-(4*k_s*(x(1)-A*sin(Om*t))*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2) - l))/ ...
(m*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2))) ];
T = 50;
x0 = [0,0];
[t,x] = ode45(f,[0,T],x0);
Response_amp = max(x(:,1)) - min(x(:,1));
plot(t,x(:,1))
xlabel('Time (s)')
ylabel('Amplitude (m)')
title('When d=-0.1', 'fontsize', 20)
set(gca,'FontSize',15)
Hi, all. This is my ode45 code. If I chagne A (Excitation Amplitude), Response_amp will have a different value. I wish to plot this relationship as A vs Response_amp.

Accepted Answer

Ameer Hamza
Ameer Hamza on 16 Mar 2020
Edited: Ameer Hamza on 16 Mar 2020
You can visualize such relation with 3D plots
l_r = 2; %Wave length of the road
v = 45; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
k_l = 26400; %Linear stiffness
m = 483; %Mass
d = -0.1; %Stretching condition
l = 0.5; %Length of the spring
k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
fig = figure();
ax = axes();
hold(ax);
view([-53 33]);
grid on
A_array = 0.05:0.05:0.3; %Excitation Amplitude
T = 15;
x0 = [0,0];
for i=1:numel(A_array)
A = A_array(i);
f = @(t,x) [ x(2); ...
-(4*k_s*(x(1)-A*sin(Om*t))*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2) - l))/ ...
(m*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2))) ];
[t, x] = ode45(f,[0,T],x0);
Response_amp = max(x(:,1)) - min(x(:,1));
plot3(t, A*ones(size(t)), x(:,1), 'LineWidth', 1);
end
xlabel('Time (s)')
ylabel('A (m)')
zlabel('Amplitude (m)')
title('When d=-0.1')

More Answers (0)

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!