Understanding continuous integrator blocks in Simulink

I'm struggling to understand how the continuous integrator block works. Let's say we are setting up a fixed step continuous case and we select a solver such as ode3. Ode3 uses the Bogacki-Shampine method. This method requires evaluating the deriviative y'=f(t,y) at specified values of t and y. And it seems Simulink evaluates at minor timesteps for those intermediate values of t.

This is quite straightforward if you have an actual function f(t,y). However, I'm don't understand how this works when the input isn't a function you can just evaluate at any point you wish (eg: an arbitrary input). And even if the block input were a function Simulink is integrating it numerically, not analytically, so you still have the same issue.

I tried building a simulink model of just an integrator block into C to see how it works, but i'm not great at C so I didn't really understand it. Can anyone help, or point me to resources you think might be useful?

Answer 1

[Sam Chak]

Hi @C

I don't what is the mathematical 'C' function in your differential equation. If the derivative function is unavailable, but you have the rate-of-change time series data, you will need to interpolate the data using the interp1() function. While I am unsure how this works in Simulink, the general idea is outlined below. This is my first attempt at implementing the Bogacki–Shampine method, so it may contain errors. Please check before using it.

%% The Script
% Parameters
c       = 1;

% Simulation time span
tStart  = 0;
h       = 0.2;   % step size
tEnd    = 10;
t       = (tStart:h:tEnd);

% Data
tdata   = 0:0.5:10;
dydata  = - exp(-c*tdata);

% Ordinary Differential Equation with Time-Dependent Terms
dy      = @(t) interp1(tdata, dydata, t);
odefcn  = @(t) dy(t);

% Initial value
y0      = 1;

% Call ode3 Solver
yode3   = ode3(@(t, y) odefcn(t), t, y0);

% True solution
ysol    = y0*exp(-c*t);

% Plot the solution
figure
hold on
plot(t, yode3)
plot(t, ysol, 'o')
hold off
grid on
xlabel('Time')
ylabel('Amplitude')
title('Bogacki–Shampine method')
legend('Numerical solution', 'True solution', 'location', 'best')

%% Bogacki–Shampine method (from your Wikipedia link)
function y = ode3(f, t, y0)
    y(:, 1) = y0;               % initial condition
    h       = t(2) - t(1);      % step size
    n       = length(t);        % number of steps
        for j  = 1 : n-1
            k1        = f(t(j), y(:, j));
            k2        = f(t(j) + 1/2*h, y(:, j) + 1/2*h*k1);
            k3        = f(t(j) + 3/4*h, y(:, j) + 3/4*h*k2);
            y(:, j+1) = y(:, j) + 2/9*h*k1 + 1/3*h*k2 + 4/9*h*k3;
        end
end

[C]

Wow, that was a very quick and thorough response! I didn't think someone would code up a whole example! Thanks so much for your time spent helping. I learned a few new things about matlab too, I haven't seen that @ function notation before!

However, I do have a followup question if you don't mind. To calculate K2 and K3 you need to evaluate the derivative function at intermediate steps of t and y. In your example y' = -exp(-c*t), the derivative is a function of only t (not y), so you appropriately just interpolate the input vector dydata at the intermediate timesteps.

If you instead chose a derivative function that was also a function of y, eg: y' = -exp(-c*t*y) you would need to also the derivative at intermediate values of not just t, but also of y (eg: y+1/2*h*k1). But y is what you're solving for, so unlike time you don't already have a vector y to interpolate. So how are you supposed to do that?

[Sam Chak]

Hi @C

It may be confusing for you that I used the solution itself to artificially generate the data for the demo previously. This time, only plain data is used. The code for solving the data-based derivative remains unchanged. The ode3 solver doesn't require the derivative equation at all, only the data.

%% The Script
% Parameters
c       = 1;

% Simulation time span
tStart  = 0;
h       = 0.2;   % step size
tEnd    = 10;
t       = (tStart:h:tEnd);

% Data
tdata   = [0.0	0.5	1.0	1.5	2.0	2.5	3.0	3.5	4.0	4.5	5.0	5.5	6.0	6.5	7.0	7.5	8.0	8.5	9.0	9.5	10.0];
dydata  = [-1	-0.60653066	-0.367879441	-0.22313016	-0.135335283	-0.082084999	-0.049787068	-0.030197383	-0.018315639	-0.011108997	-0.006737947	-0.004086771	-0.002478752	-0.001503439	-0.000911882	-0.000553084	-0.000335463	-0.000203468	-0.00012341	-7.48518E-05	-4.53999E-05];

% Ordinary Differential Equation with Time-Dependent Terms
dy      = @(t) interp1(tdata, dydata, t);
odefcn  = @(t) dy(t);

% Initial value
y0      = 1;

% Call ode3 Solver
yode3   = ode3(@(t, y) odefcn(t), t, y0);

% True solution
ysol    = y0*exp(-c*t);

% Plot the solution
figure
hold on
plot(t, yode3)
plot(t, ysol, 'o')
hold off
grid on
xlabel('Time')
ylabel('Amplitude')
title('Bogacki–Shampine method')
legend('Numerical solution', 'True solution', 'location', 'best')

●

%% Bogacki–Shampine method (from your Wikipedia link)
function y = ode3(f, t, y0)
    y(:, 1) = y0;               % initial condition
    h       = t(2) - t(1);      % step size
    n       = length(t);        % number of steps
        for j  = 1 : n-1
            k1        = f(t(j), y(:, j));
            k2        = f(t(j) + 1/2*h, y(:, j) + 1/2*h*k1);
            k3        = f(t(j) + 3/4*h, y(:, j) + 3/4*h*k2);
            y(:, j+1) = y(:, j) + 2/9*h*k1 + 1/3*h*k2 + 4/9*h*k3;
        end
end