Yeesh, Jim, I haven't looked at this newsgroup since roadrunner dropped NNTP support. I've missed the chaos and ignorance... You need to thump some of the empty heads.
Anyway, about charge conservation in circuit simulation: The problem was always with MOSFETs. BJTs and diodes were fine from the original SPICE1 code -- if you start with a q(v) equation and differentiate it to get capacitance, things are surprisingly good. Path integrals work out... The problem was the Meyer model for channel capacitance, used in the MOS level 1, 2, 3 models (and others). The model described the capacitances and did a numerical estimation of charge based on them. This was just plain wrong, since they're trying to estimate, e.g., qg (vd,vg,vs,vb) without an integrable function. A good exposition of the problem and first practical fix was in Ping Yang's 1983 paper.
So, MOS levels 1, 2, 3 are probably still wrong in SPICE -- I haven't looked lately. Any of the more recent models, BSIM3, BSIM4, EKV, PSP, all start with q(v) and differentiate to get the (non-reciprocal) capacitances, so they're more-or-less correctly set up for charge conservation by construction. For transient analysis, even if they get the differentiation wrong, as long as q(v) is correct the model will conserve charge -- if it converges.
For the picky, yes there will be some error on each cycle due to convergence tolerances. If the simulator you're using properly implements truncation error timestep control, the loss will be minimal and controllable by tightening reltol and chgtol. If you're using certain 'gold standard' simulators, you'll probably have trouble unless you're some kind of prescient in figuring the right model switches and options. My simulator does quite well.:)
--Steve