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.

[Paul]

This example considers periodic c(t) with period T = 2. The fundamental period of c(t) is c0(t) = exp(-t) for 0 <= t <= T, and 0 otherwise.

Three solutions are considered:

a) Laplace transform followed by impulse

b) ode45 using the augmented state derived by differentiating the convolution integral, as suggested by @Sam Chak and @David Goodmanson, after expanding c(t) into a Fourier series.

c) ode45 storing the solution of udot in a persistent variable via the OutputFcn and numerically evaluating the convolution integral in the odefun as suggested by @Torsten

I'd be interested if anyone has references that discuss typical c(t) functions for this type of problem.

clearvars; assume(assumptions,'clear');
m = 1; k = 2; F = 3; omega = 4;
u0 = 5; udot0 = 6; 

Method (a)

syms 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);
syms U(s) C(s)
Leqn = laplace(eqn);
Leqn = subs(Leqn,laplace([u(t),c(t)],t,s),[U(s),C(s)]);
Leqn = subs(Leqn,[u(0),du(0)],[u0,udot0]);
Leqn = isolate(Leqn,U(s));
[num,den] = numden(rhs(Leqn));
Leqn = lhs(Leqn) == num/expand(den);
[num,den] = numden(rhs(Leqn));
numU = coeffs(num,C(s),'all');
denU = coeffs(den,C(s),'all');

T = 2;
c0(t) = exp(-t);

C(s) = (laplace(c0(t)*rectangularPulse(0,T,t),t,s)/(1-exp(-s*T)))

C(s) = 

sv = tf('s');
% for exp(-t)*rect(0,T,t)
C = (exp(-2)*exp(-sv*T) - 1)/((sv+1)*(exp(-sv*T)-1));
U = ...
    (tf(double(coeffs(numU(1),s,'all')),1)*C + tf(double(coeffs(numU(2),s,'all')),1)) / ...
    (tf(double(coeffs(denU(1),s,'all')),1)*C + tf(double(coeffs(denU(2),s,'all')),1));

tfinal = 50;
[uimp,timp] = impulse(pade(U,7),0:.05:tfinal);

Method (b)

syms n integer
f(n) = fourier(c0(t)*rectangularPulse(0,T,t),t,2*sym(pi)/T*n)/T;
f0 = limit(f(n),n,0);
f(n) = simplify(f(n));
w0 = 2*sym(pi)/T;
figure
nFourier = 30; % 61-term series is a reasonable approximation.
kk = 1:nFourier;kk = [kk,-kk];
fplot([c0(t)*rectangularPulse(0,T,t),f0 + sum(f(kk).*exp(1j*w0*kk*t))],[0,10])

f = matlabFunction(f);
f0 = double(f0);
w0 = double(w0);

% consider setting the MaxStep based on nFourier*w0
[tG,uG] = ode45(@(t,u) odefun(t,u,m,k,F,omega,u0,f0,f,w0),[0,tfinal],[u0;udot0;zeros(nFourier,1)]);

Method (c)

c0 = matlabFunction(c0);
c  = @(t) c0(mod(t,T));

[tout,uout] = ode45(@(t,u) fun(t,u,m,c,k,F,omega),[0,tfinal],[u0;udot0],odeset('OutputFcn',@OutputFcn,'OutputSel',2));

Compare results

figure
plot(timp,uimp,tG,uG(:,1),tout,uout(:,1)),grid
legend('Laplace','Fourier+diffconv','persistent')

function xdot = fun(t,x,m,c,k,F,w,flag,tsol,udotsol)
persistent thist udothist nhist
if nargin == 10
    if strcmp(flag,'init')
        if ~isempty(thist)
            error('thist is not empty at initialization');
        end
        if ~isempty(udothist)
            error('udothist is not empty at initialization');
        end
        thist = zeros(1,10000);
        udothist = zeros(1,10000);
        thist(1)    = tsol(1);
        udothist(1) = udotsol;
        nhist = 1;
    elseif isempty(flag)
        if nhist + numel(tsol) > numel(thist)
            disp('Growing storage!');
            thist = [thist,zeros(1,1000)];
            udothist = [udothist,zeros(1,1000)];
        end
        thist(nhist+1:nhist+numel(tsol)) = tsol;
        udothist(nhist+1:nhist+numel(tsol)) = udotsol;
        nhist = nhist + numel(tsol);
    elseif strcmp(flag,'done')
        thist    = [];
        udothist = [];
        nhist    = [];
    end
    return
end
if isempty(nhist)
    z = 0;
else
    udotfun = @(tt) interp1([thist(1:nhist),t],[udothist(1:nhist),x(2)],tt);
    %z = integral(@(tau) c(t-tau).*udotfun(tau),0,t);
    % could also use trapz
    dt = 0.01;
    tconv = 0:dt:t;
    z = conv(c(tconv),udotfun(tconv));
    z = dt*z(numel(tconv));
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

