ABSTRACT

We investigate under what conditions transient simulation could be used to integrate backward in time so that the initial field could be recovered from later histories. In this paper we use realistic examples and find that, in long histories, traces of the initial field would be present only in the exact analytical solutions. We conclude that the recovery of initial field is possible only if the equations could be solved analytically or only short time periods are involved. In practice, it is not possible to detect those traces by measurements or observations. If numerical procedures are used, truncation and discretization errors are always present. Fine-tuning of system parameters used or transforming time into another pseudo time frame may allow numerical integration to be carried out backward in time. But numerical instability is still a problem. Large spurious increases found by numerical procedures are most likely due to numerical inaccuracy and instability.

Keywords:

Solution of Differential Equations, Transient Simulation, Numerical Instability, Mathematical Modelling, Backward Time Integration

1. Introduction

With the availability of modern computer power, there is a rapidly increased use of mathematical models to simulate real life or physical situations [1] [2] . Mathematical models use mathematical concepts and language to represent real problems and are used in natural sciences and engineering disciplines, as well as in the social sciences. Often, solutions of the mathematical equations used in the models may not be found easily. A proven technique is to employ a pseudo transient simulation [3] [4] . The word pseudo is used to indicate that the solution is found by a transient approach but the time used is not necessarily real. However, as far as the mathematical principles and methods are concerned there is no distinction between real and pseudo time.

With success of mathematical modelling, there is a tendency to replace laboratory testing of designs by numerical experiments [5] and operator control of a system by computer control [6] . It should be noted that the accuracy and effectiveness of computer-based approach are most dependent on the model or models used. For example, maximum transient stress found by the numerical solutions of a simulation based on the long cylinder approximation would be five folds smaller than the maximum stress which actually exists at the ends of the cylinder [7] . The interpretation of results obtained by simulations is also important. Often the averages over a computational grid cell are used for numerical reasons. This means that the actual maximum inside a cell could be much greater than the cell average found. This could happen such as in the case of weather forecasting; the computer program predicted wind speeds could only be halves of the actual maxima.

For saving computational time, simpler models are often used in simulation [8] . System parameters used in those models may not be the physical properties of the system but are those found to be most suitable for matching with observations in so called “fine tuning”. In such cases, it is important that the simulation is not being used far beyond the conditions employed in the calibration. In a sense this type of simulation is similar to curve fitting except that far more complicated mathematical operators are involved.

By far the more uncertain attempt involves predictions not just for what can be seen now or what would happen in the future but for predictions on what have already occurred in the past. A notable case is that astrophysicists are predicting what the universe could be billions of years ago by what they can observe now. As far as solving the field equations is concern, it is a well-established mathematical principle that solutions are possible only for the time space from zero to infinity. However, with pseudo transient treatment, the real time space from zero to negative infinity could be converted to the pseudo time space varying from zero to infinity. This would be an useful solution if it is found permissible.

In this paper, heat conduction problems are solved to give a set of analytical solutions. They are used to show why integration forward in time is always stable. Backward integration in time is permissible only if it is done analytically as any error in numerical approximates would grow exponentially. We also show that numerical backward extrapolation cannot be used to recover the initial condition even when a pseudo time is used. In another example we show that if the temperature field in details is not needed, heat conduction need not be considered. The system could be modelled by an ordinary differential equation and solved analytically. System parameters that are generally not the same as the physical properties could be found by calibration so that the system could be used for prediction, measurement or control purposes.

2. Transient Heat Conduction Equation

For a one dimensional heat conduction problem

$\frac{\partial T\left(x,t\right)}{\partial t}=\frac{{\partial }^{2}T\left(x,t\right)}{\partial x^2}+Q\left(x\right)$ (1)

where T is the temperature and Q heat source, and all are in dimensionless. The initial condition, $T\left(x,0\right)=T^0\left(x\right)$ is given and the boundary conditions are:

$T^{\prime }\left(0,t\right)=0;\frac{\partial T\left(1,t\right)}{\partial x}-hT\left(1,t\right)=0$ (2)

where h is the heat transfer coefficient.

Using N terms of chopped sine Fourier expansions

$T\left(t,x\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{a}_i{f}_i\left(t\right)\sin \left({\lambda }_{i}x+\frac{\pi }{2}\right)}$ (3)

$Q\left(x\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{q}_i\sin \left({\lambda }_{i}x+\frac{\pi }{2}\right)}$ (4)

$T^0\left(x\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{t}_i^0\sin \left({\lambda }_{i}x+\frac{\pi }{2}\right)}$ (5)

where λ_i are given by the characteristic equation derived from the boundary condition (2)

