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

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

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')

2 Comments

Omg.. Thank you sooo much. You've just saved my life! Thanks!!

Sign in to comment.

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!