convolution integral with ode45

Hi guys, can someone help me solve this equation of motion using ode45 or any other way? I dont know how to approach the convolution integral. Thank you!

Answer 1

[Star Strider]

Try this:

odefcn = @(t,y,c,F,k,m,w) [y(2); -(k.*y(1) - F.*cos(t.*w) + integral(@(tau)c(t - tau).*y(2).*tau, 0, t))./m];
c = @(t) exp(-t.^2);                                                                    % Make Something Up

[F,k,m,w] = deal(rand,rand,rand,rand);                                                  % Replace With Correct Values

tspan = [0 10];
ics = [1; 1];

[t,y] = ode45(@(t,y)odefcn(t,y,c,F,k,m,w), tspan, ics);

figure
plot(t, y)
hl = legend('$u(t)$', '$\dot{u}(t)$');
set(hl,'Interpreter','latex')

Noite that ‘c’ actually has to be defined as a function in order for this to work. Supply the correct one.

I derived it with the help of the Symbolic Math Toolbox.

[Angie]

Thank you very much for your help! I want to ask one last thing.

integral(@(tau)c(t - tau).*y(2).*tau, 0, t))

at the integral above, why is it multiplied by tau?

[Star Strider]

As always, my pleasure.

Inside the convolution integral, ‘y(2)’ should actually be a function of ‘tau’, not ‘t’, although there is apparently no way to define it as such, since that would also require ‘u’ to be a function of both ‘t’ and ‘tau’, and it does not appear to be. So I experimented, and since ‘u’ (the Symbolic Math Toolbox converts all time functions to ‘y’), is actually not a function of ‘tau’ inside the integral and is instead a constant (in terms of ‘tau’), I simply multiplied it by ‘tau’ as a way around that. I know of no other way to deal with that problem.

You can always use this version:

odefcn = @(t,y,c,F,k,m,w) [y(2); -(k.*y(1) - F.*cos(t.*w) + integral(@(tau)c(t - tau).*y(2), 0, t))./m];

It simply treats ‘y’ as a function of ‘t’.

[Angie]

Got it. Thanks a lot again! Greatly appreciated!

[Star Strider]

As always, my pleasure!

[Torsten]

I doubt that is really a correct solution to the problem. The result of the integral evaluation will always be

integral_{tau=0}^{tau=t} u'(t)*c(t-tau) dtau, 

not

integral_{tau=0}^{tau=t} u'(tau)*c(t-tau) dtau

as required.

[Angie]

Star Strider also mentioned in the previous comments that thats why he multiplied integral(@(tau)c(t - tau).*y(2).*tau with tau. But yes you are right Torsten, it gives different solution when it is multiplied with tau and when it is not and maybe neither are correct.

[JINGYU KANG]

Star Strider's answer solves the different differential equation. That is,

where the desired equation to solve is as below

One can verify that by letting c(t) = 1 with k = 0, F = 0 for simplicity

[David Goodmanson]

Whether you include a factor of tau or not, I believe both of these solutions are incorrect, and every solution along these lines will be incorrect. The term is

Int{0,t} c(t-tau) u'(tau) dtau                    % u' instead of udot

i.e. a convolution that depends on the past history of u' from 0 to t. Replacing u'(tau) with u'(t) says that u' is a constant factor in the integrand (which it isn't). This gives

Int{0,t)} c(t-tau) u'(t) dtau =  u'(t) Int{0,t)} c(t-tau) dtau.

If you attemp to fix this up by inserting tau or any function f(tau) inside the integral, and then do the integration you get

u'(t) Int{0,t} c(t-tau) f(tau) dtau  = u'(t) g(t)

for some function g, which is a pointlike expression that does not involve past history of u' at all.

I have never unaccepted someone else's answer but would encourage Star Strider to take a hard look at that answer.

Answer 2

[Paul]

If c(t) has a Laplace transform then we can take advantage of the convolution property and convert to the s-domain, solve for U(s), and then convert back to the t-domain. For simple c(t) we can actually get a closed form expression for u(t), though a numerical or vpa solution is more likely going to be needed.

Define the differential equation

syms m k F omega t tau real
syms u(t) c(t)
du(t) = diff(u(t),t);
z(t) = int(c(t-tau)*diff(u(tau),tau),tau,0,t);
eqn = m*diff(u(t),t,2) + z(t) + k*u(t) == F*cos(omega*t)

eqn = 

Take the Laplace transform

Leqn = laplace(eqn)