$-{\lambda }_i\cos \left({\lambda }_i+\frac{\pi }{2}\right)+h\sin \left({\lambda }_i+\frac{\pi }{2}\right)=0$ (6)

Substitute (3) and (4) into (1) and equate the coefficients to give

$a_i{f^{\prime }}_i\left(t\right)=-{\lambda }_i^2{a}_{i}f\left(t\right)+q_i$ (7)

The solution of the above equation that satisfies the initial condition (5) is

$a_i{f}_i\left(t\right)=\frac{q_i}{{\lambda }_i^2}\left[1-\exp \left(-{\lambda }^{2}t\right)\right]+t_i^0\exp \left(-{\lambda }^{2}t\right)$ (8)

$\therefore T\left(x,t\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[\frac{q_i}{{\lambda }_i^2}\left\{1-\exp \left(-{\lambda }^{2}t\right)\right\}+t_i^0\exp \left(-{\lambda }^{2}t\right)\right]}\sin \left({\lambda }_{i}x+\frac{\pi }{2}\right)$ (9)

and at $t\gg 0$, the time independent and stationary temperature is

$T\left(x,t\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{q_i}{{\lambda }_i^2}\sin }\left({\lambda }_{i}x+\frac{\pi }{2}\right)$ (10)

Obviously the initial condition, including any errors, as well as those differences between the current and the stationary temperature would decay exponentially in time leaving only the stationary components.

With the analytical solution, Equation (9), $T\left(0,x\right)$ can be recovered exactly from $T\left(t_N,x\right)$ even when $t_N\gg 0$. Using $T\left(t_N,x\right)$ as the initial value, $T\left(0,x\right)$ can be found by integrating backward to $t=0$. This is possible because the analytical solution contains exactly every small terms associating with the factor $\exp \left(-{\lambda }_i^2{t}_N\right)$ that combines with the very large factor $\exp \left({\lambda }_i^2{t}_N\right)$ to give the exact starting value.

Use as an example $Q\left(x\right)=\sin \left(x+\frac{\pi }{2}\right)$ and $T\left(x,0\right)=t^0\sin \left(x+\frac{\pi }{2}\right)$ to show what is needed in the analytical solution before the initial value could be recovered from that at time t_N. According to Equation (9) the solution is

$T\left(x,t_N\right)=\left[\left\{1-\exp \left(-t_N\right)\right\}+t^0\exp \left(-t_N\right)\right]\sin \left(x+\frac{\pi }{2}\right)$ (11)

As $t_N\gg 0$, the term $\exp \left(-t_N\right)$ is very small; it drops out from Equation(11) to give the steady solution

$T\left(x,t_N\right)=\sin \left(x+\frac{\pi }{2}\right)$ (12)

If Equation (12) is used as the initial value, integration backward to $t=-t_N$ will give the same temperature

$T\left(x,0\right)=\left[\left\{1-\exp \left(t_N\right)\right\}+\exp \left(t_N\right)\right]\sin \left(x+\frac{\pi }{2}\right)=\sin \left(x+\frac{\pi }{2}\right)$ (13)

This should not be a surprise as starting from the steady solution there should be no change in the temperature field whether it is marching forward or backward in time. Since t⁰ is absent in the steady solution, the arbitrary initial value could not be recovered from backward integration. However, no matter how small the term $\exp \left(-t_N\right)$ is if the full Equation (11) is used as the initial value, integration backward from time t_N to time zero does recover t⁰ exactly.

Next, consider that the starting field is slightly different to the steady field. That is

$T^0\left(x,0\right)=\left(1+\epsilon \right)\sin \left(x+\frac{\pi }{2}\right)$ (14)

where $\epsilon \approx 0$. Then according to Equation (11), the backward integration from the beginning to –t_N is

$\begin{array}{c}T\left(x,t_N\right)=\left[\left\{1-\exp \left(-t_N\right)\right\}+\left(1+\epsilon \right)\exp \left(-t_N\right)\right]\sin \left(x+\frac{\pi }{2}\right)\\ =\left[1+\epsilon \exp \left(t_N\right)\right]\sin \left(x+\frac{\pi }{2}\right)\end{array}$ (15)

It could be seen from Equation (15) that any error even as small as truncation error would growth exponentially without bound as the integration is progressing further and further away in the backward time direction. This spurious solution would sooner or later swamp the real solution. It should be noted that when starting with an analytical solution, $\epsilon =0$ and hence stability is not a problem.

3. Pseudo Transient Approach

