stratum wrote:Thanks for the explanation. Since I cannot visualize the multidimensional surface the newton's algorithm travels in this case, I can't understand the difference fully, but I guess as long as the algorithm does converge, that does not matter.

I doubt anyone can really visualise the hyper-surface, but there's no good reason to even try.

Essentially the "matrix" is just a set of linear equations. When the equations we want to solve are non-linear, we pick some estimate and find the linear tangent (hyper-)plane for the "true" non-linear surface at that point. If the non-linearities are "separate" you can do this by finding the partial derivatives one dimension at a time, but there is real reason it has to be this way (it's just convenient in many cases). Then we put the tangent plane (or projections of it) back to the equations where the non-linearities where and try again until the estimate no longer changes (which implies the tangents won't change either).

Also when it comes to transient analysis, we can usually just reuse the tangent plane from the previous time-step and skip the initial estimate completely. This generally works wonders when the circuit state is changing "slowly" in comparison to the time-step and in case of an (almost) static circuit it means we don't even need to evaluate the non-linearities at all (which is why you can get "zero iterations" with Halite in some cases). With lots of high frequency stuff going on this might not be the ideal choice, but in that case you have other problems like aliasing to deal with and the right thing to do is to decrease the time-step (ie. increasing the sampling rate).

The real difference (besides the different approach to "initial guess") between my code and the code from Urs is that he iterates on the error residue rather than the values directly. This is mathematically equivalent and potentially a bit better numerically (and in some cases can be profitable since it can save some iterations if numerical precision starts becoming an issue), but adds a little bit of overhead in the main solver, especially if you want to early-out those dimensions that have already converged, so I didn't bother. LU even with just partial pivoting is usually fairly well behaved anyway, so it's not a huge deal either way.

The only truly "weird stuff" that I do is with the capacitors, where rather than solving for the voltage or current over the capacitors, I just plug in an algebraic equation (which makes little sense in terms of the actual circuit) to solve for trapezoidal TDF2 state directly, which then saves a separate integration step, since you can just pick the new value directly (ie. the integration is done as a by-product of solving rest of the system). While the downside of this is that the integration method is now hard-coded into the stamps, I left it in to illustrate how you can plug pretty much arbitrary linear equations into the system and have it solved together with everything else.