function udot = odefun(t,u,m,k,F,omega,u0,f0,f,w0)

nFourier = numel(u) - 2;
n = (1:nFourier).';
b = f(n);
a = -1j*w0*n;

udot        = zeros(numel(u),1);
udot(1)     = u(2);
udot(2)     = (F*cos(omega*t) - k*u(1) - f0*(u(1)-u0) - 2*sum(real(u(3:end))))/m;
udot(3:end) = -a.*u(3:end) + b.*u(2);

end

[Paul]

Suppose c(t) is defined over a finite interval t1 <= t <= t2 and is 0 outside that interval. In that case (and probably under mild assumptions), c(t) can be represented by a complex Fourier series. However, unlike when c(t) is periodic, the complex exponentials are windowed over the defining interval to represent c(t). So we need to explore how to implement the solution for c(t) being a windowed exponential (which extends to sums of those). Using the OutputFcn/Persistent approach (@Torsten) doesn't change for this type of c(t), so we focus on how to use one of the built in solvers (@David Goodmanson, @Sam Chak)

clearvars; assume(assumptions,'clear');

Define c(t) as a finite-duration exponential

syms t tau real
syms t1 t2 positive
syms a b
c(t) = b*exp(-a*t)*rectangularPulse(t1,t2,t);

u(t) and its derivative

syms u(t) du(t)  % du(t) == diff(u(t),t)

G(t) is the convolution integral

syms G(t)
eqn1 = G(t) == int(c(t-tau)*du(tau),tau,0,t)

eqn1 = 

Differentiate the convolution

eqn = diff(eqn1,t);

And simplify it

eqn = lhs(eqn) == simplify(rhs(eqn),10)

eqn = 

ccc = children(rhs(eqn));
eqn = lhs(eqn) ==  subs(expand(ccc{1}),-a*expand(rhs(eqn1)),-a*lhs(eqn1)) + rewrite(ccc{2},'heaviside') + rewrite(ccc{3},'heaviside')

eqn = 

The derivative of the convolution contains two delayed terms. So we can use ode45 in three steps (0 -> t1, t1->t2, t2 -> t3) and adjust the odefun over those intervals, and we'd have to store udot persistently as well. Alternatively, we can use dde23, at least I think it applies to this type of problem. I've never used dde23 before, so thought I'd give it try.

clearvars
% Parameters of the system
m = 1; k = 2; F = 3; omega = 4;
u0 = 5; udot0 = 6;
a = 1/2; b = 5; t1 = 1; t2 = 3;  % c(t) = b*exp(-a*t) t1 <= t <= t2, 0 otherwise
lags = [t1,t2];
tspan = [0,10];
c0 = 5; % Additional, constant damping in the system; keeps it stable.

% solve with dde23
sol = dde23(@(t,y,Z) ddefun(t,y,Z,m,k,F,omega,t1,t2,a,b,c0),lags,@(t)[u0,udot0,0],tspan,ddeset('Jumps',[t1,t2],'RelTol',1e-6,'AbsTol',1e-9));

function ydot = ddefun(t,y,Z,m,k,F,omega,t1,t2,a,b,c0)

ylag1 = Z(:,1);
ylag2 = Z(:,2);

ydot    = zeros(3,1);
ydot(1) = y(2);
ydot(2) = (F*cos(omega*t) - k*y(1) - c0*y(2) - y(3))/m;
ydot(3) = -a*y(3) + b*ylag1(2)*exp(-a*t1)*(t>=t1) - b*ylag2(2)*exp(-a*t2)*(t>=t2);

end

% function for c(t) and udot(t)
c    = @(t) b*exp(-a*t).*(t >= t1 & t <= t2);
udot = @(t) deval(sol,t,2);

% extract the solution from the mesh output
dtsol = 0.01;
tsol = 0:dtsol:10;
[ysol,ypsol] = deval(sol,tsol);

% plot the solution and its derivative
figure
plot(tsol,ysol),legend('y(t)','ydot(t)','G(t)');

figure
plot(tsol,ypsol),legend('ydot(t)','ydotdot(t)','Gdot(t)');

% compute the convolution G(t) from the solution
G = nan(1,numel(tsol));
for ii = 1:numel(tsol)
    G(ii) = integral(@(tau) c(tsol(ii)-tau).*udot(tau),0,tsol(ii));
end

% verify the solution
figure
plot(tsol,m*ypsol(2,:) + k*ysol(1,:) + c0*ysol(2,:) + G - F*cos(omega*tsol)),grid

The verification seems reasonable, but I really have no idea what to expect from this approach.

Would be interested in any feedback on this use of dde23.

Answer 3

[David Goodmanson]

