Include convolution integral inside ode function

I have the following differential equation to describe free decay roll motion (as mass-spring-damper system) which I solved in Matlab:

With initial conditions:

A44 = 1.0e+07; % added mass 
Ixx = 3.0e+07; % inertia 
B44 =  2.00e+06; % damping
C44 = 1.0e+07; % spring 
dt = 0.2; % time increment [s]
t_end = 400; %time end
tspan = [0.2:dt:t_end]; % simulation time [s] - start,dt, end simulation time.
IC = [1.0;0]; % roll amplitude and velocity at 0s. 
odefun = @(t,z)ODE_function_case1( t,z, Ixx,A44,B44,C44)
[t_out,z_out] = odeRK4(odefun,tspan,IC); % Here I use RK4 method. ODE45 also works. 

My ODE function looks as such:

function [ zdot ] = ODE_function_case1( t,z, Ixx,A44,B44,C44)
x1= z(1); % motion
x2= z(2); % velocity
A_system= [1 0 ; 0 Ixx+A44];
B_system= [x2; -B44*x2-C44*x1];
zdot= A_system\B_system; % matrix inversion
end

My RK solver looks like:

function [t_out,Y ]= odeRK4(odefun,tspan,y0)
h = diff(tspan);
try
  f0 = feval(odefun,tspan(1),y0);
catch
  msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
  error(msg);  
end  
y0 = y0(:);   % Make a column vector.
if ~isequal(size(y0),size(f0))
  error('Inconsistent sizes of Y0 and f(t0,y0).');
end  
neq = length(y0);
F = zeros(neq,4);
N = length(tspan);
Y = zeros(neq,N);
Y(:,1) = y0;
for i = 2:N
  ti = tspan(i-1);
  hi = h(i-1);
  yi = Y(:,i-1);
  F(:,1) = feval(odefun,ti,yi);
  F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1));
  F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2));
  F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3));
  Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
  t_out(i,1) = ti;
end
Y = Y.';
end

This all works. Now I want to include a convolution integral into the differential equation:

The retardation function h(t) is known and is a vector length(tspan) * 1.

I am stuck how to solve this differential equation, the time history of the velocity should be used in this integral.

I am not sure how and where to do this. Does anyone know how to solve this challenging question?

How would I get the time history as calculated in odeRK4 to ODE_function_case1?

Help is greatly appreciated!

[Paul]

Can h(t) be expressed in a closed form?

[JvdS]

h(t) is another integral equation. it is a fourier transform of frequency dependent damping in this case.

With damping known values of size: length(omega) x 1

[Paul]

Can B(w) be expressed in closed form?

Answer 1

[Torsten]

Suppose you know the solution for velocity of your system of ODEs in

0*dt, 1*dt, 2*dt,..., n*dt.

To calculate the solution in (n+1)*dt, you have to evaluate the convolution integral

for t = n*dt with a value for velocity vndt,

for t = (n+1/2)*dt with two different values for velocity vn+1/2dt and

for t=(n+1)*dt with one value for velocity vn+1dt.

Approximation of the integral I for t=n*dt and velocities v(0),v(dt),...,v((n-1)*dt) and v(n*dt):

I = 0.5*dt*h(n*dt)*v(0) + dt*sum_{i=1}^{n-1} [h((n-i)*dt)*v(i*dt)] + 0.5*dt*h(0)*v(n*dt)

Approximation of the integral for t=(n+1/2)*dt and velocities v(0),v(dt),...,v(n*dt) and vn+1/2dt:

I = 0.5*dt*h((n+1/2)*dt)*v(0) + dt*sum_{i=1}^{n-1} [h((n-i+1/2)*dt)*v(i*dt)] + 0.75*dt*h(dt/2)*v(n*dt) + 0.25*dt*h(0)*vn+1/2dt

Approximation of the integral for t=(n+1)*dt and velocities v(0),v(dt),...,v(n*dt) and vn+1dt:

I = 0.5*dt*h((n+1)*dt)*v(0) + dt*sum_{i=1}^{n} [h((n+1-i)*dt)*v(i*dt)] + 0.5*dt*h(0)*vn+1dt

That's all.

Dividing your equation by 1e7 will make things easier, I guess.

To get the time history for the velocities, just call odeRK4 in a loop from 0 to dt, from dt to 2*dt,...

[JvdS]

I have to use the time history of velocities in ODE_function_case1 to calculate forces due to this convolution.

These velocities are calculated in odeRK4 as Y(2,2:N), then how exactly do I get this result of Y(2,2:N) in ODE_function_case1? I tried to declare it as global variable, but this screws up to solution after some time.

It looks like I have to write a custom output function with persistent variable, but cannot make this work somehow.

The working example without convoluition has the time history of velocities in z_out(:,2), but now their intermediate result during calculation has to be used in the convolution.

[Torsten]

In odeRK4, you will have to pass tspan(1:i-1) and Y(2,1:i-1) to odefun as third and fourth arguments and a flag whether you want to evaluate the convolution integral at tspan(i-1),tspan(i-1)+0.5*dt or tspan(i) as fifth argument. Then you can approximate the convolution integral in function ODE_function_case1 by the formulae from above.

By the way:

Why don't you simply write

zdot=[z(2);-0.05*z(2)-0.25*z(1)]

?

[JvdS]

Thanks for thinking along with me. Now I have stored the time history of velocity in ODE_function_case1 as:

persistent vel
vel = [vel;z(2)]; % time history of the velocity. 
nt = length(vel); % the number of the time step in the integration. 

I checked this with the solution and this works all well. I think I should only evaluate the convolution integral at timespan(nt), right? The odeRK4 already passes the correct velocity to ODE_function_case1.

[Torsten]

You have to evaluate the convolution integral at the times where you evaluate the function ODE_function_case1 according to the runge-kutta scheme, i.e. at times t = tspan(i) and t = tspan(i) + 0.5*dt.

[JvdS]

Yes, but how is this actually evaluated looking at the time history? Can you give an actual example code where you put the evaluation?

[Torsten]

I supplied the three approximating formulae for the convolution integral I depending on whether odeRK4 evaluates the right-hand side of your ODE at t=tspan(n), t=tspan(n)+0.5*dt and t=tspan(n+1).