Most effective way to solve nonhomogeneous linear ODE problem

조회 수: 18(최근 30일)
Michal
Michal 2022년 5월 10일
편집: Michal 2022년 5월 23일
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
  댓글 수: 2
Michal
Michal 2022년 5월 10일
@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.

댓글을 달려면 로그인하십시오.

답변(2개)

Paul
Paul 2022년 5월 19일
편집: Paul 2022년 5월 19일
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().
  댓글 수: 8
Michal
Michal 2022년 5월 23일
@Paul OK ... my cumulative integration is obviously completely wrong!!! Thanks for help...

댓글을 달려면 로그인하십시오.


Star Strider
Star Strider 2022년 5월 10일
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.
  댓글 수: 38
Michal
Michal 2022년 5월 12일
Everything solved and functional ... Plenty of time is good for solving of any problem :)

댓글을 달려면 로그인하십시오.

제품


릴리스

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by