# Calculation of consistent initial conditions for ODE

2 views (last 30 days)
Fabio Freschi on 2 Jun 2022
Commented: Fabio Freschi on 7 Jun 2022
(minor edit)
Hi everyone,
I have written wrapper functions calling MATLAB's ODE solvers. Specifically, my problems are DAEs with large and sparse matrices, similar to matrices coming from finite element analysis. For this reason I use ode15s, ode23t and sometimes ode15i.
When calculating consistent initial conditions for ode15i using decic, everything is fine with small problems, but for even medium problems decic takes ages and may run out of memory. Profiling the code and looking at the implementation of decic I have found that the bottleneck is obviously the QR factorization. However I also have found that the matrix passed to the qr function is converted to full, as shown in the snippet below
function [rank,Q,R,E] = qrank(A)
% QR decomposition with column pivoting and rank determination.
% Have to make A full so that pivoting can be done.
[Q,R,E] = qr(full(A))
The reason of the conversion is explained, however, removing the full function, the code still works.
Is this a legacy option that now is no more required? I don't like to make changes to the built-in files, so in case I would locally copy the dicic function, however, if I remove the full command, am I doing something dangerous?
Below you can find a code that can be used to test decic
% dimensions
nx = 25;
n = 3*nx^2+4*nx+1;
% mass matrix
M = gallery('wathen',nx,nx);
% stiffness matrix
K = gallery('tridiag',n);
% rhs
b = rand(n,1);
% problem: f(t,x,dx) = b(t)-M*dx(t)-K*x(t)
f = @(t,x,dx)b-M*dx-K*x;
% jacobian
options = odeset('Jacobian',{-K,-M});
% consistent initial conditions
x0 = zeros(n,1);
dx0 = zeros(n,1);
[x0,dx0] = decic(f,0,x0,[],dx0,[],options);
% check
norm(f(0,x0,dx0))
Christine Tobler on 7 Jun 2022
@Fabio Freschi Quick response on QR and R2017a vs. now:
"Checking the documentation of qr for release 2017a (here), I can find the following
[Q,R,E] = qr(A)
If A is full: [...] The column permutation E is chosen so that abs(diag(R)) is decreasing.
If A is sparse: [...] The column permutation E is chosen to reduce fill-in in R.
There is no such a distinction in the current documentation (here)"
The change in the more recent documentation is purely one of documentation; on rewriting that page in the more modern template, we had preferred making the doc page simpler by removing those additions. I will bring up the possibility of changing this back.
The behavior is still the same as in R2017a: In dense column-permuted QR, you can rely on abs(diag(R)) being decreasing. In the sparse case, that is not the true. However, usually low-rank matrices are detected there in the sense that columns that have norm near some small tolerance are usually dropped completely; this kind of behavior is used in backslash to get an estimated rank of a sparse matrix.
To the larger point of this post, we have received a request to improve performance of decic for the sparse case, and will look into it - thank you for pointing out this issue.
Fabio Freschi on 7 Jun 2022
@Christine Tobler: Thanks for the comment. It helps in the understanding of the behavior of QR.
The request you mentioned is mine!

Bjorn Gustavsson on 3 Jun 2022
@Fabio Freschi, OK, now we're at least in the same general region of worrying. If I read decic correctly the use of the QR-output is in lines like these:
d = - Q' * res;
dy = E * (R \ d);
Where E is the permutation-matrix. To the best of my understanding the QR-decomposition is a well-defined matrix decomposition, and it would be "very peculiar" if we would get a completely different solution of a linear system of equations if we used
[Q,R,E] = qr(A);
compared to
[Q,R,E] = qr(full(A));
The difference in the permutation might/will give normal variations on the level of floating-point accuracy, but surely nothing worse than that (which is very easy for me to say since you're the one working with the time-bomb). If you have a perticularly ill-posed problem that might be a problem, but that would be a difficult problem anyways.

### Categories

Find more on Programming 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!