[Qucs-devel] Offer to help concerning the numerical difficulties

[Soeren D. S. <soe...@gm...>]

Hi everyone,

I'm a student of Industrial Mathematics at the University of Bremen. 
I've always liked Qucs, but I've noticed the convergence problems in 
transient simulation.  Those are one reason why I'm currently writing my 
Bachelor's thesis on exactly this subject.  Due to months of study of 
the numerical mathematics involved, I think I now have quite some 
understanding of what's going wrong in Qucs.

Last summer, I played a bit with the source code of Qucs.  This is a 
quick comparison of my insights in the implementation and the 
theoretical knowledge that I've developed in the meantime.

In general, the problems that arise in circuit simulation are called 
"Differential-Algebraic Equations" (DAEs).  Those are equations which 
look like:

M dx/dt = f(t,x)

with a possibly singular matrix M.
Compared to Ordinary Differential Equations (ODEs)

dx/dt = f(t,x),

they expose some very nasty numerical behavior, and Qucs doesn't really 
address this.

1. First of all, M and f are never really accessible in the source code. 
  This makes it difficult to apply different numerical methods.

2. I realize that multi-step methods are very popular, but concerning 
DAEs, they aren't always the best choice.  The problem about DAEs is the 
decline of the condition number of the Jacobian with small stepsize. 
Upon stability problems, multi-step methods often require low order and 
small stepsize, which may be fatal in the case of circuit simulation. 
The "fallback" method "Implicit Euler and minimum stepsize" is about the 
worst thing you can do.  There are equations on which the Implicit Euler 
method doesn't even converge in theory, and even less so in practice. 
To me, implicit one-step methods are far more promising because they 
show excellent stability even at arbitrarily high order.  They're harder 
to implement, but good implementations have been in existence since the 
1990s (most prominently radau5.f).

3. The implementation Newton's method is *very* unfortunate because it 
effectively factors out the Jacobian, squaring the condition of the 
problem.  As the condition number is often bad enough anyway (see 
above), this can be fatal, too.  It's very laudable that the Jacobian is 
calculated analytically (though not necessary with an otherwise good 
implementation), but running a new decomposition on every Newton 
iteration is just a waste of calculation time.

4. In the context of my Bachelor's thesis, I have developed a method 
which improves the condition problems, sometimes fixing them entirely. 
I have proven the theory about it, and now I'm just personally curious 
how well it performs in practice.  The methods employs Singular Value 
Decomposition, which is already implemented in Qucs, but it applies it 
an a different way.


Unfortunately, the general numerical approach is quite hardcoded in 
Qucs, which makes it non-trivial to improve anything.  As a first step, 
(1) has to be fixed, so it will be possible at all to link Qucs to some 
pre-existing code.  I see that Qucs implements a huge number of methods 
on its own, but I don't really understand the reason, as there exist 
highly optimized implementations in Netlib and similar libraries.

I'm currently still busy with my thesis and some other things, but I 
will have time in August and September.  Is there anyone who currently 
maintains the back-end?  I'd rather communicate my changes to someone 
who knows the current code and who will be around for some time in the 
future (which I can't promise for myself once the next semester starts).


Sören