# Trapezoidal approximation of inductor current from moment to moment

I'm working on a circuit simulation program.

I've gotten passive circuits to work, even ones that have capacitors and inductors, but the simulation has to have a very small timestep when inductors are present.

I've found experimentally that when capacitors are present the largest possible timestep is x10 the value of the smallest capacitor (although I wouldn't mind slightly larger timesteps, on the whole, the timestep scale required by capacitors isn't too bad). For example, if I have a 1uF capacitor, the greatest timestep to produce reliable results is 10us, which is reasonable.

But it's much worse for inductors - to produce reliable results the timestep has to be as low as x1000 the lowest inductor. So for a 1 uH inductor, the timestep can only be 1 ns.

The issue is that I'm not using any approximations, but only the instantaneous formula. For example, the current at time t2 for an inductor is computed by:

I(t2) = I(t1) + ((V(t1) / L) * (t2 - t1))

It works nicely, but requires a very small timestep unfortunately.

I've found an article that explains how to use a trapezoidal approximation so you could reliably estimate for relatively larger timesteps (where applying the instantaneous formula would produce exaggerated results): http://circsimproj.blogspot.com/2009/07/companion-models.html

I'm rather daft in the abstract, so I'd really appreciate if you could show me how to apply what the article shows to one concrete example, and I can generalise from there.

• Let's say the current at t1 through an inductor is 1 mA, and the voltage 1 V. The inductor is 1 uH.

• A later time t2 is: t2 = t1 + 1 us (a microsecond later).

• I'd like to use the simplest possible way (through trapezoidal approximation or something similar to what the article above shows) to find I(t2) or the delta between I(t1) to I(t2).

• I don't know the voltage at t2 (I haven't solved the circuit for t2 yet, only for t1, and from my solution at t1, I need to find the inductor current at t2).

A step by step solution to this example would be amazingly helpful (because in honesty I can a little daft with abstraction).

Please note that all I know is:

1. The inductance of the inductor
2. The current of the inductor at t1
3. The voltage of the inductor at t1
4. The time gap between t1 and t2

And from these I'd like an approximation of the current at t2 (I don't care about the voltage at t2 - that I will find when I simulate the whole circuit at t2). The approximation has to be better than the basic formula that I already use (namely: I(t2) = I(t1) + ((V(t1) / L) * (t2 - t1)))

Even now it works surprisingly well, usually the difference with the Falstad circuit simulator is within 1 mV (when the input voltages sum to 2 V). But the speed of simulation is much slower, in part because I'm forced to use very small time steps..

