Most effective way to solve nonhomogeneous linear ODE problem

What is the most effective way to solve following "small" linear 1st order ODEs problem:

x'(t) = Ax(t) + Bu(t)
x(t0) = x0

where A, B are (2x2) real matrices with constant coefficients , and u(t) is represented by discrete (measured) real signals.

u(t_i) = [u_1(t_i);u_2(t_i)].

The requirements:

• fast as possible ... repeated numerical solution for different signals u(t) with the same A and B.
• signals u(t) may contains some additive noise
• coefficients at matrices A and B may differs by many magnitudes of order ... stiff problem

[Torsten]

I think you don't have a big choice between methods.

I only see ode15s or ode45 in combination with interp1 to interpolate your signal to the actual time during integration.

[Michal]

@Torsten ode15s or ode45 are good in general case, but in my specific case are too slow.

What about the direct use of general integral form of the system solution with matrix exponentials, which are possible to precompute very fast for (2x2) matrix A.

Answer 1

[Paul]

Hello Michal,

Following up on this comment, I'd be concerned about the solution of the "reference" method if that's what's being used to compare to the actual implementation.

Define the parameters of the problem as in the linked comment, but extend the time vector for a few more points and include an intial condition.

rng(100);
A  = rand(2);    % homogeneous matrix A
t  = (0:20);     % time vector
Nt = length(t);  % number of samples
u  = rand(1,Nt); % signal vector
B  = rand(2,1);  % non-homogeneous vector B 
x0 = rand(2,1);  % initial condition

Note that A has both eigenvalues in the RHP, so the system is unstable. Not sure how/if that impacts the accuracy of any of the solutions.

Solve the problem by the reference method in the linked comment:

xa = zeros(2,Nt);
xa(:,1) = x0;
for ii = 2:Nt
    % standard numerical integration
    xa(:,ii) = expm(A*(t(ii) - t(1)))*x0 + integral(@(s) expm(A*(t(ii)-s))*B*interp1(t,u,s), t(1), t(ii), 'ArrayValued',true);
end

My concern here is that integral() doesn't know about the breakpoints in the interp1() of u and so could miss changes of direction in u at the breakpoints.

Here is the reference method from the linked comment, but tell integral() about the breakpoints

xb = zeros(2,Nt);
xb(:,1) = x0;
for ii = 2:Nt
    xb(:,ii) = expm(A*t(ii) - t(1))*x0 + integral(@(s) expm(A*(t(ii)-s))*B*interp1(t,u,s), t(1), t(ii), 'ArrayValued',true,'Waypoints',t);
end

Here, build up the solution incrementally so that integral() only ever needs to integrate between two breakpoints in u.

xc = zeros(2,Nt);
xc(:,1) = x0;
for ii = 2:Nt
    xc(:,ii) = expm(A*(t(ii)-t(ii-1)))*xc(:,ii-1) + integral(@(s) expm(A*(t(ii)-s))*B*interp1(t,u,s), t(ii-1), t(ii), 'ArrayValued',true);
end

If the input data is equally spaced in time, the above can be made more efficient because the first call to expm() and be hoisted out of the loop because it's only a function of dt.