In a real and physical situation, it is impossible to run a system backward in time. But in mathematics, integration can be carried out forward as well as backward in time. For an initial value problem, however, it is generally accepted that solutions can be found only for the time space from zero to infinity. In order to comply with this condition, a pseudo time, $\theta =t_{end}-t$ is used, where t_end is a nominal positive value in the real time space. Then the pseudo time space used is from zero forward to t_end. It also corresponds to the real time space from t_end to zero. That is backward in real time space. With this pseudo time, Equation (1) becomes

$\frac{\partial T\left(x,\theta \right)}{\partial \theta }=-\frac{{\partial }^{2}T\left(x,\theta \right)}{\partial x^2}-Q\left(x\right)$ (14)

The solution based on chopped sine Fourier expansion in Equations (3)-(5) is

$T\left(\theta ,x\right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[\frac{q_i}{{\lambda }_i^2}\left\{1-\exp \left({\lambda }^2\theta \right)\right\}+t_i^0\exp \left({\lambda }^2\theta \right)\right]}\sin \left({\lambda }_{i}x\right)$ (15)

It is apparent that numerical integration forward in pseudo time would have a numerical stability problem as any error present at the beginning would growth exponentially and swamp the solution itself. However, integration backward in time would be stable, confirming the condition that for this type of field equation numerical integration in real time should only be carried out in the forward direction.

To investigate the situation when integration is done by finite difference, suppose that Equation (14) is to be solved with $Q\left(x\right)=q\sin \left(2x+\frac{\pi }{2}\right)$ and $T\left(x,0\right)=t^0\sin \left(2x+\frac{\pi }{2}\right)$. Then, only one term in the Fourier series needs to be considered. Let the left hand side be replaced by Crank-Nicolson finite difference.

At $t=t^i$, $T\left(x,t=t^i\right)=T^i\sin \left(2x+\frac{\pi }{2}\right)$ with superscript i denotes the time step. Then $\frac{T^{i+1}-T^i}{\Delta \theta }=2\left(T^{i+1}+T^i\right)-q$

$\therefore T^{i+1}=\frac{1+2\Delta \theta }{1-2\Delta \theta }{T}^i-q\Delta \theta \simeq \left(1+4\Delta \theta \right)T^i-q\Delta \theta$ (16)

If an inaccurate ${\tilde{T}}^i=T^i+\epsilon$ is used in Equation (16), the estimated

${\tilde{T}}^{i+1}=\left(1+4\Delta \theta \right)\left(T^i+\epsilon \right)-q\Delta \theta$ and

${\tilde{T}}^{i+1}-T^{i+1}=\left(1+\Delta \theta \right)\epsilon$ (17)

It can be seen that for a numerically unstable system, any error will grow in the forward time integration even when Crank-Nicolson difference scheme has been used.

4. A Point Model Approach

Consider as an example temperature in a heat generating body. If heat conduction is not considered, the mathematical equation representing such a model may not involve any boundary condition. Then it is not an initial boundary value problem, and there could be no restriction for forward or backward integration in time. But, there are still conditions that need to be fulfilled before a solution could be used to predict future or past events.

When electric current is passing through a heating element, due to the resistance heat is generated in the element. At the beginning, some of this heat is loss to the surrounding and the remainder causes a raise in the temperature. The heat loss can be calculated by knowing the heat transfer coefficient and the temperature difference between the heating element and the surrounding. The temperature rise can be calculated if the heat capacity of the element is known. A stationary temperature will be reached when the rate of heat generation is equal to that of heat loss. Using a point model approach, temperature is represented by a single value and could be found by the following ordinary differential equation:

$\frac{\text{d}T}{\text{d}t}=\frac{1}{c}\left[Q-hT\right]=c_1-c_{2}T$ (18)

where T is the temperature relative to the surrounding temperature, Q the heat generation rate, h the heat transfer coefficient, c the heat capacity and t is time. In this equation, $c_1=\frac{Q}{c}$ and $c_2=\frac{h}{c}$ are parameters characteristic of the system. Equation (18) can be cast in a normalized form where all the quantities involved are dimensionless. Although c₁ and c₂, are related to the physical properties such as conductivity, density, heat capacity as well as to the system dimensions and how cooling is being carried out, they could be directly determined from the system by calibration. The analytical solution for Equation (18) is:

$T=\frac{c_1}{c_2}+\left[T_0-\frac{c_1}{c_2}\right]\exp \left(-c_{2}t\right)$ (19)

where T₀ is the initial temperature at t = 0. The solutions for different starting temperature are shown in Figure 1.

