Include convolution integral inside ode function

15 次查看（过去 30 天）
JvdS
JvdS 2021-12-16
评论： Torsten 2021-12-19
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!
  3 个评论
显示 1更早的评论
JvdS
JvdS 2021-12-17
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

请先登录，再进行评论。

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

回答（1 个）

Torsten
Torsten 2021-12-16
编辑：Torsten 2021-12-16
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,...
  6 个评论
显示 4更早的评论
JvdS
JvdS 2021-12-19
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
Torsten 2021-12-19
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).

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by