Leqn = 

Sub in U(s) and C(s)

syms U(s) C(s)
Leqn = subs(Leqn,laplace([u(t),c(t)],t,s),[U(s),C(s)])

Leqn = 

Solve for U(s)

Leqn = isolate(Leqn,U(s))

Leqn = 

At this point we need the expresion for C(s). Assume a simple form of c(t) and sub C(s) into the expression for U(s)

c(t) = exp(-t);
Leqn = subs(Leqn,C,laplace(c(t),t,s))

Leqn = 

Simplify the expression for U(s)

[num,den] = numden(rhs(Leqn));
Leqn = lhs(Leqn) == num/den

Leqn = 

Choose some arbitrary parameters for illustration and sub in

mval = 1; kval = 2; Fval = 3; omegaval = 4;
Leqn = subs(Leqn,[m,k,F,omega,u(0),du(0)],[mval,kval,Fval,omegaval,5,6])

Leqn = 

If we only want a numerical evaluation of the solution, we can just use impulse

figure
[num,den] = numden(rhs(Leqn));
impulse(tf(double(coeffs(num,'all')),double(coeffs(den,'all'))),0:.01:20);grid

●

Get the solution for u(t)

u(t) = ilaplace(rhs(Leqn),s,t)

u(t) = 

u(t) = rewrite(rewrite(u(t),'expandsum'),'expandroot');

The solution includes some terms in i = sqrt(-1), but we know the solution has to be real. It proved too difficult to simplify u(t) into an expression with only real terms.

Derivative of u(t)

du(t) = diff(u(t),t);

Plot the real and imaginary parts of u(t) to verify the latter is zero, and plot the derivative of u(t)

figure;
fplot([real(u(t)),imag(u(t)),real(du(t))],[0,20]);grid
legend('u(t)','imag(u(t))','du(t)');

Note that u(0) and du(0) satisfy the initial conditions assumed above.

Now verify the solution satisfies the differential equation.

z(t) actually has a closed form expression

z(t) = int(c(t-tau)*du(tau),tau,0,t)

z(t) = 

Subtract the rhs of the differential equation from the lhs and sub in the parameters

check = subs(lhs(eqn) - rhs(eqn),[m,k,F,omega],[mval,kval,Fval,omegaval]);

Sub in the expressions for u(t) and c(t), note we could just set u(t) = real(u(t)) and z(t) = real(z(t)).

check = subs(check);

Check the result at some points in time, the answer should be zero

vpa(subs(check,t,0:10)).'

ans = 

[Paul]

Inspired by a conversation with someone else ... I suppose that one way to tackle this problem with an ode solver would be to use a fixed step with a method that only requires derivatives up to the previous time step, i.e., something like

w(t) = w(t - dt) + wdot(t-dt)*dt

Then the convoultion is always one step behind and it can be computed using stored values of xdot and the known function c(t) at every derivative calculation

t = 0:.01:20;
N = numel(t); dt = t(2);
m = 1; k = 2; F = 3; omega = 4;
w = zeros(2,N); % w(1) = x; w(2) = xdot
w(:,1) = [5;6];  
wdot = zeros(2,1);
for ii = 2:N
    wdot(1) = w(2,ii-1);
    z = conv(exp(-t(1:ii-1)),w(2,1:ii-1));
    z = z(ii-1)*dt;
    wdot(2) = (F*cos(omega*t(ii-1)) - k*w(1,ii-1) - z)/m;
    w(:,ii) = w(:,ii-1) + dt*wdot;
end
figure
plot(t,w),grid

Compare to the solution derived above,

y = impulse(tf([5,11,94,179,176],conv([1,0,16],[1,1,3,2])),t);
hold on
plot(t,y,'o')
legend('x','xdot','x from original answer')

[Torsten]

Looks the same:

t = 0:.01:20;
N = numel(t); dt = t(2);
m = 1; k = 2; F = 3; omega = 4; c = exp(-t);
w = zeros(2,N); % w(1) = x; w(2) = xdot
w(:,1) = [5;6];  
for j = 2:N
  w(1,j) = w(1,j-1) + dt*w(2,j-1);
  w(2,j) = w(2,j-1) + dt*(F*cos(omega*t(j-1)) - k*w(1,j-1) - trapz(t(1:j-1),fliplr(c(1:j-1)).*w(2,1:j-1)))/m;
end
plot(t,w)
legend('x','xdot')
grid on

Do you have an idea how to solve it with standard adaptive MATLAB ODE-integrators (dde23, ode45, ...) ?