If the input data is equally spaced in time, use lsim() with the foh method, which is built for inputs that vary linearly with time between the input samples. BTW, on my local system this method is the fastest by far. It's just iterataing a difference equation and doesn't need a bunch of calls to expm() (it might need one or two to get the difference equation coefficient matrices, I'm not sure)

xd = lsim(ss(A,B,eye(2),0),u,t,x0,'foh').';

Finally, something like ode45 may be consdidered as the "ground truth".

xe = zeros(2,Nt);
xe(:,1) = x0;
for ii = 2:Nt
    temp = ode45(@(s,x) (A*x + B*interp1(t,u,s)),[t(ii-1) t(ii)],xe(:,ii-1),odeset('MaxStep',1e-4));
    xe(:,ii) = temp.y(:,end);
end

Compare the errors between the methods and the ode45 solution

format short e
[sum(vecnorm(xe-xa,2,1)) sum(vecnorm(xe-xb,2,1)) sum(vecnorm(xe-xc,2,1)) sum(vecnorm(xe-xd,2,1))]

ans = 1×4

6.0091e+00 2.8455e-05 1.6259e-05 1.4903e-05

Also, because A is just a 2x2, it would be very straightforward to compute the matrix exponential symbolically once ahead of time, and therefore avoid all the calls to expm().

[Michal]

@Paul Thanks for interesting comment.

Regarding reference solution of my problem, I am afraid, that the ode45 or any other solver based on any kind of approximation (finite-difference) can not be considered as "ground truth" at all. ode45 is only approximate solver (of course, very well working in many cases, but still the approximation).

From my point of view, the real reference solution should be based on general analytic integral form solution of system x'(t) = A*x(t) + B*u(t) , with respect of linearly interpolated segments between sample points of discrete signal u(t).

Now I am using as reference solution the following code:

xr = zeros(2,Nt);

I_i = zeros(size(A));
xr(:,1) = x0;

for i = 2:Nt
    I_i     = I_i + integral(@(s) expm(A*(t(i)-s))*(a(i)*s+b(i)), t(i-1), t(i), 'ArrayValued',true);;
    xr(:,i) = expm(A*(t(i)-t(1)))*x0 + I_i*B;     
end

where each linear 2-point segment is represented as u(s) = ai*s+bi, where ai is the slope and bi is the intercept of the i-th linear segment.

a = diff(u)/dt;                             % slopes
b = u(1:end-1) - a.*(0:(length(u)-2))*dt;    % intercepts

Of course expm is based on some approximations, too. But it is far more accurate, than any finite-diference approximation, like ode45.

Moreover, integral can be splited on two parts (slope and intercept part) and then is possible to find analytic solutions in closed form, at least for A (2 x 2) case.

[Torsten]

Regarding reference solution of my problem, I am afraid, that the ode45 or any other solver based on any kind of approximation (finite-difference) can not be considered as "ground truth" at all. ode45 is only approximate solver (of course, very well working in many cases, but still the approximation).

This is ridiculous (sorry) in view of a noisy u vector and its piecewise approximation in the integral.

[Michal]

So, you mean that some FD approximation with 2-point linear interpolation of noisy signal u (like ode45 solver) is always better than analytic solution with the same signal u interpolation??? Who is ridiculous (sorry) here?

I am wondering how do you know, that in this special case is ode45 solver better than analytical solution, regardless of whether the signal U is noisy or not.

[Torsten]

Due to your problem description, both methods - the mix of analytical and numerical (analytical integration + interp1) as well as the pure numerical (ode45 + interp1) - are infeasible.

I just found it funny that you argued an analytical method for a noisy data vector on the right-hand side yields better results because a finite difference method in integration is only an approximation. Using an analytical method together with noisy data cannot save the result.

[Michal]

So, there is no "ground truth" solution (at least from numerical implementation point of view)!

ode45 + interp1 vs. analytical integration + interp1 both methods suffer from the same problem ... interpolation between samples. But there is no reason to prefer ode45 as a more accurate "ground truth" method. Moreover ode45 needs sub-sampling (by interp1). On the other hand analytical integration with linear 2-point segment interpolation can be performed without any sub-sampling via possibility to find analytical closed form integrals.

Finally, I am asking one's again: Is there any feasible method in a case of discrete signal u with a noise (real measurements), which is very common in practice?

[Bruno Luong]

I answeredyou already : apply suitable filtering or regularization.

[Paul]

@Michal

Let me preface this comment by saying that the options offered in my Answer were all predicated on the assumption that you are willing to approximate the noisy input with linear interpolation between the input samples. I made this assumption because it was exactly the assumption you made in what was then called the reference code in this comment. It is in this context, i.e., with this assumption, that I stated that the ode45 solution may be considered as "ground truth". The quotation marks mean not to take that statement literally. The reason I included the ode45 method is because it was clearer to me how that solution works. Another reason to consider the ode45 approach is because it can be used as comparison to the various expm() methods in this thread, as will be done below. Whether or not integral() or ode45() is more accurate for this particular problem, under the stated assumption, is not something I can comment on. Probably depends on the ode45() and integral() options.

I hate to say it, but now I have a concern about the current reference solution in this comment. I've added that approach to the code to illustrate the concern, corrected it, and added a symbolic approach for your consideration.

Here is the code from the Answer:

rng(100);
A  = rand(2);    % homogeneous matrix A
t  = (0:20);     % time vector
Nt = length(t);  % number of samples
u  = rand(1,Nt); % signal vector
B  = rand(2,1);  % non-homogeneous vector B 
x0 = rand(2,1);  % initial condition

xa = zeros(2,Nt);
xa(:,1) = x0;
for ii = 2:Nt
    % standard numerical integration
    xa(:,ii) = expm(A*(t(ii) - t(1)))*x0 + integral(@(s) expm(A*(t(ii)-s))*B*interp1(t,u,s), t(1), t(ii), 'ArrayValued',true);
end

xb = zeros(2,Nt);
xb(:,1) = x0;
for ii = 2:Nt
    xb(:,ii) = expm(A*t(ii) - t(1))*x0 + integral(@(s) expm(A*(t(ii)-s))*B*interp1(t,u,s), t(1), t(ii), 'ArrayValued',true,'Waypoints',t);
end

xc = zeros(2,Nt);
xc(:,1) = x0;
for ii = 2:Nt
    xc(:,ii) = expm(A*(t(ii)-t(ii-1)))*xc(:,ii-1) + integral(@(s) expm(A*(t(ii)-s))*B*interp1(t,u,s), t(ii-1), t(ii), 'ArrayValued',true);
end

xd = lsim(ss(A,B,eye(2),0),u,t,x0,'foh').';

xe = zeros(2,Nt);
xe(:,1) = x0;
for ii = 2:Nt
    temp = ode45(@(s,x) (A*x + B*interp1(t,u,s)),[t(ii-1) t(ii)],xe(:,ii-1),odeset('MaxStep',1e-4));
    xe(:,ii) = temp.y(:,end);
end

Here is the new code.

First, the current reference code. I had to make some mods to make it work, hopefully it runs the way you intended

xf = zeros(2,Nt);
I_i = zeros(size(A));
dt = t(2) - t(1);
a = diff(u)/dt;                              % slopes
b = u(1:end-1) - a.*(0:(length(u)-2))*dt;    % intercepts
a = [0 a];  % zero pad a and b, so that a(ii) and b(ii) define the input from t(ii-1) to t(ii)
b = [0 b];
xf(:,1) = x0;
for ii = 2:Nt
    I_i     = I_i + integral(@(s) expm(A*(t(ii)-s))*(a(ii)*s+b(ii)), t(ii-1), t(ii), 'ArrayValued',true);
    xf(:,ii) = expm(A*(t(ii)-t(1)))*x0 + I_i*B;     
end

Compare the first few and last points of xf to xe (which was basically the same as xb, xc, and xd).

format short e
xe(:,[1:5 end])

ans = 2×6

1.0e+00 * 2.5243e-01 1.1847e+00 4.1890e+00 1.3160e+01 3.9055e+01 1.0619e+09 7.9566e-01 2.0998e+00 5.8547e+00 1.6821e+01 4.8831e+01 1.3158e+09

xf(:,[1:5 end])

ans = 2×6

1.0e+00 * 2.5243e-01 1.1847e+00 4.1614e+00 1.2882e+01 3.7755e+01 1.0221e+09 7.9566e-01 2.0998e+00 5.8325e+00 1.6582e+01 4.7605e+01 1.2666e+09

We see that the solutions at t(1) and t(2) are the same, but diverge after that. In short, I'm pretty sure you can't break up the integration of I like that because t is the upper limit of integration and it is also in the integrand.

The correction is to integrate over each step separately (not cumulatively), and use the final condition from the previous step as the initial condition of the current step, as is done in the solution for xc

xg = zeros(2,Nt);
I_i = zeros(size(A));
dt = t(2) - t(1);
a = diff(u)/dt;                              % slopes
b = u(1:end-1) - a.*(0:(length(u)-2))*dt;    % intercepts
a = [0 a];
b = [0 b];
xg(:,1) = x0;
for ii = 2:Nt
    I_i     = integral(@(s) expm(A*(t(ii)-s))*(a(ii)*s+b(ii)), t(ii-1), t(ii), 'ArrayValued',true);
    xg(:,ii) = expm(A*(t(ii)-t(ii-1)))*xg(:,ii-1) + I_i*B;     
end

Interestingly (or unsurprisingly?), xc and xg are identical.

isequal(xc,xg)

ans = logical
   0

Finally, here's an approach that uses the exact solution for the integral (evaluated in double rather than carrying symbolic math all the way through). Maybe this should be considered "ground truth," under the stated assumption, insofar as there is no need for any approximate numerical solution via integral() or ode45().

syms tsym ssym asym bsym t1sym t2sym
E(tsym) = expm(sym(A)*tsym);
I = int(E(t2sym - ssym)*sym(B)*(asym*ssym + bsym),ssym,t1sym,t2sym);
Ifunc = matlabFunction(vpa(I));
Efunc = matlabFunction(vpa(E));

xh = zeros(2,Nt);
xh(:,1) = x0;
for ii = 2:Nt
    I_i     = Ifunc(a(ii),b(ii),t(ii-1),t(ii));
    xh(:,ii) = Efunc(t(ii) - t(ii-1))*xh(:,ii-1) + I_i;     
end

The same exp() terms show up multiple times in E and I, so some efficiency could be gained by taking the extra step of only computing them once and then using those results multiple times.

Compare the errors between the methods and the symbolic solution

format short e
[sum(vecnorm(xh-xa,2,1)) sum(vecnorm(xh-xb,2,1)) sum(vecnorm(xh-xc,2,1)) sum(vecnorm(xh-xd,2,1)) sum(vecnorm(xh-xe,2,1)) sum(vecnorm(xh-xf,2,1)) sum(vecnorm(xh-xg,2,1))]

ans = 1×7

1.0e+00 * 6.0091e+00 3.3342e-05 3.3634e-05 3.5038e-05 4.9316e-05 9.6336e+07 3.3634e-05

If all you have are samples of the input to a continuous system, then some assumption will be needed on the behavior of the input between the samples or on the properties of the input itself, like its bandwidth, or some other action will be needed, as has been suggested by @Bruno Luong.

[Michal]

@Paul OK ... my cumulative integration is obviously completely wrong!!! Thanks for help...

Answer 2

[Star Strider]

I would use the Control System Toolbox lsim function with the system object created with the ss function and your (A,B,C,D) matrices.

[Michal]

Thanks for the hint. I will try this method.

But, in a case of small matrices A and B ... (2x2)!!!, is simply possible to compute matrix exponentials by its analytical form (very fast) and general solution is then in the form

x(t) = expm(A*(t-t0))*x0 + Integral(expm(A*(t-s))*B*u(t)*ds),

where

s = t0:dt:t  

I am wondering if is this approach more effective and robust than other more general methods, for this specific case.

[Bruno Luong]

expm matrix A, for constant is exp of of the eigen space component.

So you should call eig on A and works on the eigen space.

[Michal]

@Bruno Luong Yes, you are right!!! But how to perform this task really effectively? Could you please elaborate your recomendation in a bit more details?

You propose to use:

[V,D] = eig(A)
expmAt = V*diag(exp(diag(D)*t))/V % is equal to expm(A*t)

Where V and D for constant matrix A are precomputed? Am I right?

[Bruno Luong]

Exactly. You can even simplify more

A=rand(2)

A = 2×2

0.1128 0.5315 0.0007 0.4671

t=rand;

% precompute V and row vector D of eigen values
[V,D]=eig(A,'vector');
D = reshape(D,1,[]);

(V.*exp(D*t))/V

ans = 2×2

1.0225 0.1112 0.0002 1.0966

% which is
expm(A*t)

ans = 2×2

1.0225 0.1112 0.0002 1.0966

If u(t) is "nice", the integration of exponential would also close formula, you wouldn't need to do it numerically.

[Star Strider]

If I remember correctly,

.

[Bruno Luong]

More detail of solving using eigen vectors in here

There is an integral in the solution if u(t) depends on t. But the integral becomes 1D integrals if written in the eigen space of A, IIRC.

[Michal]

OK, I prepared two solution of the same problem:

1. your method (general ... works with any size of A), with vector of times

function Y = expmct(X,t)
Nt = length(t);
% precompute V and row vector D of eigen values
[V,D]=eig(X,'vector');
D = reshape(D,1,[]);
Y = zeros(2,2,Nt);
for i = 1:Nt
    Y(:,:,i) = (V.*exp(D*t(i)))/V;  
end
end

2. my method for special case (2x2) constant matrix A with analytical form of matrix exponential (see here)

 function Y = expmt_22(X,t)
%

% Check input matrix size
if any(size(X) ~= [2,2])
    error('Input matrix X size is not 2x2 !!!')
end

% Matrix X is real or complex (2 x 2)
g = (X(1,1) - X(2,2))/2;
d = (X(1,1) + X(2,2))/2;
D = 4*X(1,2)*X(2,1) + (X(1,1) - X(2,2))^2;

% first auxilliary term
aux1 = exp(d.*t);
if any(isinf(aux1))
    warning('exp overflow ... indefinite result')
end

% D special cases
if (D ~= 0 && ~isreal(X)) || D > 0
    D  = sqrt(D)/2;
    Dt = D.*t;
    s_Dt = sinh(Dt);
    c_Dt = cosh(Dt);
    if any(isinf(s_Dt)) || any(isinf(c_Dt))
        warning('sinh or cosh overflow ...  indefinite result')
    end
elseif D < 0 && isreal(X)
    D  = sqrt(abs(D))/2;
    Dt = D.*t;
    s_Dt = sin(Dt);
    c_Dt = cos(Dt);
else % D == 0 case
    D = 1;
    s_Dt = 1;
    c_Dt = 1;
end

% auxiliaries
aux2 = aux1.*s_Dt./D;
aux3 = g.*aux2;
aux4 = aux1.*c_Dt;

% output 
Y = zeros(2,2,length(t));
Y(1,1,:) = aux3 + aux4;
Y(1,2,:) = X(1,2).*aux2;
Y(2,1,:) = X(2,1).*aux2;
Y(2,2,:) = aux4 - aux3;

end

The finel test:

t = 0:0.0001:1;
X = rand(2);

tic;Y  = expmct(X,t);toc
Y_ = expmt_22(X,t(1)); % warming up
tic;Y  = expmt_22(X,t);toc % run

max(abs(Y-Y_),[],'all')

The problem is how to perform effectively the integration of exponentials (Non-homogeneous term) in a case of time vector t = t0:dt:t1? How to effectively formulate this task via proper vectorized code with multi-arrays?

[Michal]

@Bruno Luong What exactly you mean by: "If u(t) is "nice", the integration of exponential would also close formula, you wouldn't need to do it numerically."

Cloud you clarify this specific task? u(t) are noised measurements so it is difficult to say if they are "nice" or not :)