Thanks in advance for figuring this out for me! (Only if you could please make the math as concrete as possible so I'd have a recipe in one example, which I can then generalise to other initial current \ initial voltage \ inductance values).

My simulator:

The Falstad circuit simulator (which is really great for education by the way):

Edit: The article linked above may not be as informative as I'd like, but I'm quite sure the trapezoidal rule is the way to go because Paul Falstad uses it in his simulator (and he's a true genius).

Conclusion: read Numerical Methods for Engineers by Steven Chapra (Jonk suggested it). Looks like a particularly well written book that matches my level. In the meantime, I'm glad that the simulation works very nicely when I reduce the timestep enough.

• Are you trying to solve "delta" "recurrent" equations? Commented Oct 16, 2022 at 6:46
• I did this in my youngest time when simulators did not "exist" or were not too cheap. If I remember well, one should write as many equations as the number of reactive components. In the example, there is an inductor and a capacitor, thus one must write 2 equations (one with the current inductor increment - delta I, and one with voltage increment of capacitor - delta Vc). Solve two equations to find the delta IL and the delta Vc. I did this for solving buck smps and others. Will try to "recover" these files ... Commented Oct 16, 2022 at 7:05
• What I did to find the "next" point is a bit complicated ... because it is an "estimated value" calculated with the "n" last values of the variable which are used to calculate the nearest "polynomial" function. Commented Oct 16, 2022 at 7:21
• EE&O. Changing the time step "automatically" is used in the simulator (in mine). It is called the auto-adaptative step (however max step time is fixed by default in parameters). One can choose also the "fixed" step if needed (used in convergence problems). This is fixed also by some parameters in the simulator. Commented Oct 16, 2022 at 9:36
• @ee_student I just copied this from LTspice's Help page on the topic. It's a professional grade simulator, so there may be some helpful thoughts in it. I completely agree with Hugh's comment that you should get yourself a good numerical methods book. Within it, there will be a couple of chapters you'll want to focus on more closely. I have an earlier edition of Numerical Methods for Engineers, which I've found useful. You also may want to look at Numerical Methods in Engineering and Science: C, C++, and MATLAB.
– jonk
Commented Oct 17, 2022 at 3:36

Your calculation of the current increment during one timestep by multiplying the derivative of the current at t=t1 i.e. V(t1)/L with the timestep t2-t1 is known as Euler's integration method. It's computationally ineffective, because it does not take into the account how the voltage evolves during a timestep. To get in practice accurate enough results you may need very small timestep, as you have already noticed.

Just for information: Mathematicians (I'm not one) have worked hundreds of years to develop methods which give more accurate results with less timesteps. Their job was not only to invent new calculation formulas, but also to give valid estimates how much error they cause. To understand the thing properly you really should get some book of numerical integration and numerical methods to solve differential equations. I'm afraid proper learning needs background knowledge worth one year of university level math studies.

As you already have heard and also got as an answer, the simplest method to take the evolution of the voltage into the account is the trapezoidal method. It simply uses the average (V(t2)+V(t1))/2 instead of V(t1) to calculate the current increment during time interval t1...t2. That's no problem if the voltage comes from a known source which can be evaluated beforehand for any moment of time.

(BTW. The name trapezoidal comes from the geometric interpretation. To get the new increment one calculates the area of a trapezoid instead of a rectangle.)

In solving circuits (or any state variable differential equations) it happens often that one cannot look one step forward to get the V(t2) because V(t2) can be got only by calculating at first the current at t=t2. If looking one step forward was possible (=known voltage source), it of course would be used, but if V(t2) needs I(t2), you can try to guess V(t2) from the preceding values of V. The older answer suggests you to find the first estimate for V(t2) by calculating at first with the simple Euler method I(t2) and then V(t2) with the circuit theory. Then you use the estimated V(t2) in the trapezoidal rule.

A little simpler idea is to think that V changes during t1...t2 as much as it changed during the previous timestep. If we write V(to) for the V in the beginning of the previous timestep you can calculate the current increment as

I(t2)-I(t1) = ((3V(t1)/2 - V(to)/2)/L)(t2-t1)

Of course, in the starting point there's no previous calculated V(t), but you can calculate the first step with Euler's method and continue with the formula I wrote above. Unfortunately I do not know if it has any name. But I'm sure it helps substantially.

In my own informal numerical tests I have found that to reduce the cumulated error 50% one must shorten the timestep to 50% if the integration is done with Euler's method. By using the formula above shortening the timestep to 70% i.e. by multiplying the timestep with the squareroot of 0,5 gives the same 50% reduction of the cumulated error. To test I integrated exponentially growing function e^t.

As said, my test was only informal tinkering. To get scientifically acceptable error estimates you must bite the bullet and study the subject from some book of numerical integration.

• Thanks kindly for writing out the answer. For me to understand it a bit better, if you could please write, how do you get I(t2) when you only know I(t0) \ V(t0). (ie two iterations of your proposed method in sequence). No worries about plugging in actual numbers, just with params for numbers would be great. Commented Oct 16, 2022 at 20:08
• I'm just thinking a sequence of two timesteps of your method where every operation \ equation is written down would be great. By the way, is the 3 times V(t1) in the last equation intentional? Commented Oct 16, 2022 at 20:16
• It sounds to me that your method is similar to what I wrote in the last set of numbered list points (1-5) at the very bottom of my question (as an edit). Might there be a difference or would these be the same methods? Commented Oct 16, 2022 at 20:25
• I've added a code snippet (to the question). Perhaps this might make my question more concrete. If you get a chance, thanks for any suggestions & such.. Commented Oct 17, 2022 at 5:26

Using your terminology, the trapezoidal rule is:

$$I(t2) = I(t1) + \frac{v(t1) + v(t2)}{2\times L} \times \Delta t$$

Note that in your example, $$\ \Delta t\$$ is the timestep (1 $$\ \mu \$$s).

So the main difference compared to what you already have is that you use the average of V(t1) and V(t2) instead of just using V(t1). I think you may find that sometimes (depending on the circuit) in order to solve for V(t2) you need to first find I(t2) so that could be a problem.

In that case you might have to go through several rounds of approximation until V(t2) converges on a final consistent value compatible with I(t2). But I am not sure.

• No no no. You don't need to know V(t2) to solve I(t1). You need to know V(t2) to solve for I(t2). What you can do is solve for I(t2) using your existing equation, then solve for V(t2), then plug that in to the trapezoidal rule. Then you get an improved estimate of I(t2), then use that to solve for V(t2) again. Etc. You only have to do it until the solution converges and neither V(t2) nor I(t2) are changing appreciably with further iterations. This is NOT my area of expertise. But convergence in simulation is an academic topic. Commented Oct 16, 2022 at 7:51
• Yes. I am not sure if it is better to iterate multiple times at each time step or just have a much smaller timestep. But those seem to be the choices. You are traveling a well-worn path. I don't have a map, but I do know that others have traveled down this path before and learned much about the terrain. Commented Oct 16, 2022 at 7:59
• I wonder if using ngspice as a computation engine for your program might make sense? Commented Oct 16, 2022 at 8:00
• I have wrote the code for some circuit solver. I used Euler's integration method, trapezoidal and Runge Kutta and get similar results. In my view, at the end of the day, it's always a matter of compromise between the time step and convergence. In my case, when the time constant (L/R) gets too small, I need to reduce the time step, but then I increase the step again if the constant increases. I don't keep calculating the constant, I just look at the number of 'inner' steps to get convergence and decide whether to increase or decrease the time step. It's not an answer, but it might help. Commented Oct 16, 2022 at 12:10
• The theory documentation about how ATP (based on EMTP developed at Bonneville Power Administration by Herman Dommel and others) works goes into great detail on doing just this. See here where are you can get a free license and access to that documentation. Commented Oct 16, 2022 at 13:56

As pointed out by @Hugh Harper, I (re)make a Maple sheet showing
the result of recurrence (Euler ?) solving versus really "mathematically" solving it.

A simple case of serial Resistor (1 Ohm) and inductor (1 H), "calculating" current (start of curve).
Here, the number of points is fixed (dt = 0.1s, np=100 -> for 10 seconds).

One can see that recurrence solving (square points) is not the "real curve" (line).
This show clearly the dispersion of the points.
One can see also that, in this case, recurrent points are "above" real points.
What if ... number of points is bigger ... or dt is lower.

EDIT:
I made also a "test" with some kind of ponderation, modifying procedure upd.
Added a line that "interpolates" the value of one point by the two neighbors (2,2,1).

And here is the result (?) ... See only the left half part of the picture.
Error is now negative at starting ... and "positive" after ... then ~ 0.