[Paul]

yes, trapz would be another way to compute the convolution integral. I suppose one could even use integral.

No, I can't think of way to use ode45 (et. al.) (I've tried). I just don't see a way around the need to properly maintain the history of xdot(t) in those solvers. Would be very interested if anyone else can.

[Paul]

Here is a solution using the fixed-step RK4 method. At each minor step (evaluation of ki), the derivative function stores and uses only those values of udot that are evaluated at previous major steps of the solution and the current value of udot at each minor step. The convolution in the derivative calculation is computed with integral (admittedly might be overkill) and assumes linear interpolation in the udot variable.

Still don't see how this can adapted to the stock ode solvers (i.e. ode45, et. al.)

m = 1;
k = 2;
F = 3;
wF = 4;
a = 1;              % c = exp(-a*t)
h = 1e-2;
tspan = 0:h:20;
u = nan(numel(tspan),2);
icond = [5,6];
deriv = @(t,u,rk) fun(t,u,m,k,a,F,wF,rk);

u(1,:) = icond;

for ii = 2:numel(tspan)
    n = ii - 1;
    tn = tspan(n);
    un = u(n,:);
    k1 = deriv(tn    ,un         ,1);
    k2 = deriv(tn+h/2,un + h*k1/2,2);
    k3 = deriv(tn+h/2,un + h*k2/2,3);
    k4 = deriv(tn+h  ,un + h*k3  ,4);
    u(ii,:) = un + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
% Clear the persistent variables for subsequent run.
% Or take other action inside fun to initialize on first call.
clear fun 
figure
plot(tspan,u),grid

function udot = fun(t,u,m,k,a,F,wf,rk)
persistent thist udothist
if rk == 1
    thist    = [thist, t];
    udothist = [udothist, u(2)];
end

c = @(t) exp(-a*t);
if rk == 1
    udotfun = @(tt) interp1(thist,udothist,tt);
else
    udotfun = @(tt) interp1([thist,t],[udothist,u(2)],tt);
end

if t == 0
    z = 0;
else
    z = integral(@(tau) c(t-tau).*udotfun(tau),0,t);
end

udot(1) = u(2);
udot(2) = (F*cos(wf*t) - k*u(1) - z)/m;

end

[Paul]

Solution using ode45 with an OutputFcn as insightfully suggested by @Torsten, a 2-element tspan, and the 2-argument output form. In this formulation, the elements of udot are stored in the derivative function at time steps determined by Refine (default is Refine = 4 for ode45).

If tspan is specified with more than two elements (all else the same) then the points in tspan need to be spaced close enough together (whatever that means) to ensure an accurate solution because the values of udot will be stored in the derivative function only at the values in tspan.

Futher investigation needed if using the single output form of ode45, in which case Refine "does not apply," which I assume means it's effectively Refine = 1 as far as the OutputFcn is concerned.

The code could be improved by using a preallocation-fill scheme for the persistent variables instead of growing them dynamically.

m = 1;
k = 2;
F = 3;
w = 4;
a = 1;              % c = exp(-a*t)
c = @(t) exp(-a*t);
icond = [5;6];
tspan = [0,20];
[t,x] = ode45(@(t,x) fun(t,x,m,c,k,F,w), tspan, icond,odeset('OutputFcn',@OutputFcn,'OutputSel',2));

figure
plot(t,x(:,1))

[y,ty] = impulse(tf([5,11,94,179,176],conv([1,0,16],[1,1,3,2])),tspan(2));
hold on
plot(ty,y)
legend('ode45','impulse')

function xdot = fun(t,x,m,c,k,F,w,flag,tsol,udotsol)
persistent thist udothist
if nargin == 10
    if strcmp(flag,'init')
        thist    = tsol(1);
        udothist = udotsol;
    elseif isempty(flag)
        thist    = [thist,tsol];
        udothist = [udothist,udotsol];
    elseif strcmp(flag,'done')
        thist    = [];
        udothist = [];
    end
    return
end

if numel(thist) <= 1
    z = 0;
else
    udotfun = @(tt) interp1([thist,t],[udothist,x(2)],tt);
    z = integral(@(tau) c(t-tau).*udotfun(tau),0,t);
end

xdot    = zeros(2,1);
xdot(1) = x(2);
xdot(2) = (F*cos(w*t) - k*x(1) - z)/m;

end

