# Most effective way to solve nonhomogeneous linear ODE problem

64 views (last 30 days)
Michal on 10 May 2022
Edited: Michal on 23 May 2022
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 on 10 May 2022
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 on 10 May 2022
@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.

Paul on 19 May 2022
Edited: Paul on 19 May 2022
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().
Paul on 20 May 2022
Edited: Paul on 21 May 2022
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 on 23 May 2022
Edited: Michal on 23 May 2022
@Paul OK ... my cumulative integration is obviously completely wrong!!! Thanks for help...

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2022a

### Community Treasure Hunt

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

Start Hunting!