Table of Contents

Numerical Integration of Differential Equations

Most models are on a deeper level described as differential equations. The reason is that most physical phenomena are described as variations of state variables or other quantities. For example, the derivatives in time of voltages can characterize an electronic circuit completely. The derivatives in time and place of concentrations are the source of modelling particle transport phenomena. A constellation of phenomena descriptions (component equations) and associated relations (network equations) result into differential equations. If we can not solve the equation to an explicit form analytically, there are basic methods to simulate a response curve. The essence of this method is the Newton integration.

Differential equations

Almost every engineering problem / model can be stated in a differential equation. For a free falling body we find for example that the forces on the body are gravitation (mass $M$ times gravitational constant $g$), and the friction which is a function of velocity speed $v$. As a result, the experienced force on the body (equal to acceleration a times the mass $M$) can implicitly be expressed as

Now it is a minor step to express the position $x$ of the body, by realizing that acceleration a equals $\frac{\partial^2 x}{\partial t^2}$ and velocity $v$ equals $\frac{\partial x}{\partial t}$. We find the differential equation

With this equation we can try to find a solution for $x$ as a function of time. What is needed, however, are boundary conditions. What are the velocity $v$ and position $x$ at $t = 0$? Only then, the solution is defined.

Another example is an electronic RC circuit that is charged as shown in figure 1. The input voltage is the sum of the voltage over the resistor and the capacitor, the output voltage the voltage over the capacitor.

This one essentially differs from the free falling body because there is an input signal. In this case it is not sufficient to know the voltage over the capacitor $U_C$ as a function of time only. We also need the signal of $U_{in}$, which can be a signal in time as well.

Euler integration

For readability, we will now use a different way of writing:
\begin{equation}
u'=\frac{\partial u}{\partial t},
u''=\frac{\partial^2 u}{\partial t^2}.
\end{equation}

Consider having a system in a state where it has a known output $u \left ( t_{0} \right )$. If we make a small step $\Delta t$ in time, and we want to know the corresponding output $u \left ( t_{0}+\Delta t \right )$, we could make the approximation

which is in fact the first order Taylor series expansion. What does this equation mean? In fact, it approximates the next value $u \left ( t_{0}+\Delta t \right )$ by using the known current value $u \left ( t_{0} \right )$ and drawing a straight line to the supposed next point. The straight line is constructed from the slope times the step size $\Delta t$. The slope is found from the gradient, which is the first derivative in time of $u \left ( t_{0} \right )$ indicated as $u' \left ( t_{0} \right )$. So, to calculate the next value, we need $u \left ( t_{0} \right )$ and $u' \left ( t_{0} \right )$ only, because it is a linear approximation. This is indicated in figure 2.

Fig. 2: The principle of Euler integration along a curve u(t)

Once knowing $u \left ( t_{0}+\Delta t \right )$, we can calculate $u \left ( t_{0}+\Delta t+\Delta t \right )$ in the same way starting at $u \left ( t_{0}+\Delta t \right )$. This iterative method can be generalized by

which is called first order Euler integration, or the first order Euler method.

What is needed to apply this method for reconstructing a curve? First of all, we have to take a step size $\Delta t$ (not necessarily equal for all steps) which is small enough so that we do not loose track of the real curve. Next, we have to know the initial position (or the start criteria) and the values at any $t_{n}$ of all independent varables. The start criteria and independent variables are called the boundary conditions. Finally, we must be able to know the first derivative at any time in order to calculate the next step. However, when we have a model in differential equation form, we know the systemic behavior, and we do know the first derivative.

An example of how to express the first derivative is based on the RC-circuit. The differential equation was given in equation \eqref{eq:RCDiff}. We can re-write

Higher order Euler integration

Another example is based on the free falling body which is modelled with the differential equation \eqref{eq:FreeFallingDiff}. At first sight, we may think we have a problem because there is a second derivative included. The mathematical solution is in breaking the equation into multiple steps. Start with the substitution

Analysis in Simulink

The Euler integration of differential equations \eqref{eq:FreeFallingDiff} and \eqref{eq:RCDiff} can be done in block-wise simulators like Simulink. In that case, we don't have to separate the orders as described in the previous section. To convert a differential equation to a block model, there are two methods.

Bringing the highest order to the left