function status = OutputFcn(tsol,udotsol,flag)
status = 0;
fun([],[],[],[],[],[],[],flag,tsol,udotsol)
end

[Torsten]

Your solution works well. Here is a suggestion where ode45 is called in a loop.

In some sense, it's again a solution with fixed time stepping - only within each "outer" time step, the "inner" solution is adaptive.

m = 1;
k = 2;
F = 3;
w = 4;
a = 1;              % c = exp(-a*t)
c = @(t) exp(-a*t);
icond = [5;6];
tspan = [0,20];

% Length of history vector 
nt = 100;

% Define vector of output times
T = linspace(tspan(1),tspan(2),nt);

% Initialize history vectors
th = tspan(1);
udoth = icond(2);

% Initialize solution vector at output times
U = icond;

% Call ode45 in loop
tstart = T(1);
for i = 1:nt-1
  tend = T(i+1);
  [t,u] = ode45(@(t,y)fun(t,y,m,k,F,w,c,th,udoth),[tstart tend],icond);
  % Update history vector
  th = [th,tend];
  udoth = [udoth,u(end,2)];
  % Update inital conditions
  icond = [u(end,1);u(end,2)];
  % Update solution vector
  U = [U,icond];
  % Update tspan
  tstart = tend;
end

% Compare with impulse response
[y,ty] = impulse(tf([5,11,94,179,176],conv([1,0,16],[1,1,3,2])),tspan(2));

% Plot results
hold on
plot(T,U(1,:))
plot(ty,y)
hold off
grid on
legend('ode45','impulse')

function dudt = fun(t,u,m,k,F,w,c,th,udoth)
  % Don't take into account incomplete last time interval
  %ch = c(th(end)-th);
  %z = trapz(th,udoth.*ch);
  % Take into account incomplete last time interval
  th = [th,t];
  udoth = [udoth,u(2)];
  ch = c(t-th);
  z = trapz(th,udoth.*ch);
  
  dudt    = zeros(2,1);
  dudt(1) = u(2);
  dudt(2) = (F*cos(w*t) - k*u(1) - z)/m;
end

[Paul]

I had been thinking about this approach but wanted to avoid having to define the vector of output times a-priori. I'd rather let the solver figure it out (in general, but there could be exceptions). Of course, the approach I illustrated has the same concern if the user wants to specify a full tspan vector. In that case, maybe better to use 2-element tspan and the structure output from ode45 with subesquent calls to deval, but I didn't check into that option.

In the numerical solutions posted in this thread, the convolution in the derivative function assumes linear interpolation of something. Do you think that more accuracy could be obtained by storing the history of u(t), its derivatives, and any other information and somehow using all of that to get better estimates of udot(t) between the stored points for the convolution operation?

[Torsten]

Do you think that more accuracy could be obtained by storing the history of u(t), its derivatives, and any other information and somehow using all of that to get better estimates of udot(t) between the stored points for the convolution operation?

"trapz" on an equidistant grid has integration order 2.

Maybe "interp1" with the "spline" option if you want to use "integral" for the convolution integral ?

And of course the absolute and relative tolerances for ode45 will influence the number of output points, thus the size of the history vector and thus the accuracy of the convolution integral.

Answer 3

[David Goodmanson]

Hi Paul, Torsten et. al.

The code below saves the history of udot in ode45 and then does the convolution. I used delt = 1e-2 as Matt did and the plot is pretty similar. Taking delt down to 1e-4 does not change things very much.

A more accurate method might be available and I will add it to this answer if I can get it to work.

One thing I don't like about the ode solvers is that for e.g a second order ode you get back y and y' but not y'' which you have calculated at great expense in the odefun but can't get to. With this method you could save up the calculated y'' values. The code below just used the saved stuff internally, so accss wasn't a problem. But to get the y'' vales out I can't think of a way other than using a global variable, which is not great, but it would be possible to have a wrapper function that contains the global and then passes the values out from there, so you wouln't have to have a global mucking up the command window.

m = 1;
k = 2;
F = 3;
wF = 4;
a = 1;              % c = exp(-a*t)
delt = 1e-2;        % step size for the convolution integration
tspan = [0 20];
icond = [5;6];
[t u] = ode45(@(t,u) fun(t,u,m,k,a,delt,F,wF,tspan,icond), tspan, icond);
figure(3); grid on
plot(t,u)

function udot = fun(t,u,m,k,a,delt,F,wF,tspan,icond)
persistent udotstore
persistent tstore
if t == tspan(1)
    udotstore = icond(2);
    tstore = tspan(1);
    G = 0;