******************** MODIFIED ******(**************

because the method suggested here does not work. I took out some useless text but left the basic idea as an example of what not to do.

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 it works.

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.

[David Goodmanson]

I see that persistent is not going to work. Thanks to Torsten and especially Paul for taking the time to point out the flaws in that idea.

Going back somewhere on this thread I mentioned Matt in conndection with delta t = .01 but I meant to say Torsten.

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 ?

[Sam Chak]

Hi @Torsten

This high-school integration technique clearly does not work for arbitrary functions. If this technique cannot be applied, then the convolution must be evaluated numerically while integrating the ODE function, as demonstrated by @Paul using the OutputFcn method. In theory, this should be analogous to solving a Volterra integral equation. Are there any examples available in the MathWorks Help Center or official documentation?

Answer 5

[David Goodmanson]

This is a very interesting problem! Convolution of the loss term for the harmonic oscillator. I had been working on some extensions of the same approach that Sam posted so I will post these. If

c = b e^(-at)
G = Int{0,t} b e^(-a(t-tau)) u'(tau) dtau

then

(d/dt + a) G(t) = b u'(t)

and the system consists of that result combined with

mu'' + G + ku = F cos(wt)

Let

y = [u u' G]'              % two different uses of '

The odefun is

ydot = [ 0   1   0
       -k/m  0 -1/m
         0   b  -a ] y + (F/m) cos(wt) [0 1 0]'

with initial condition G(0) = 0, all of which was covered by Sam. An advantage of this method is that as well as u and u' you get the convolution integral as a function of time.

In everything below, the initial condition for every G is zero.

Some generalizations of c are possible. These are sums of functions of the form

% b t^n e^(-at)

which include e.g.

% b1 e^(-a1 t) + b2 e^(-a2 t) + b3 e^(-a3 t) .... sum of exponentials   b)
% ( ... b3 t^3 + b2 t^2 + b1 t + b0 ) e^(-a t)    polymomials           a)

and mixtures of those. If additional easily calculated c functions don't come around, then one task would be to fit an arbitrary c function to sums of these functions There are more possibilites for success there than it might seem. For example in addition to [ t e^(-at) ], the function [ e^(-a1 t) - e^(-a2 t) ] also equals 0 at the origin and has one peak.

a) For polynomials you can define the convolutions

 G_n = Int{0,t} b (t-tau)^n e^(-a(t-tau)) u'(tau) dtau

in which case

(d/dt + a) G_n) = n G_(n-1)
(d/dt + a) G_0) = b u'                  % as before

For a fourth degree polynomial with

y = [u u' G4 G3 G2 G1 G0]'

then after a straightforward process the odefun looks like

ydot=  [ 0     1    0   0  0  0  0
       -k/m    0  -1/m  0  0  0  0
         0     b4  -a   1  0  0  0
         0     b3   0  -a  1  0  0
         0    2b2   0   0 -a  1  0
         0    6b1   0   0  0 -a  1
         0   24b0   0   0  0  0 -a] y  + (F/m) cos(wt) [0 1 0 0 0 0 0]'

A better loooking result is obtained by rescaling the G's with

Gnew_n = Gold_n / (4-n)!

in which case

ydot=  [ 0    1    0   0  0  0  0
       -k/m   0  -1/m  0  0  0  0
         0    b4  -a   1  0  0  0
         0    b3   0  -a  2  0  0
         0    b2   0   0 -a  3  0
         0    b1   0   0  0 -a  4
         0    b0   0   0  0  0 -a] y  + (F/m) cos(wt) [0 1 0 0 0 0 0]'

b) An easier process is the sums of terms of the same type. For e.g. three exponentials the odefun is

ydot = [ 0    1    0     0     0
       -k/m   0  -1/m  -1/m  -1/m
         0    b1  -a1    0     0  ]
         0    b2   0    -a2    0  ]
         0    b3   0     0    -a3 ] y  + (F/m) cos(wt) [0 1 0 0 0]'

[Paul]

Hi David,

In case there's any interest, Matlab can help with the differentiation of the convolution as shown by you and @Sam Chak

syms t tau real
syms n integer
syms u(t) c(t) G(t,n)
du(t) = diff(u(t),t);
eqn1 = G(t,n) == int(c(t-tau)*du(tau),tau,0,t)

eqn1 = 

eqn2 = diff(eqn1,t)

eqn2 = 

syms a b
c(t) = t^n*b*exp(-a*t);
eqn3 = subs(eqn2)

eqn3 = 

eqn4 = lhs(eqn3) == expand(rhs(eqn3))

eqn4 = 

The first term on the RHS is -a*G(t,n), the second terms is n*G(t,n-1). We could do a little more work to get the software to express the RHS explicitly in terms of G(t,n) and G(t,n-1). Note that the third term is 0 unless n == 0, in which case (with Matlab convention that 0^0 = 1)

