convolution integral with ode45

Question

Angie 2019-6-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/466124-convolution-integral-with-ode45

编辑： Paul about 11 hours 前

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2019-6-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/466124-convolution-integral-with-ode45#answer_378375

在 MATLAB Online 中打开

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

JINGYU KANG 2023-9-12

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 2025-12-11，1:16

编辑：David Goodmanson 2025-12-11，7:35

在 MATLAB Online 中打开

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 2025-12-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/466124-convolution-integral-with-ode45#answer_1572498

编辑：Paul 2025-12-11，2:18

在 MATLAB Online 中打开

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 =

10 个评论
显示 8更早的评论隐藏 8更早的评论

Paul 2025-12-13，14:34

编辑：Paul 2025-12-13，15:09

在 MATLAB Online 中打开

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 2025-12-14，15:56

在 MATLAB Online 中打开

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 2025-12-15，13:31

编辑：Torsten 2025-12-15，13:48

在 MATLAB Online 中打开

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 2025-12-20，18:25

在 MATLAB Online 中打开

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 2025-12-31，19:20

编辑：Paul about 11 hours 前

在 MATLAB Online 中打开

Suppose c(t) is defined over a finite interval t1 <= t <= t2 and is 0 outside that interval. The OutputFcn/Persistent approach (@Torsten) and Laplace transform approach do not change for this type of c(t), so we focus on how to use one of the built in solvers w/o the need for persistent storage (@David Goodmanson, @Sam Chak).

For this class of c(t) (and probably under mild assumptions), c(t) can be represented by a complex exponential 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).

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 step 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 2025-12-13，10:37

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/466124-convolution-integral-with-ode45#answer_1572658

编辑：David Goodmanson 2025-12-16，9:19

在 MATLAB Online 中打开

******************** 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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Paul 2025-12-13，13:47

编辑：Paul 2025-12-14，19:49

在 MATLAB Online 中打开

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

Paul 2025-12-13，22:17

编辑：Paul 2025-12-13，22:46

@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 2025-12-16，9:11

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 2025-12-15，18:56

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/466124-convolution-integral-with-ode45#answer_1572702

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Torsten 2025-12-15，20:05

编辑：Torsten 2025-12-15，20:17

在 MATLAB Online 中打开

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 2025-12-17，13:38

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 2025-12-16，9:04

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/466124-convolution-integral-with-ode45#answer_1572709

编辑：David Goodmanson about 18 hours 前

在 MATLAB Online 中打开

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. To begin with I think it's better to redefine constants in the convolution loss term.

MODIFIED below to redfine constants so as to eliminate factors of (1/m) that were floating around

For the standard harmonic oscillator,

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

and defining R/m = g (gamma) k/m = w0^2 then

u'' + g u' + w0^2 u = (F/m) cos(wt) % g has units of freq^(-1)

For the convolution version, start with a convolution integral involving a function R(t-tau) and define R/m = c. Then

u'' + G(t) + w0^2 u = (F/m) cos(wt)
with 
G(t) = Int{0,t} c(t-tau) u'(tau) dtau        % c has units of freq^(-2).

For the case at hand

c = b e^(-at)

then

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

and the system consists of that result combined with

u'' + G + w0^2 u = (F/m) cos(wt)

Let

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

The odefun is

ydot = [ 0   1   0
       -w0^2 0  -1
         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
       -w0^2   0   -1   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
       -w0^2  0   -1   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
       -w0^2  0   -1    -1   - 1
         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]'

6 个评论
显示 4更早的评论隐藏 4更早的评论

Paul 2025-12-16，23:40

编辑：Paul 2025-12-17，2:48

在 MATLAB Online 中打开

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 2025-12-18，12:48

在 MATLAB Online 中打开

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 2025-12-19，9:28

编辑：David Goodmanson 2025-12-19，9:43

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 2025-12-19，13:13

编辑：Paul 2025-12-19，20:29

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.

请先登录，再进行评论。

convolution integral with ode45

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（4 个）

10 个评论
显示 8更早的评论隐藏 8更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

2 个评论
显示无隐藏无

6 个评论
显示 4更早的评论隐藏 4更早的评论

另请参阅

类别

标签

Community Treasure Hunt

convolution integral with ode45

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（4 个）

10 个评论 显示 8更早的评论隐藏 8更早的评论

5 个评论 显示 3更早的评论隐藏 3更早的评论

2 个评论 显示 无隐藏 无

6 个评论 显示 4更早的评论隐藏 4更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

2 个评论
显示无隐藏无

6 个评论
显示 4更早的评论隐藏 4更早的评论