else
    udotstore = [udotstore u(2)];
    tstore = [tstore t];
    % sort the times, toss out near-duplcate times
    [t1 ind] = sort(tstore);
    udot1 = udotstore(ind);
    fin = find(diff(t1)<1e-7);      % judgment call here
    t1(fin) = [];
    udot1(fin) = [];
    % set up equally spaced time array, do the convolution G
    t2 = t1(1):delt:t1(end);
    udot2 = interp1(t1,udot1,t2);
    c = exp(-a*t2);
    G = trapz(t2,flip(c).*(udot2));
end
udot = [0 1; (-k/m) 0]*u -(G/m)*[0;1] + (F/m)*cos(wF*t)*[0;1];
end

[Torsten]

@David Goodmanson

Storing the values for preceding times is a good idea.

The problem is that within function "fun", u is not the final solution at time t, but only serves to evaluate the right-hand side of the ODE system within the integration process.

So I think that either you have to call "ode45" in a loop for a fixed vector of output times and save the results when it returns the respective solutions or to use an OutputFcn function. This function is only called after a successful integration step and can be used to save the results so far as you did in "fun" above.

[Paul]

Hi David,

I'm afraid I have an issue with this approach at a fundamental level, even though the result seems to pass the eye test for the example problem.

As I undersand it, the code is storing all values of u(2) ( = udot) on every derivative evaluation by the ode45 solver and treating each of those values (as long as they are separated by 1e-7) as actual points in the solution of the underlying ode. However, they are not. As I know that you know, these solvers make a lot of calls to the derivative function for the full state vector and then the outputs of the derivative function are combined to form the full step of the state. At most, only one of those derivative evaluations is actually evaluated at a valid point in the solution of the ode, e.g., the K1 term in the classical fixed step RK4 algorithm.

To put it more concretely, consider the K2 and K3 terms in the RK4 algorithms. They are both evaluated at the same time, but they will have different values (neither of which are an actual point in the solution). Which one would should be used in the convolution?