eqn5 = subs(eqn4,n,0)

eqn5 = 

eqn4 and eqn5 are what you derived in the middle of your post.

This statement "then one task would be to fit an arbitrary c function to sums of these functions" got me thinking. I have no idea what types of c(t) functions are of interest for these types of problems. In my example, I chose c(t) = exp(-t) just for simplicity. Now we have your type (a) and (b) functions and, as you say, mixtures of those. I'm not familiar with curve fitting with exponentials and so don't know what the possibilities are there. But some other functions might be of interest.

a) If c(t) is sinusoidal (or sums thereof), then it can expressed as a sum of complex sinusoids, which gets us back to your type (b). I've never tried use an ode solver for a system with complex coefficients; is there any reason to think it wouldn't work (doing appropriate things to maintain a real solution)?

b) If c(t) is periodic, then we can expand it in a complex exponential Fourier series. Hopefully, we can truncate to a finite number of terms in the series to maintain a good approximation to c(t), and we're back to type (b).

I have no idea if any these types of c functions are applicable to these types of problems, nor if the proposed solutions can work, but I found the thought experiment interesting.

[Paul]

I ran an example with c(t) = 1 + cos(t) and it worked just fine after expanding cos(t) into its complex exponential form.

[David Goodmanson]

Adding the 1 to make c positive seems like a good idea since if c has some negative regions one might run the risk that the convolution goes negative and sets the oscillator into exponential growth.

I wasn't sure if you mean you used the complex form and combined terms within the odefun, or whether you ran the original ode with a force of (F/m)*exp(i*w*t) and took the real part afterwards, which is legit since all the operations in the set of equations are real.

[Paul]

I used the complex form of c(t) = 1 + cos(t). Even though c(t) has two exponential terms that correspond to two G states, those G states are conjugates so the ode only needs to be augmented with one additional (complex-valued) state variable. The constant in c(t) can be handled explicitly.

Code below is for c(t) = 1 + cos(2*t) and includes commented lines for an earlier version that computed both G states.

m = 1; k = 2; F = 3; omega = 4;
u0 = 5; udot0 = 6; 
syms t tau
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);
syms U(s) C(s)
Leqn = laplace(eqn);
Leqn = subs(Leqn,laplace([u(t),c(t)],t,s),[U(s),C(s)]);
Leqn = subs(Leqn,[u(0),du(0)],[u0,udot0]);
Leqn = isolate(Leqn,U(s));
[num,den] = numden(rhs(Leqn));
Leqn = lhs(Leqn) == num/expand(den);
c(t) = (1 + cos(2*t));
C(s) = laplace(c(t),t,s);
Leqn = subs(Leqn);

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

[num,den] = numden(rhs(Leqn));
[uimp,timp] = impulse(tf(double(coeffs(num,'all')),double(coeffs(den,'all'))),20);

[t,u] = ode45(@(t,u) odefun(t,u,m,k,F,omega,u0),[0,20],[u0;udot0;0]);

figure
plot(timp,uimp,t,u(:,1),'.'),grid
legend('Laplace','ode45')

function udot = odefun(t,u,m,k,F,omega,u0)

% c(t) = 1 + cos(2*t)
% exponential form of cos(2*t)
% cos(2*t) = exp(-t*2i)/2 + exp(t*2i)/2
a1 = -2i; b1 = 1/2;
%a2 = 2i;  b2 = 1/2;

% convolution with constant term
% int(1*udot(tau),tau,0,t) = u(t) - u(0)

udot = zeros(3,1);
udot(1) = u(2);
%udot(2) = (F*cos(omega*t) - k*u(1) - (u(1)-u0) - real(u(3) + u(4)))/m;
udot(2) = (F*cos(omega*t) - k*u(1) - (u(1)-u0) - 2*real(u(3)))/m;
udot(3) = -a1*u(3) + b1*u(2);
%udot(4) = -a2*u(4) + b2*u(2);

end

[David Goodmanson]

Hi Paul

Impressive agreement. I tried running your code (I have the Signal Processing toolbox) and got the following error message:

---------

Incorrect number or types of inputs or outputs for function 'tf'.

Error in ode_example1 (line 21)

[uimp,timp] = impulse(tf(double(coeffs(num,'all')),double(coeffs(den,'all'))),20);

--------------

and I guess it's saying it wants a digital filter object, I dunno.

[Paul]

Hi David,

I used the tf and impulse functions from the Control System Toolbox.

AFAIK, the Signal Processing Toolbox does not have a function to compute the impulse response of a continuous-time system. I guess the next best thing would be to extract the continuous-time transfer function coefficents as above and then use impinvar or bilinear to convert to discrete-time and then use impz for the impulse response.