The generic form of a differential equation is
\begin{equation}
A_{3}\frac{\partial^3 u_{out}}{\partial t^3}+A_{2}\frac{\partial^2 u_{out}}{\partial t^2}+A_{1}\frac{\partial u_{out}}{\partial t}+A_{0}u_{out}+C=B_{0} u_{in}
\label{eq:GenericDiff}
\end{equation}
which is an example of a third order system. Write the differential equation as
\begin{equation}
\frac{\partial^3 u_{out}}{\partial t^3}=\frac{B_{0} u_{in}-A_{2}\frac{\partial^2 u_{out}}{\partial t^2}-A_{1}\frac{\partial u_{out}}{\partial t}-A_{0}u_{out}-C }{A_{3}}
\label{eq:GenericDiff2}
\end{equation}

with the highest derivative at the left hand side. In this case, that is the third derivative. Now the core of the block scheme can be drawn as shown in figure 3. Each integrator represents one order of integration. The signal path contains $u_{out}$ at the right hand side, and from right to left there is $\frac{\partial u_{out}}{\partial t}$, $\frac{\partial^2 u_{out}}{\partial t^2}$ and $\frac{\partial^3 u_{out}}{\partial t^3}$ respectively.

Fig. 3: Building the Simulink model step 1

Realize that equation \eqref{eq:GenericDiff2} tells us exactly what $\frac{\partial^3 u_{out}}{\partial t^3}$ is. It is the sum of five terms with a division by $A_{3}$. The first term of the addition is $B_{0}u_{in}$. This is added in figure 4 to the Simulink model.

Fig. 4: Building the Simulink model step 2

Finally, we add the other terms as done in figure 5. This is the model that can be used: $u_{in}$ is an input and $u_{out}$ is an output as it should be.

Fig. 5: Building the Simulink model step 3

Directly from the system

An alternative is to convert the system to the simulink model directly from the component equations and the network equations. For example, take the RLC network of figure 6.

Fig. 6: An RLC network

We can see the components as elements that convert a current into a voltage, or vice versa. This is

\begin{equation}
u=Ri \leftrightarrow i=\frac{1}{R}u
\label{eq:Resistor_iu}
\end{equation}
for a resistor,
\begin{equation}
u=\frac{1}{C}\int i dt \leftrightarrow i=C \frac{\partial u}{\partial t}
\label{eq:Capacitor_iu}
\end{equation}
for a capacitor, and
\begin{equation}
u=L \frac{\partial i}{\partial t} \leftrightarrow i=\frac{1}{L}\int u dt
\label{eq:Inductor_iu}
\end{equation}
for an inductor. Of the last two, we prefer the integral versions because they can be translated directly into Simulink operations according to figure 7.

\begin{equation}
i_{C}=i_{L}+i_{R}.
\label{eq:Kirchhoff4}
\end{equation}
Equations \eqref{eq:Resistor_iu}, \eqref{eq:Capacitor_iu}, \eqref{eq:Inductor_iu}, \eqref{eq:Kirchhoff1}, \eqref{eq:Kirchhoff2}, \eqref{eq:Kirchhoff3}, and \eqref{eq:Kirchhoff4} can directly be copied into a Simulink model. First start with the notion that $u_{out}=u_{C}$ from \eqref{eq:Kirchhoff1} and obtain the $u_{C}$ from the capacitor component equation \eqref{eq:Capacitor_iu}. This is converted to Simulink in figure 8.

Fig. 8: Simulink partial model for the RLC circuit

The next unknown is the current through the capacitor. Network equation \eqref{eq:Kirchhoff4} shows it is the sum of $i_{L}$ and $i_{R}$, so the outputs of the component equations \eqref{eq:Resistor_iu} and \eqref{eq:Capacitor_iu}. Now the Simulink model has grown towards figure 9.

Fig. 9: Simulink partial model for the RLC circuit

Because $u_{L}=u_{R}$ according to the network equation \eqref{eq:Kirchhoff2}, we can simply draw a line in figure 10.

Fig. 10: Simulink partial model for the RLC circuit

This new line represents the voltage over either the resistor or the capacitor, so network equation \eqref{eq:Kirchhoff3} closes the loop. This means that all network and component equations are used now, and figure 11 is the finished model.

Fig. 11: Simulink partial model for the RLC circuit