And that's just for a fixed step solver. I'm sure this issue becomes more convoluted (sorry, couldn't resist) with these adaptive solvers.

To illustrate, I changed the model to a standard, second order model with c = 0.1, but used your code to store the udot history in the derivative function and then plot it at the end and show the final value of udot1. We see that:

@fun is called twice with t = tspan(end)

The values of udot1 at those times are not the same; so at least one of them is incorrect (as compared the actual solution).

The stored history values of udot in @fun do not match the actual solution output from ode45. The stored values have little bits and bobs that deviate from the actual solution.

In the problem as originally posed, we don't really know what happens because the effect that causes these deviations in this fixed-coefficient problem actually influences the solution in the original problem as posed. I'm not sure I've stated that clearly, hopefully the gist is clear.

The issues that I've raised and illustrated here are what I was refering to when I said " .. properly maintain the history of xdot(t)." (emphasis added). At least my opinion of "properly".

Having said all of this, and giving it a bit more thought, I think there may be a path forward for a fixed-step, multi-step solver, such as RK4, that assuages my concerns. I might come back later with an update to my answer. I believe the approach I'm thinking of could be used if one were to write one's own adaptive step solver, but I don't see how it could be used with the built-in ode suite of functions.

m = 1;
k = 2;
F = 3;
wF = 4;
a = 1;              % c = exp(-a*t)
delt = 1e-2;        % step size for the convolution integration
tspan = [0 20];
icond = [5;6];
[t u] = ode45(@(t,u) fun(t,u,m,k,a,delt,F,wF,tspan,icond), tspan, icond);

ans = "udot1 at tspan = 20 is -2.9486"

ans = "udot1 at tspan = 20 is -3.0678"

["u(2) at tspan = " + tspan(end) " is: " + u(end,2)]

ans = 1×2 string array
    "u(2) at tspan = 20"    " is: -3.0678"

hold on
plot(t,u(:,2))

●

function udot = fun(t,u,m,k,a,delt,F,wF,tspan,icond)
persistent udotstore
persistent tstore
if t == tspan(1)
    udotstore = icond(2);
    tstore = tspan(1);
    G = 0;
else
    udotstore = [udotstore u(2)];
    tstore = [tstore t];
    % sort the times, toss out near-duplcate times
    [t1 ind] = sort(tstore);
    udot1 = udotstore(ind);
    fin = find(diff(t1)<1e-7);      % judgment call here
    t1(fin) = [];
    udot1(fin) = [];
%    % set up equally spaced time array, do the convolution G
%    t2 = t1(1):delt:t1(end);
%    udot2 = interp1(t1,udot1,t2);
%    c = exp(-a*t2);
%     G = trapz(t2,flip(c).*(udot2));
end
%udot = [0 1; (-k/m) 0]*u -(G/m)*[0;1] + (F/m)*cos(wF*t)*[0;1];
c = 0.1;
udot = [0 1; (-k/m) -c/m]*u + (F/m)*cos(wF*t)*[0;1];

if t == tspan(end)
    figure
    plot(t1,udot1)
    ["udot1 at tspan = " + tspan(end) + " is " + udot1(end)]
end
end

[Torsten]

As said: not the ODE function, but an OutputFcn function should be used to store the values for u.

OutputFcn is called for increasing values of t and only after a basic integration steps has successfully finished (thus when u really contains the solution of the ODE at time t).

[Paul]

@Torsten

I didn't see this comment before I prepared and submitted mine. I see that we both have the same fundamental concern (I think).

I did not think of using an OutputFcn.

I think the derivative function can include additional arguments, which won't be used by the ode solver, but can be used in a call from the OutputFcn to pass data from the OutputFcn to the derivative function for the derivative function to store persistently. Don't see why that wouldn't work.

Or the state data could be stored persistently in the OutputFcn and it could be called from the derivative function with an optional second argument and specific value of flag or a fourth argument to return the stored state data. I think that should work too.

Answer 4

[Sam Chak]

In the past few days, I have been closely following the discussions on this intriguing topic of solving the integro-differential equation that involves a convolution integral term using the ode45 solver. It seems challenging to define the time-shifted dummy variable τ, which serves as a temporary placeholder for integration in the ODE function, as pointed out by @Star Strider.

I commend everyone for their contributions, including @JINGYU KANG, who revived this unsolved problem, and @David Goodmanson, who initially implemented storing the historical solutions from ode45 needed to perform the convolution. I would also like to especially recognize @Paul for developing solutions using the Laplace transform method, the fixed-step RK4 method, and the OutputFcn method (under @Torsten's insightful suggestion).

Here is my take: I initially wondered if the product of functions in the convolution integral term could be decomposed using integration by parts, a common mathematical technique that we all learned in high school calculus.

The idea eventually evolved into a mathematical technique inspired by integration by substitution.

Methodology:

The mass-spring-damper system is given by:

where is the convolution integral term.

If the damping function, is differentiable or Laplace transformable, we can introduce an auxiliary variable y to define that

Next, applying Leibniz integral rule for differentiating , we obtain

where the kernel function is

.

Evaluating the two terms:

and

.

Since

it can be rewritten simply as a coupled ODE:

.

Essentially, the integro-differential equation is converted into a coupled ODE system, allowing the ode45 solver to adaptively solve the system.

% Parameters
m       = 1;            % mass
k       = 2;            % stiffness
F       = 3;            % amplitude of the sinusoidal Force
w       = 4;            % natural angular frequency

% Call ode45 solver
tspan   = [0, 20];      % simulation time
icond   = [5            % u(0)  = 5
           6            % u'(0) = 6
           0];          % y(0)  = 0
[tu, u] = ode45(@(t, u) odefcn(t, u, m, k, F, w), tspan, icond);

% Compare with impulse response
[v, tv] = impulse(tf([5, 11, 94, 179, 176], conv([1, 0, 16], [1, 1, 3, 2])), tspan(2));

% Plot results
figure
hold on
    plot(tu, u(:,1), '.')
    plot(tv, v)
hold off
grid on
xlabel('Time')
ylabel('Amplitude')
legend('Adaptive ode45', 'Impulse response')

●

%% Transform the integro-differential equation into a coupled ODE system
function du = odefcn(t, u, m, k, F, w)
    du      = zeros(3, 1);
    
    du(1)   = u(2);
    du(2)   = 1/m*(F*cos(w*t) - k*u(1) - u(3));     % convolution integral is substituted by u(3) = y
    du(3)   = u(2) - u(3);                          % dy/dt = du/dt - y, after applying Leibniz rule
end

[Torsten]

Very nice solution.

What a luck that

d/dt (c(t-tau)) = -c(t-tau) 

in this specific case :-)

Can you think of a similar approach for a general convolution integral

int_{tau = 0}^{tau = t} c(t-tau) * u(tau) d(tau)

with an arbitrary function c ?