[Bruno Luong]

This is the integration term of your solution

A=rand(2);
t = rand;
Btu=rand(2,1); % <= EDIT; B(t)*u(t)

[V,D]=eig(A,'vector');
D = reshape(D,1,[]);

% Numerical integration term Integral(expm(A*(t-s))*B*u(t)*ds)
I =  integral(@(s) expm(A*(t-s))*Btu, 0, t, 'ArrayValued',true)

I = 2×1

0.2050 0.3494

% Eigen space, direct formula (no integration needed)
I =   (V.*((exp(D*t)-1)./D)) * (V \ Btu)

I = 2×1

0.2050 0.3494

[Bruno Luong]

Noise are not nice. You have stochastic ode and not ordinary ode. I'm not expert on stochatic ode, but it seems to me that you have to pay attention to integration of such "noise" and what represents when you discretize the noise.But don't ask me more question, it's not my domain of expertise.

[Michal]

@Bruno Luong integration term should be a vector not a matrix.

Btu (vector) = B*u (matrix * vector)

[Bruno Luong]

B is the matrix no? If it's a vector then replace B with diag(Btu) in the last formula of my code.

EDIT, simply enter Btu as (2 x 1) vector.

[Michal]

Ok

Did you mention that my method expmt_22 is significantly faster then a metod you propose (eigen space)? Is possible to optimize your method even further? How to vectorized your method for time vector to compute all time steps of complete solution more effectively?