Note that this second method gives a different form in Simulink. If we would have used the previous method with the highest derivative to the left, then the two integrators would be connected directly together without a multiplier in between. As a result, the multipliers will be at different locations. When using the first method on the electrical circuit of figure 6, we would probably have found a model with combined multipliers $LC$ and $RC$, so the capacitor value C would occur two times. It is the most intuitive when all physical components (R's, C's) occur only once in the Simulink model.

Example: A heat transfer system

Consider a cooling plate like figure 12 on an active heat dissipator. Such a heat dissipator can be a high speed processor for example. We want to calculate/simulate the temperature behaviour inside the dissipator, because that disspitor should not over-heat when it is a semiconductor device.

Fig. 12: Cooling plate

The power developed in the dissipator is $P_{in}$, the power radiated to the environment is $P_{out}$. As a result, the heat accumulation $\Delta Q$ per period $\Delta t$ is equal to
\begin{equation}
\Delta Q=\left ( P_{in} - P_{out}\right ) \Delta t .
\label{eq:HeatAccumulation}
\end{equation}
This is indicated in figure 13. The dissipated power $P_{in}$ is normally known because it is an active source.

Fig. 13: Heat transfer model of a dissipator with a cooling plate

The heat accumulation $\Delta Q$ results into an increase in temperature $\Delta T$ with a proportionality factor $C$ like
\begin{equation}
C = \frac{\Delta Q}{\Delta T}
\label{eq:ChangeTemperature}
\end{equation}
with $C$ a constant that is defined by $c_{m}m$ with heat capacity $c_{m}$ and $m$ the mass of the material. Combination of \eqref{eq:HeatAccumulation} and \eqref{eq:ChangeTemperature} gives the differential equation
\begin{equation}
C \frac{\partial T}{\partial t}=P_{in}-P_{out}
\label{eq:ThermalDifferential1}
\end{equation}
which relates the input power to the temperature as a function of time. Only the output power $P_{out}$ is unknown in this equation, but can be found by defining a thermal resistance $R_{th}$:
\begin{equation}
P_{out}=\frac{T-T_{o}}{R_{th}}
\label{eq:ThermalResistance}
\end{equation}
with $T_{o}$ the outside temperature. This equation holds for $T \approx T_{0}$. Differential equation \eqref{eq:ThermalDifferential1} now becomes:
\begin{equation}
C \frac{\partial T}{\partial t}=P_{in}-\frac{T-T_{o}}{R_{th}}.
\label{eq:ThermalDifferential2}
\end{equation}
So this is our model for the temperature $T$ of a dissipator, dissipating the amount of power $P_{in}$ in an environment of outside temperature $T_{o}$. The constants $R_{th}$ and $C$ are determined by the material and geometry.

However, in most cases we can not find an analytical solution because the model is too complex, or the input signal is not a simple function. In that case we may use a numerical integration method, for example in Simulink. Since we have the differential equation already, the method of bringing the highest order to the left from the previous paragraph is the most appropriate.

Start with \eqref{eq:ThermalDifferential2} as shown in figure 3: draw the base line with the temperature $T$ and its first derivative $\partial T/ \partial t$ first.

This model can be simulated for a step response by setting the initial condition to $T(t=0)=T_{o}$ in the integrator block and your will see that the response has the same exponential shape as the analytical solution \eqref{eq:AnalyticalRC}.

For those with an electrical engineering background, you may wonder why the thermal system has a similar response as an electrical RC network. In fact, we have a thermal resistance $R_{th}$ and a heat capacitance $C$ which insinuates “RC behaviour”. The equivalence is indeed present. If a thermal resistance is seen as the equivalence of an electrical resistor, and heat capacitance is seen as the equivalence of an electrical capacitor, then the similarity is completed by taking temperature equivalent to electrical voltage and power equal to current. It can be understood that Power is equivalent to current, because Power if the flow of heat ($\partial Q / \partial t$). The equivalent model is represented in figure 17 by an LTSpice screenshot. Mathematiclaly seen, this is the same as the Simulink model of figure 6 and differential equation \eqref{eq:ThermalDifferential3}. The theory behind physical analogies is explained on the page with Lumped element models.

Fig. 17: Equivalent model of the heater system expressed as an electric circuit