As there is a one to one relationship between T and t, this arrangement could be used as a time-measuring device by using two temperature readings to find out the time between them. One of the essential conditions is that the temperature at the starting must be known so that the appropriate curve may be used. For example, depending on which curve is used a temperature reading of 85 will give different times since the system is turned on. Only knowing the initial temperature, the curve could be selected correctly.

For good measurement sensibility, the device must be designed to operate over the range $t<5$. When $t>10$ there is very little change of temperature with time and the device is no longer sensitive enough for measurement. On the other hand, as a controlling device, the chosen operating parameters should be such that $t>10$ could be reached within the design specification for the change to attain equilibrium quickly. Also needed is that the operation conditions must be close to those used in the calibration.

Perhaps one feature unique to this device as compared with other clocks is the fact that time lapse of a past event could be found. That is to measure time backward. This principle is used in radiometric dating that is described by the same Equation (19) with c₁ = 0. The accuracy relies on the fact that the half-life of a radioactive isotope is a constant and could be easily determined. However, to measure the time accurately the initial conditions must be known. That is only the time lapse between two known events can be found.

5. Linearization of a Non-Linear Problem

There are not many cases where non-linear problems have analytical solutions. A numerical approach would be needed to replace the problem by a linearized approximation at each forward time step. In non-linear form, Equation (18) becomes

$\frac{\text{d}T\left(t\right)}{\text{d}t}=\frac{1}{c\left(T,t\right)}\left[Q\left(T,t\right)-h\left(T,t\right)T\left(t\right)\right]=c_1\left(T,t\right)-c_2\left(T,t\right)T\left(t\right)$ (20)

The temperature and time dependency is replaced by the averages at the i^th and (i + 1)^th time step so that in linearized form

$\frac{\text{d}T\left(t\right)}{\text{d}t}={\tilde{c}}_1-{\tilde{c}}_{2}T\left(\bar{t}\right)$ (21)

where ${\tilde{c}}_1=\frac{1}{2}\left[c_1\left(T^{i+1},t^{i+1}\right)+c_1\left(T^i,t^i\right)\right]$, ${\tilde{c}}_2=\frac{1}{2}\left[c_2\left(T^{i+1},t^{i+1}\right)+c_2\left(T^i,t^i\right)\right]$, and $\bar{t}=\frac{1}{2}\left(t^i+t^{i+1}\right)$. As $T^{i+1}$ is yet to be determined, iteration is needed in the forward integration procedure. If a stationary solution does exist, this scheme should be stable. Solution of Equation (21) is in the same form as Equation (19) and the integration could be carried from tⁱ to tⁱ⁺¹ smoothly. Should there be any error such as that due to numerical truncation; the error would decay exponentially and eventually disappear from the final solution. However, if a backward integration in time is carried out, an error would grow exponentially, meaning that the backward time integration would be unstable also for a non-linear problem. Therefore, $T\left(t\gg 0\right)$ could not be used to recover $T\left(t=0\right)$ by backward integration due to two reasons: 1) due to exponential decay the contribution of the initial value has already become insignificant in $T\left(t\gg 0\right)$, and 2) the initial value so obtained could have come from very small error in the starting input.

The same approach could be applied to more complex field, for example, dynamic study of a rotor in a turbine based on a point model [9] , propagation of complex electromagnetic waves in optical fibres [10] and matter-wave laser [11] . In every case, providing the system is operating within parametric stability domain, stationary solutions could be found numerically from arbitrarily chosen initial guesses. In all those cases, iterations are needed to solve the linearized equations. It is not possible to back track from the stationary solutions and recover the initial values because they are arbitrarily chosen.

6. Conclusions

Every physical problem that could be described by a well-posed transient field equation starts with an initial value. For a real and conservative system this initial value would decay and eventually disappear leaving only the stationary solution. It would be normally impossible to recover the initial value from the stationary solution by working backward in time. If it was possible to identify the very small remnant of the initial value after a long time history, the recovery relies on tracing back exactly its histories at every step. This is certainly not possible in a numerical procedure.

Forward numerical integration in a conservative transient field equation is normally stable. But backward time integration of the same system is numerically unstable.

System parameters could be fine-tuned to allow backward time integration. However, problems about the recovery of the initial value and numerical stability also apply to this fine-tuned system.

No matter how well-tuned a system could be the time interval between two events cannot be accurately measured unless the fields for the two events are known.

Observations in a conservative system in equilibrium should not change in time either forward or backward. If a numerical model predicts the field intensity, for example temperature is very much larger in the past, such a prediction can only come from numerical instability and inaccuracy in the measurements.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

Cite this paper

Chen, P. (2019) On Transient Simulation of Field Equations. Applied Mathematics, 10, 719-727. https://doi.org/10.4236/am.2019.109051

References