[Bruno Luong]

I don't know I did not read your code, and I don't know how exactly you implement mine (it can be fully vectorized (the matrix multiplication would be performed outside, NOT repeatly inside on vectors). I won't involve in the implementation.

[Torsten]

@Michal Kvasnicka

All your methods don't include the B*u-term with discrete u-values. But this will be the crucial point with the "analytical" solution method.

[Michal]

Could you please perform the final test run (see my previous post). I will be really happy, if you could be so kind to help me find optimal solution of my problem.

[Michal]

@Torsten In what way? So, is there any way how to solve my problem for discrete signal u(t). Of course, the sampling period is small as possible, so I hope that the any kind of noise does not disturb the final solution significantly. You propose any kind of pre-smoothing of signal u(t)?

[Torsten]

You propose any kind of pre-smoothing of signal u(t)?

I have no experience in this field, so I don't know. As Bruno said, it's also difficult to decide whether the usual integration methods for ODEs make sense in your case as we don't know the strength of the noise component. But I have a feeling that strict "analytical" integration and noise somehow contradict each other.

[Michal]

@Torsten Ok... Last question: Application of ode15s or ode45 is more reliable in a case of discrete signal u(t) with noise? I am afraid that the problem is still the same.

[Torsten]

Application of ode15s or ode45 is more reliable in a case of discrete signal u(t) with noise?

No, but (at least for me) the ODE integrators are easier to handle than the analytical solution formula.

[Bruno Luong]

IMO using ode solver are all wrong approach.

One should rather use sde

[Torsten]

Don't know how such problems are handled in signal processing and control engineering.

Maybe Star Strider has a better insight than we have.

[Bruno Luong]

I don't know the suitable numerical technique, however I know how the integral term ends up to be if u(t) is a Guassian process. IMO none of the ode solver would provide the correct solution that I expect.

[Star Strider]

It just looks like a linear control problem to me, at least as originally described (even though the C and D matrices are not defined) . That is the only reason I suggested that approach.

[Michal]

Bruno is very probably right. ODE solvers of any kind are not suitable tools for this kind of problem with discrete noised control signal u (t). The SDE approach is, of course, much more difficult, but fully adequate. So far, I have no idea how to translate my relatively simple problem into the correct SDE formulation. In addition, any performance requirements are now out of interest until the problem is reformulated into SDE.

Thank you for clarifying my problem.

[Michal]

@Bruno Luong Bruno, I just re-check your integration terms

I =  integral(@(s) expm(A*(t-s))*Btu, 0, t, 'ArrayValued',true)
I =   (V.*((exp(D*t)-1)./D)) * (V \ Btu)

I am afraid, that are not correct, because But should be integrated over s. In your codes is this term constant?!

Right formulation is:

time = [t1,t2, ... , tN]; % sorted samples time vector
N = length(time);
t = (time(end)-time(1))*rand + time(1); % single random time sample  (t1 <= t <= tN)
u = [u1,u2, ...,uN]; % u(t) samples    
I = integral(@(s) expm(A*(t-s))*B*interp1(time,u,s), 0, t, 'ArrayValued',true);

What impact has this fact on your eigenspace formulation?

[Bruno Luong]

@Michal Kvasnicka Sorry firstly because YOU wrote above (comment just bellow Star answer) u(t) and NOT u(s) inside the integration:

x(t) = expm(A*(t-t0))*x0 + Integral(expm(A*(t-s))*B*u(t)*ds)

[Michal]

@Bruno Luong Sorry ...my mistake (typing error)!!! General solution of linear ODE system

x'(t) = A*x(t) + B*u(t)
x(t0) = x0

looks like:

x(t) = expm(A*(t-t0))*x0 + Integral(expm(A*(t-s))*B*u(s)*ds)

where integration is over t0 <= s <= t

[Bruno Luong]

Change are not large if integram over u(s), the integral term

I = Integral(expm(A*(t-s))*B*u(s)*ds)

can be computes as

I(t) = (V.*((exp(D*t)-1)./D)) * (V \ B) * U(t),

where U(t) = integral on (0,t) of u(s) ds.

You simply integrate u() using MATLAB trapz or cumtrapz or whatever. The rest of the term can be computed efficiently.

[Michal]

@Bruno Luong See the following test code:

A=rand(2); % homogeneous matrix A
t=(0:3); % time vector
Nt = length(t); % number of samples
u = rand(1,Nt); % signal vector
B=rand(2,1); % non-homogeneous vector B 

[V,D]=eig(A,'vector');
D = reshape(D,1,[]);

I_ = zeros(2,Nt);
for i = 2:Nt
    % standard numerical integration
    I_(:,i) = integral(@(s) expm(A*(t(i)-s))*B*interp1(t,u,s), t(1), t(i), 'ArrayValued',true);
end
I = zeros(2,Nt);
Ut = cumtrapz(t,u);
for i = 2:Nt
    % Eigen space approach
    I(:,i) = (V.*((exp(D*t(i))-1)./D)) * (V \ B) * Ut(i);
end

and compare I and I_ array. There are huge differences! Standard numerical integration is my reference method.

Any idea what is wrong?

I am not sure, that is possible in principle to exclude expm(A(t-s)) from the numerical integration by any "trick".

[Bruno Luong]

Sorry I'm not going to debug, don't have time?

[Michal]

I have a plenty of time :), but I still do not fully understand if is your eigenspace approach mathematically correct or not. I am trying to follow all your recommendation, but so far, does not work

