# Calculation of consistent initial conditions for ODE

3 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))
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.

R2022a

### Community Treasure Hunt

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

Start Hunting!