# Get all graphs from ODE23

3 views (last 30 days)
Ali Awada on 29 Dec 2020
Commented: Ali Awada on 30 Dec 2020
Hello all,
I am trying to write a code using the ode23 to solve for a free vibration damped SDOF. I wrote the following code without using the ode23 function:
function f=free_viscous(m,c,k,x0,v0,n,tspan);
t=[0:tspan:tspan*n];
wn=sqrt(k/m);
cc=2*m*wn;
zeta=c/cc;
phi=atan((v0+zeta*wn*x0)/(x0*wn*sqrt(1-zeta^2)));
X=sqrt(x0^2*wn^2+v0^2+2*x0*v0*zeta*wn)/(wn*sqrt(1-zeta^2))
syms d;
x=X*cos(wn*d*sqrt(1-zeta^2)-phi)*exp(-zeta*wn*d);
v=diff(x,d)
a=diff(v,d)
for i=1:101
x(i)=X*cos(wn*t(i)*sqrt(1-zeta^2)-phi)*exp(-zeta*wn*t(i));
v(i)=X*wn*exp(-t(i)*wn*zeta)*sin(phi - t(i)*wn*(1 - zeta^2)^(1/2))*(1 - zeta^2)^(1/2) - X*wn*zeta*exp(-t(i)*wn*zeta)*cos(phi - t(i)*wn*(1 - zeta^2)^(1/2));
a(i)=X*wn^2*exp(-t(i)*wn*zeta)*cos(phi - t(i)*wn*(1 - zeta^2)^(1/2))*(zeta^2 - 1) + X*wn^2*zeta^2*exp(-t(i)*wn*zeta)*cos(phi - t(i)*wn*(1 - zeta^2)^(1/2)) - 2*X*wn^2*zeta*exp(-t(i)*wn*zeta)*sin(phi - t(i)*wn*(1 - zeta^2)^(1/2))*(1 - zeta^2)^(1/2);
end
plot(t,x,t,v,t,a);
Running code for the previous is(i was able to get all of my graphs):
clc
clear all
free_viscous(450,1000,26519.2,0.539657,1,100,0.025)
I was able to write the ode23 code but i didnt know how to use MATLAB to plot the acceleration. I managed only the displacement and velocity:
function f=dfunc2(t,x);
m=450;
c=1000;
k=26519.2;
f=zeros(2,1);
f(1)=x(2);
f(2)=(-c/m)*x(2)-(k/m)*x(1);
the main code is:
clc
clear all
tspan=[0:0.025:2.5];
x0=[0.539657;1];
[t,x]=ode23('dfunc2',tspan,x0);
plot(t,x(:,1),t,x(:,2));
Thanks!

Star Strider on 30 Dec 2020
Try this approach:
function f=dfunc2(t,x);
m=450;
c=1000;
k=26519.2;
f=zeros(2,1);
f(1)=x(2);
f(2)=(-c/m)*x(2)-(k/m)*x(1);
end
tspan=[0:0.025:2.5];
x0=[0.539657;1];
[t,x]=ode23(@dfunc2,tspan,x0);
acc = gradient(x(:,2))./gradient(t); % Calculate Acceleration By Taking The Numerical Derivative Of Velocity
figure
plot(t,x(:,1),t,x(:,2), t,acc);
grid
.

Jan on 30 Dec 2020
Vectorize the function to be integrated: Then it accepts a matrix with different X values also:
tspan = 0:0.025:2.5;
x0 = [0.539657; 1];
[t, x] = ode23(@dfunc2, tspan, x0);
% Call dfunc2 to get the accelerations:
a = dfunc2(t, x.').'; % x must be transposed twice
figure;
plot(t, x(:,1), t, x(:,2), t, a(:, 2));
end
function f = dfunc2(t, x)
m = 450;
c = 1000;
k = 26519.2;
f = [x(2, :); ...
(-c/m) * x(2, :) - (k/m) * x(1, :)];
end
a(:, 1) is the same as x(:, 2) as expected.
Ali Awada on 30 Dec 2020
Thanks for the detailed explanation. Definitely a better approach!