[Bruno Luong]

You are right the integral can not be separate the exponetial term and the u term.

However if you interpolate u(s) term as linear between 2 adjacent poins (as you do with you reference code), in analytic way we end up to compute the two atomic integral,

integral exp(-d*s) ds
integral exp(-d*s) s ds

Those have close formula. Then you combine on each interval using the linear coefficients of u(s)

[Michal]

Yes, this is the right approach how to solve analytically my ODE system. Linear 2-point interpolation is completely sufficient. Thanks

But still is open general question if it is mathematically correct method for noisy signals u.

[Bruno Luong]

It all depends on what you want to simulate: if you want to simulate stochastic process, no, if you have noisy data where you want to simulate with minimal noise impact, you have to study the sensitivity of the your system to noise you have to do some suitable filtering or regularization. Only you know what is the goal. The ODE with linear system and damping behavior should natuarally regularize the noise. But who knows about the characteristic of your A and B matrices (the real part of the eigenvalues reveal the damping characteristic of your system).

[Michal]

@Bruno Luong Just a last question regarding application of atomic integrals on 2-point linear segments of u(s).

If I uderstand well your recomendation, you propose the following algorithm:

1. each linear segment of u(s) is approximated as line a_i*s+b_i (where a_i - slope and b_i - intercept of i-th linear segment) on interval [0,dt], where dt is time step.
2. the integral term is then I = Integral (expm(A(t-s))*B*u(s) ds), where is integrated over [0,t] interval
3. at eigen space is this integral in the form: a_i*Integral(exp(D*(t-s))*s ds) + b_i*Integral(exp(D*(t-s)) ds), where is integrated over [0,dt] interval correponding to each linear segment (in this case intercept b_i equal to u(s) at 1st point of linear segment, and slope a_i is difference u(s) at end - start point of i-th linear segment).
4. in the close form: A_i = a_i*(exp(D*dt)-D*dt-1)/D^2 + b_i*(exp(D*dt)-1)/D)
5. Fot each i-th consecutive linear segment of signal u(s) we have increment A_i
6. The final step is then I = V.*sum(A_i)/V * B

Am I right? So far I am not able to create implementation which correspond to the reference solution

I = integral(@(s) expm(A*(t-s))*B*interp1(t,u,s), 0, t, 'ArrayValued',true)

[Bruno Luong]

I don't think it's completely right the sum must involve exp(ti) somewhere I don't see it in your formula, where the integral needs do be written for each segment (ti,ti+1), upi might change ine variable of integration s to ri := s-ti for each inteval to avoid a big value of intercept at 0.

I don't have time to give you the exact formula and algorithm. As you have plenty of time you just need to do step by step and check the formula.

[Michal]

Everything solved and functional ... Plenty of time is good for solving of any problem :)