To test accuracy of discrete approximation in time the equation of diffusion has been examined. To find P(x,y,z,t) particles distribution at tn moment in time the Taylor expansion with β = 0 (the Euler ,,forward'' scheme) has been applied. The diffusion equation with the following initial condition g(., 0) = -x(x - π)y(y - π)z(z - π) for the cubic domain and g(., 0) = |(r - R0)z(z - π)| where R0 denotes a radius of the cylindrical domain has been solved. The boundary value of the P(x,y,z,t) is set as 0. Results for both domains: cubic (with uniform elements - see Fig.3a) and cylindrical (see Fig.1), this one tuned to the designed element volume set as 0.5 with the Metropolis recipe, are shown in Figs 5c, 6a and 6b. The diffusion coefficient is set as 1.0 [units]. In the case of cylindrical domain mean values of the ratio P(x, ti)/P(x, ti + 10) calculated at each point of domain for times i = 0, 10, …, 100 have the average value 1.2699 ± 0.0049. Discrepancy between analytic and numerical solutions computed at the following moments in time equals in average -0.002 ± 0.005 for t = 0, -0.03 ± 0.04 for t = 0.5, -0.05 ± 0.06 for t = 1.0 (see Fig.5c). The analytic solution has been computed on the basis of the equation where the maximum number of used Bessel functions n is 20 (which corresponds to the radius R0 = 2.22 of the domain), the upper limit of m for their zeros knm is set as 100 for each n, and eigenvalues kz of the eigenfunctions sin(kzz) varies from 1 to 100. Note that in the case of the initial-boundary value problems the FEM solution at time t = 0 is exact (i. e. states as the initial condition). Therefore, an accuracy of an analytic approximation of the series at t = 0 is checked (see Fig.6a). Comparison between a FEM result and an analytic approximation at t = 1.0 with D = 1 [units] on meshes with h0 = 0.5 and h0 = 0.3 shows greater discrepancy than at the initial state (see Figs 6b and 6c) having the mean value -0.047 ± 0.061 and -0.029 ± 0.031, respectively. It is caused by sizes of the FEM elements h0 = 0.5 and h0 = 0.3. Results are presented in arbitrary units.
Fig. 5. The picture shows a) FEM approximation (linear) of the initial-boundary value problem in the center of the cubic domain i. e. at z = π/2 and at the time t = 0.19 with Δt = 0.01 [units] vs. the analytic result together with b) a distribution of differences between numerical and analytic solutions obtained for each node in Ω domain; c) FEM approximation (linear) of the diffusion equation for the cylindrical domain at z=π/2 and at the following times: t0 = 0 (blue), tmid = 0.5 (black) and tend = 1.0 (red) with Δt = 0.01 [units]; d) volume profile of elements within the cylindrical domain where V0 = 0.015 and h0 = 0.5; mesh quality were enhanced with help of the Metropolis algorithm.
Fig. 6. The picture shows FEM approximation (linear) of the initial-boundary value problem vs. the analytic solution depicted in the center of the cylindrical domain i. e. at z = π/2 ± 0.1 and at the times t = 0 a) and t = 1.0 b) and c). Computations have been done with Δt = 0.01 [units], the diffusion coefficient is set as D = 1.0 [units], and the element size as h0 = 0.5 b) and h0 = 0.3 c); mesh quality has been enhanced with help of the Metropolis algorithm. To better clarity of the illustration that both solutions (i.e. FEM and anal.) are getting closer each other with decaying h0, in the c) picture the number of presented nodes has been chosen smaller than in the two other cases.