finding the value at a specific time from second-order ODE

I'm new to MATLAB

I've been working on solving ODE and I now have learnt two ways to plot the 2-order ODE

1. using ode45
2. using conv() fucntion

Here's the problem : I want to know the value at a specific time for example f(2.5) when t=2.5

it can be displayed in any way like output in the command window or on the figure

please teach a way to output the value

I'll attach the files I've been working on

if you think the codes are tltr you could just give suggestions though

[time,y] = ode45('ode45_function' , [0 10] , [0,0]) ;
figure ;
subplot(212);plot(time,y(:,1));xlabel('t');ylabel('y(t)');

m = 5;
k = 18;
c = 1.2;
f0 = 100;
wn = sqrt(k/m);
cc = 2*m*wn;
Dr = c/cc;
wd = wn*sqrt(1-Dr^2);
dt = 0.001; % assumed dt
time = 0:dt:10;
f = zeros(1, length(time));
f(0001:3000) = (f0/3) * (0:dt:3-dt); % Linear increase from 0 to 3 seconds
f(3001:5000) = f0; % Constant from 3 to 5 seconds
% f remains 0 after 5 seconds
g = (1/(m*wd))*exp(-Dr*wn*time).*sin(wd*time);
x = conv(f, g)*dt;
x = x(1:length(time)); % Ensure x is the same length as time

% x is the solution function for the 2-order ODE　

figure;
subplot(221); plot(time, f, 'r'); xlabel('t(s)'); ylabel('f(t)');
subplot(222); plot(time, g); xlabel('t(s)'); ylabel('g(t)');

subplot(223); plot(time, x); xlabel('t(s)'); ylabel('位移(m)');

% this is the function which I want to extract the specific value from 
% it for ex: f(3), f(2.5) ...

function ydot=ode45_function(t,y)
    m=5;
    k=18;
    c=1.2;
    
   
    if t<3
        f= 100*t/3;
    elseif t<5
        f=100;
    else
        f=0;
    end

   ydot(1)=y(2);
   ydot(2)=(1/m)*(f-c*y(2)-k*y(1));
   ydot=ydot';

end

Answer 1

[Steven Lord]

You could call ode45 with just one output argument and then call deval, or you could call ode45 with the desired time as one element of the timespan input argument, or (if you're using a sufficiently recent release of MATLAB) you could use the ode object instead of just ode45.

I also recommend you pass a function handle into ode45 rather than passing the name of the function.

One output with deval

sol = ode45(@ode45_function , [0 10] , [0,0]) ;
solutionAtTime = deval(sol, 2.5)

solutionAtTime = 2×1

5.2273 1.9079

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Timespan input contains the desired time

[t, y] = ode45(@ode45_function , [0 2.5 10] , [0,0]);
solutionAtTime2 = y(2, :).'

solutionAtTime2 = 2×1

5.2273 1.9079

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ODE object

odeobj = ode;
odeobj.ODEFcn = @ode45_function;
odeobj.InitialValue = [0 0];
results = solve(odeobj, 2.5); % Just solve at t = 2.5
solutionAtTime3 = results.Solution

solutionAtTime3 = 2×1

5.2267 1.9092

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Tien Yin]

APPRECIATE your hard work! Thanks alot!

[Tien Yin]

BY the way the two answers stand for {y} and {y'} is that correct ?

[Steven Lord]

The deval solution and the solution returned by calling solve on the ODE object are returned as columns of the matrix returned from those functions.

The solution returned by calling ode45 with two outputs and including the desired time in the time span returns the time vector as a column vector and the solution at each of those times in the corresponding row of y, so yes I did need to transpose the solution in that case to make it a column vector (so it looked like the solutions returned by the other two approaches.)

[Tien Yin]

thanks , while what　I meant by {y} and {y'} is position and velocity , I use quotation marks as a meaning of derivative

[Steven Lord]

Yes, if the first component of y represents position then the second component represents velocity.

[Tien Yin]

I get it , thanks for answering

Answer 2

[Torsten]

The easiest (maybe slightly inaccurate) way is to use interpolation:

f_conv = interp1(time,x,2.5)
f_ode45 = interp1(time,y(:,1),2.5)

[Tien Yin]

thanks now i learnt a new way to calculate the answer!

[Sam Chak]

@Tien Yin, Also try learning to "vote" other helpful answers or solutions.

[Tien Yin]

sure!