3.7 Summary
4.2.1 The Heat equation
The one dimensional heat equation (4.8) describes the evolution of the temperature distribution of a finite rod, where the cross section of the rod is sufficiently small, the rod is totally insulated and both ends are held at zero temperature.
du
Wherel is the length of the rod, kis the heat conductivity, u(x, t) is the temperature distribution of the rod at time t while f(x) is the initial temperature distribution.
Equation (4.8) can be solved analytically and the solution is the following [21]:
u(x, t) =
This infinite series can be used only when the initial conditions are simple and the integral in the computation of the an values can be determined. For simplicity the length of the rod and the heat conductivity is set to 1 and the initial temperature distribution is a reversed parabolic curve described by the following equation:
f(x) = 4x(1−x) (4.10)
After substituting (4.10) into (4.9) and performing the integration the following equa-tion is obtained:
an =−16 cos (nπ)
n3π3 −8 sin (nπ)
n2π2 + 16
n3π3 (4.11)
Figure 4.8: Analytical solution of the heat equation
where sin(nπ) is zero if n is integer and cos(nπ) is -1 ifn is odd while it is +1 if n is even. Thus an is zero for every even value ofn and (4.11) can be simplified:
a2n = 0 a2n+1 = 32
n3π3
(4.12)
After computing the an values (4.8) can be used to approximate the analytical solu-tion. Fortunately only the first 30 elements should be computed for every point of the rod to achieve 10−14accuracy. The analytical solution is computed on 1025 points of the rod in every 2−8s the results of the first 0.5s are shown in Figure 4.8.
Equation (4.8) must be spatially discretized to solve it on CNN architecture.
We used the following well-known approximation of the Laplace operator and the corresponding CNN template:
∂2u
∂x2 = k
∆x2 (ui+1−2ui+ui−1) + O ∆x2
(4.13a) A= k
∆x2
1 −2 1
(4.13b) where ui is the state of theith finite element and ∆x is the spatial distance between the finite elements. The discretized version of (4.8) is solved by the numerical methods
Euler
Figure 4.9: Error of the numerical methods using different spatial discretizations examined in the previous section. Both methods were run by using different spatial discretizations and timesteps and the results were compared to the analytical solution.
The error values are shown in Figure 4.9.
The results show that in this case the main source of errors is the spatial dis-cretization. If the number of elements is small, the accuracy of the solution cannot be increased by using smaller timesteps. Moreover the error cannot be decreased by using the higher order 2nd and 4th order Runge-Kutta methods. The simple forward Euler method seems to be suitable to solve the heat equation thus only this method will be examined by using different fixed-point precision. The absolute maximum difference between the analytical and the fixed-point computations are shown in Figure 4.10.
The characteristics of the error functions are very similar to the errors of the solu-tion of the mechanical vibrating system. If the fixed-point precision is high enough, its error is equal to the error of the floating-point computation and an optimal bit width can be found for every grid size. If the number of grid points is increased without increasing the computing precision, the accuracy of the solution decreases
(a)
Figure 4.10: Error of the forward Euler method
due to the rounding errors. The distribution of the errors in space and time is shown in Figure 4.11.
The error value is increasing very quickly to the maximum error value at the beginning of the simulation. Later the error is slowly decreasing and increasing again at the end of the simulation. This increase is the result of the inaccuracy of the fixed-point number system when the state values are very small. It is also possible that the solution is ’frozen’ when the computed derivatives are smaller than the smallest representable number.
In our case simple initial conditions are selected and the an coefficients in (4.8) can be computed easily. Unfortunately in most cases computation of the analytical solution is very hard or impossible. Therefore the fixed-point results must be com-pared to the floating-point results using the same grid size and timestep value. The result of this comparison with the 64 and 32 bit floating-point results is shown in Figure 4.12. and Figure 4.13. along with the difference between the analytical and the floating-point solution.
The error functions do not follow the error of the 64 bit floating-point computation but cross it similarly to the error functions of the mechanical vibrating system which was obtained in the previous section. In the 32 bit floating-point case the difference cannot be decreased by increasing the state width over 32 bit. Therefore the heuristic methods described in the previous section can be used to determine the optimal bit width for the fixed-point computations. If the accuracy of the solution can be estimated or it is given a-priori, this value can be used in the qualification of the fixed-point solution. If the accuracy of the fixed-point and floating-point solution should be equal, the break-point of the error function should be determined.
Figure 4.11: Distribution of the error of the 44 bit fixed-point solution on a 1025 element grid
(a)
1.0E-15 1.0E-14 1.0E-13 1.0E-12 1.0E-11 1.0E-10 1.0E-09 1.0E-08 1.0E-07 1.0E-06 1.0E-05 1.0E-04 1.0E-03 1.0E-02
17 33 65 129 257 513 1025
Number of elements
Error
FP 64 8 bit 12 bit 16 bit 20 bit 24 bit 28 bit 32 bit 36 bit 40 bit 44 bit 48 bit 52 bit 56 bit
(b)
1.0E-15 1.0E-14 1.0E-13 1.0E-12 1.0E-11 1.0E-10 1.0E-09 1.0E-08 1.0E-07 1.0E-06 1.0E-05 1.0E-04 1.0E-03 1.0E-02
8 12 16 20 24 28 32 36 40 44 48 52 56
State width (bit)
Error
17 33 65 129 257 513 1025
Figure 4.12: Difference between the 64 bit floating-point and the fixed-point solutions
(a)
Figure 4.13: Difference between the 32 bit floating-point and the fixed-point solutions