From: Soeren D. Schulze <soeren.d.schulze@gm...>  20130826 18:29:29

Hi again, I created a new branch called "numerical_improvements" which includes some changes to improve the convergence behavior of transient analysis. Summary of the main changes: * In the corrector and the nonlinear solver, the implementation of Newton's method was changed to an "update" style rather than computing the solution directly. This leads to a number of advantages. First, the achievable accuracy is higher. That's because the updates between the integration steps are rather small, so a specific relative error (which is usually determined by the condition of the Jacobian) has a smaller impact when it's on the update rather than on the solution itself. The program is less likely to experience Newton rejections upon small step size. Additionally, the change enables the use of the socalled "simplified Newton method" which may save some computation time. * Specifically, energystorage components are no longer modeled as "conductance + current source", but their values are directly written into the socalled "F matrix" which is later combined with the A matrix according to the step size. * Source stepping is now used as a convergence helper for the initial DC analysis. All other convergence helpers are temporarily disabled. (It's possible to reenable them, but I currently don't see any need because Newton's method converges locally.) It was a bit surprising to me that source stepping works the way it is implemented, because the code looks erroneous: After a Newton rejection, the rejected guess is used as a starting value for the new iteration. If I "fix" it, I get into a local minimum. If I leave it this way, the method seems to generate "random" guesses which eventually converge after a few 100 steps. (Maybe it's better to use actually random starting values?) * An additional rejection criterion for the Newton iteration was introduced: The iteration is rejected if the update of the solution in one step is greater than the update in the previous Newton step. A similar criterion is thinkable for the RHS, but we probably don't need two. * The Gear predictor is used for *all* integration methods. As I explained earlier, the others can never exceed order 1 because there are no precise values for the derivative of the solution available. This boosts the Trapezoidal method quite a bit, while Euler and Gear are unaffected (they coincide for order 1 anyway). Higherorder Adams methods are always bad, so a better predictor guess doesn't do too much. * e_trsolver was temporarily disabled. What exactly does it do? The changes were done in a conservative way, so old models remain working. Their performance may vary a bit due to some changes in the convergence criteria. On the other hand, almost everything was updated to the new way, except for one problem: * Nonlinear capacitors (and inductors) should to be reconsidered. I fail to find any mathematical convergence analysis for what Qucs currently does. There is some literature on BDF concerning the code DASSL, but the exact algorithm is a bit different. So while the current method may work, I think what saves a lot of trouble is switching to the "chargeoriented MNA", which introduces an extra equation and an extra variable for the charge of each nonlinear capacitor. In this case, the numerical properties are much betterknown, and we can also plot the charge easily if it's interesting to the user. Unfortunately, such a change is not completely transparent, and it requires changes in VerilogAgenerated code. Is there anyone who knows this code well and can estimate if it's easy to add an extra column and an extra row for each nonlinear capacitor? What's left on my TODO list: * Clean up the code. * Document it in the technical papers. * Speed it up. I'm currently following "premature optimization is the root of all evil". My code introduces a bit of CPU time and memory overhead. While matrix decomposition is still slow, it doesn't really matter, but it shouldn't be there. On the other hand, implementing the simplified Newton method might reduce the number of necessary decompositions by a factor of 2 immediately. * Implement the Radau5 onestep method. Just implementing the corrector isn't hard, but making the method work well is. Error estimation is based on an embedded order 4 method, the predictor extrapolates an internal polynomial ("collocation polynomial") from the previous step, and an *efficient* solver has to solve both a realvalued and a complexvalued matrix. This is a bit hard to realize using the current abstractions in "nasolver". (Without this optimization, LU decomposition is slower by a factor of 9.) In the original code, a Gustaffson controller is used for step size control. * Implement the method for improving the condition of the Jacobian that I developed in my Bachelor's thesis. If you want to help, please do some tests (I already tested everything except for the nary inductor myself). Tell me if the code works as I described, or if things stop working for you. As a benchmark problem, I attached the circuit of the "Blameless Amplifier" from "Douglas Self: Audio Power Amplifier Design Handbook". It doesn't work in the old code primarily because no source stepping is used in the Initial DC analysis. With the new code, Euler and Trapezoidal work better, but the difference for Gear is small. Circuits that work well with the new code and badly with the old code (or vice versa) would be very interesting for me. Sören 