User Tools

Site Tools


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

\begin{equation} Ma+\frac{C_{friction}}{M}v=Mg. \label{eq:FreeFalling} \end{equation}

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

\begin{equation} M\frac{\partial^2 x}{\partial t^2}+\frac{C_{friction}}{M}\frac{\partial x}{\partial t}=Mg. \label{eq:FreeFallingDiff} \end{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.

Fig. 1: An electrical RC network

We find

\begin{equation} U_{in}=U_{R}+U_{C}=I_{R}R+U_{C}=I_{C}R+U_{C}=C\frac{\partial U_{C}}{\partial t}R+U_{C} \label{eq:RC} \end{equation}

which can be written as the differential equation

\begin{equation} RC\frac{\partial U_{C}}{\partial t} +U_{C}=U_{in}. \label{eq:RCDiff} \end{equation}

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

\begin{equation} u \left ( t_{0}+\Delta t \right ) = u \left ( t_{0} \right ) + u'\left ( t_{0} \right ) \Delta t \label{eq:EulerSimple} \end{equation}

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

\begin{equation} t_{n+1} = t_{n} + \Delta t \\ u \left ( t_{n+1} \right ) = u \left ( t_{n} \right ) + u' \left ( t_{n} \right ) \Delta t \label{eq:EulerIterative} \end{equation}

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

\begin{equation} U_{C}' \left ( t \right ) = \frac{U_{in}\left ( t \right ) - U_{C}\left ( t \right )}{RC} \label{eq:RCDiff2} \end{equation}

which gives the first derivative. In time-discrete form this becomes \begin{equation} U_{C}' \left ( t_{n} \right ) = \frac{U_{in}\left ( t_{n} \right ) - U_{C}\left ( t_{n} \right )}{RC} \label{eq:RCDiff3} \end{equation}

which looks like an obvious step, but now we can see, in combination with equation \eqref{eq:EulerIterative}, that we can indeed calculate $t_{n+1}$ from $t_{n}$ using

\begin{equation} U_{C} \left ( t_{n+1} \right ) = U_{C} \left ( t_{n} \right ) + U_{C}' \left ( t_{n} \right ) \Delta t . \label{eq:EulerIterativeRC} \end{equation}

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

\begin{equation} z_{1}=x \\ z_{2}=x' \\ z_{3}=x'' \label{eq:substitutes} \end{equation} and subsequently we can write down \begin{equation} z'_{1}=z_{2} \\ z'_{2}=z_{3}. \label{eq:substitutes2} \end{equation}

We know $z_{3}$ from equation \eqref{eq:FreeFallingDiff} and can be written as

\begin{equation} z_{3}=x''\left ( t \right )=g-\frac{C_{friction}}{M^{2}}x'\left ( t \right ). \label{eq:FreeFallingDiff2} \end{equation}

So now, equivalent to equation \eqref{eq:RCDiff3}, we can write the system of first derivatives \eqref{eq:substitutes2}

\begin{equation} z'_{1}=z_{2} \\ z'_{2}=g-\frac{C_{friction}}{M^{2}}z_{2} \label{eq:substitutesModel} \end{equation} and subsequently the model becomes \begin{equation} z_{1} \left ( t_{n+1} \right ) = z_{1}\left ( t_{n} \right ) + z_{1}' \left ( t_{n} \right ) \cdot \Delta t \\ z_{2} \left ( t_{n+1} \right ) = z_{2}\left ( t_{n} \right ) + z_{2}' \left ( t_{n} \right ) \cdot \Delta t \label{eq:EulerIterativeFallingBody} \end{equation} which is the iterative description of the model for the free falling body as first derivatives. When we write it like this \begin{equation} x \left ( t_{n+1} \right ) = x \left ( t_{n} \right ) + z_{2} \left ( t_{n} \right ) \cdot \Delta t \\ z_{2} \left ( t_{n+1} \right ) = z_{2}\left ( t_{n} \right ) + z_{2}' \left ( t_{n} \right ) \cdot \Delta t \label{eq:EulerIterativeFallingBody2} \end{equation} it can be understood better how it will give us the position x at any time, using the helper variable $z_{2}$. This $z_{2}$ has the meaning of velocity (the first derivative of x) and its first derivative needed for the Euler integration can be found from equation \eqref{eq:substitutesModel}.What is needed to start iteration is the start position $x \left (t_{0}\right )$ and the start velocity $z_{2} \left ( t_{0} \right ) = x' \left (t_{0}\right ) = v \left (t_{0}\right )$.

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.

Fig. 7: Simulink partial models for a capacitor and an inductor

Next there are four network laws according to Kirchhoff:

\begin{equation} u_{out}=u_{C} \label{eq:Kirchhoff1} \end{equation}

\begin{equation} u_{L}=u_{R} \label{eq:Kirchhoff2} \end{equation}

\begin{equation} u_{L}=u_{in}-u_{out} \label{eq:Kirchhoff3} \end{equation}

\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.

For an analytical solution, we better re-write equation \eqref{eq:ThermalDifferential2} with all terms including $T$ at the left hand side: \begin{equation} R_{th}C \frac{\partial T}{\partial t} +T=R_{th}P_{in}+T_{o}. \label{eq:ThermalDifferential3} \end{equation} and next try a generic solution like $Ae^{Bt}$. Together with some boundary conditions for $T \left ( t =0 \right )$ and an input signal for $P_{in}$, we will find a solution. Assume the input signal $P_{in}$ is a step function: \begin{equation} T\left ( 0 \right )=T_{o} \\ P_{in}\left ( t \right )=\left\{\begin{matrix} 0 & t < 0 \\ P_{c} & t\geqslant 0 \end{matrix}\right. \label{eq:StepFunction} \end{equation} the solution becomes \begin{equation} T=T_{o}+P_{c}R_{th}\left ( 1-e^{-\frac{t}{R_{th}C}} \right ). \label{eq:AnalyticalRC} \end{equation}

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.

heat0.jpg Fig. 3: First step of heat model in Simulink

Next, the differential equation shows us that $\partial T/ \partial t$ is equal to $\left ( P_{in} - P_{out} \right ) / C$. This is added in figure 5.

heat2.jpg Fig. 5: Third step of heat model in Simulink

Finally, the model is completed in figure 6 by noting that $P_{out} = \left ( T-T_{0} \right )/R_{th}$.

heat3.jpg Fig. 6: Complete heat model in Simulink

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

Sensor Technology TOC

These are the chapters for the Sensor Technology course:

theory/sensor_technology/st14_differential_equation_numerical_models.txt · Last modified: 2018/10/10 00:25 by