The Monte Carlo method simulates neutron transport with minimal assumptions, providing high accuracy but at significant computational cost. The improved deterministic Truncation of Monte Carlo (iDTMC) method addresses this limitation by accelerating the convergence of the fission source distribution using partial-current based Coarse Mesh Finite Difference (p-CMFD) and obtaining a pin-wise reactor solution from a partial-current based Fine Mesh Finite Difference (p-FMFD) equation. The parameters used to construct the p-FMFD problem are tallied during transport, and thus the solution itself also contains uncertainty. To quantify this uncertainty, multiple parameter sets are sampled and the variance of the corresponding solutions are used as an estimate of the real variance. However, these parameters exhibit strong node-wise and parameter-wise correlations, that must be preserved for accurate estimation. To achieve this, we introduce the implicit Correlated Sampling (iCS) method using Dirichlet distribution that implicitly conserves these correlations when generating parameter sets. With this approach, the variance of iDTMC solutions can be reliably estimated without the need for independent batch calculations. This development enhances the reliability of the iDTMC method in reactor analysis.

This article proposes a dynamic model of a coupled thermoelasticity problem under stress. A coupled boundary value problem is formulated, consisting of three differen tial equations for the stress and temperature tensor components, as well as a heat flux equation with corresponding initial and boundary conditions. Explicit and implicit finite difference equations are developed, solved by successive application of the sweep method along the coordinate axes and recurrence relations, respectively. The coupled dynamic thermoelasticity problem for an anisotropic rectangle under stress is solved numerically. A comparison of the results of explicit and implicit difference schemes demonstrates the reliability of the obtained numerical results and the validity of the proposed coupled dynamic boundary value problem of thermoelasticity under stress.

Harmonic average method has been widely utilized to deal with heterogeneous coefficients in solving differential equations. One remarkable advantage of the harmonic averaging method is that no derivative of the coefficient is needed. Furthermore, the coefficient matrix of the finite difference equations is an M-matrix which guarantees the stability of the algorithm. It has been numerically observed but not theoretically proved that the method produces second order pointwise accuracy when the solution and flux are continuous even if the coefficient has finite discontinuities for which the method is inconsistent ($O$(1) in the local truncation errors). It has been believed that there are some fortunate error cancellations. The harmonic average method does not converge when the solution or the flux has finite discontinuities. In this paper, not only we rigorously prove the second order convergence of the harmonic averaging method for one-dimensional interface problem when the coefficient has a finite discontinuities and the solution and the flux are continuous, but also proposed an improved harmonic average method that is also second order accurate (in the $L^∞$ norm), which allows discontinuous solutions and fluxes along with the discontinuous coefficients. The key in the convergence proof is the construction of the Green’s function. The proof shows how the error cancellations occur in a subtle way. Numerical experiments in both 1D and 2D confirmed the theoretical proof of the improved harmonic average method.

A variety of methods that provide approximate piecewise- analytical solutions to initial-value problems governed by scalar, nonlinear, first-order, ordinary differential equations is presented. The methods are based on fixing the independent variable in the right-hand side of these equations and approximating the resulting term by either its first- or second-order Taylor series expansion. It is shown that the second-order Taylor series approximation results in Riccati equations with constant coefficients, whereas the first-order one results in first-order, linear, ordinary differential equations. Both approximations are shown to result in explicit finite difference equations that are unconditionally linearly stable, and their local truncation errors are determined. It is shown that, for three of the nonlinear, first-order, ordinary differential equations studied in this paper that are characterized by growing or decaying solutions, as well as by solutions that first grow and then decrease, a second-order Taylor series expansion of the right-hand side of the differential equation evaluated at each interval’s midpoint results in the most accurate method; however, the accuracy of this method degrades substantially for problems that exhibit either blowup in finite time or quadratic approximations characterized by a negative radicand. It is also shown that methods based on either first- or second-order Taylor series expansion of the right-hand side of the differential equation evaluated at either the left or the right points of each interval have similar accuracy, except for one of the examples that exhibits blowup in finite time. It is also shown that both the linear and the quadratic approximation methods that use the midpoint for the independent variable in each interval exhibits the same trends as and have errors comparable to the second-order trapezoidal technique.

Usually, the solution of a plane problem of the theory of elasticity in stresses is reduced to solving a biharmonic equation for the Airy stress function. In this paper, two (A and B) variants of plane boundary value problems of the theory of elasticity are formulated directly in terms of stresses. In the first case (A), the boundary value problem consists of two equilibrium equations and one Beltrami-Michell equation with the corresponding boundary and additional boundary conditions. In the formulation of the second plane boundary value problem (B), in contrast to the first, the equations of equilibrium differentiated with respect to x and y, respectively, are used. Symmetric finite-difference equations are constructed and the known Timoshenko-Goodier problem of stretching a rectangular plate with a parabolic load is solved for comparison. The discrete analogs of boundary value problems A and B are composed by the finite-difference method and the iterative method and the marching method are used to solve them. By comparing the numerical results of boundary value problems, A and B obtained by two methods, the validity of the formulated boundary value problems and the reliability of the obtained results are ensured.

In this paper, the theory and simulation results are presented for 3D vector cylindrical rotationally symmetric electromagnetic wave propagation in an isotropic nonlinear medium using a modified finite-difference time-domain general vector auxiliary differential equation method for nonlinear polarization vectors, the nonlinear dipole moment per unit volume, including both an isotropic part and an anisotropic part. The theory is presented for both the transverse magnetic and transverse electric cases, along with the combined equations. The simulation results for transverse electric spatial soliton propagation in BK7, and the simulation results for transverse electric and transverse magnetic spatial soliton propagation in fused silica, assuming a purely isotropic polarization vector, are shown. The simulation results for both transverse magnetic and transverse electric spatial soliton propagation in carbon disulfide, including the anisotropic part of the polarization vector, are also shown.

The transient process of growing axisymmetric crystals by the liquid-encapsulated Czochralski method is considered. The mathematical model accounts for heat transfer in the crucible-crystal-melt-encapsulant, formation of the melt/encapsulant meniscus, crystallization interface movement, and changes in crystal radius. A new algorithm was developed to determine the lateral surface shape of the crystal during the process. The proposed numerical approach utilizes a geometrically conservative difference scheme that guarantees the fulfillment of conservation laws of energy and mass. The special splitting technique is used to solve the corresponding set of finite difference equations. The proposed approach ensures the consistency of crystal shape evolution with conservation laws. The designed numerical procedure is used to evaluate the impact of external thermal regime on the shape of the growing crystal. When the heater temperature is maintained constant, the crystal radius gradually decreases over time. To study the influence of the external temperature field on the shape of the lateral surface, the mathematical model is supplemented with a proportional-integral temperature controller equation that links the change in heater temperature to the radius of the growing crystal. In a general case, the application of an integral temperature controller leads to fluctuations in the crystal radius around a set value, with both the frequency and amplitude of fluctuations increasing progressively. Based on the results of the numerical experiments, the parameter values for a proportional-integral temperature controller that ensure the growth of crystals with a nearly constant radius are determined. The results of transient numerical simulations are compared with the results obtained using the quasi-steady-state model of crystal growth.

This article is devoted to modeling of the gas filtration process in a dynamically interconnected multilayer porous medium. The article is devoted to modeling the process of gas filtration in a dynamically interconnected multilayer porous medium. In the article, the process of gas filtration in a heterogeneous three layer porous medium with low-permeability intermediate layers and the dynamic interaction between the layers are described by a mathematical model based on a system of differential equations of parabolic type. This mathematical model is numerically simulated using finite difference methods, i.e. explicit and implicit schemes. Since the resulting system of finite difference equations is nonlinear with respect to the pressure function, a quasilinear method was used. The dynamics of the pressure function over time was analyzed for time intervals of 360, 720 and 1080 days, and during this period the pressure distribution in the layers, the rate of pressure drop around the well and the dynamics of interlayer interaction were studied. The calculation results are presented in numerical and graphical form, which accurately reflect how the interlayer movement of the gas flow occurs. Using graphical analysis, the time step limit was determined, ensuring the stability of the computational process in the explicit scheme: a stable calculation is carried out only with a dimensionless time step Δt ≤ 1.7e-4. Also, the calculations carried out using the implicit scheme showed that this method has more stable stability compared to the explicit scheme. The results show that with a large permeability coefficient of the formation, the pressure distribution accelerates, and in the wells the pressure drop slows down. At the same time, in directly connected multilayer porous media, the permeability coefficient of the layers plays an important role. Based on the obtained results, it is possible to carry out calculations for various parameters to improve the efficiency of gas field development. It is also possible to analyze and forecast oil and gas deposits using software created on the basis of numerical models and algorithms developed in the article.

A computer program has been developed using PYTHON for the analysis of multi turn CLPHP made up of copper with capillary dimensions such as 2mm and 3.1 mm inner and outer diameters. The central difference finite difference equations are used to determine the temperatures at various sections by using numerical methods. The heat transfer equations for each cell solved iteratively, and the temperature values are updated at each time step until the solution converged. The performance of the heat pipe was analysed by varying the heat inputs and working fluid conditions. PYTHON code is developed for different fluids by using libraries such as NumPy, SciPy, and Matplotlib to obtain temperature matrix. Thermal resistance as a performance parameter is calculated with the help of temperatures. The obtained temperatures by running Python code are compared with the existing experimental results.

Abstract Slug calorimeters are commonly applied to surface heat flux measurements in hypersonic flight and ground tests. Based on measurement principle of slug calorimeters, analytical solution for one-dimensional heat conduction is applied for obtaining surface heat fluxes. However, due to the differences of thermophysical properties between substrate, slug and insulation sleeve, lateral heat conduction occurs that lead to measurement error of surface heat flux by calorimeter. Therefore, quantitative analysis of lateral heat transfer should be carried out in slug calorimeter and the interface between insulation sleeve and model wall. In this paper, the multi-dimensional heat conduction between six materials of model wall and four materials of insulation sleeve under long-term heating were investigated by numerical simulations. Then the structural designs were implemented for slug calorimeter and model wall considering the lateral heat transfer characteristics. And the derivation of axisymmetric finite difference equation was obtained by improved Crank-Nicolson discrete scheme. We found the best model wall materials for thermal matching with slug calorimeter is aluminum. And some structural optimizations can maintain a stable inverted heat flux curve.The conclusions of this study offer a practical guidance for selection of model wall material and design optimization of slug calorimeters.

The article deals with the plane stress state of elastic thin orthotropic plates in the form of irregular triangles; orthotropy is assumed to be both physical and constructive. The research method used is the numerical finite difference method using a grid of scalene triangles. For such a grid, the authors have obtained resolving finite difference equations; in difference form they have obtained correct records of boundary conditions on the plate edges through the stress function based on the frame analogy taking into account the material orthotropy; typical finite difference equations are presented that allow solving problems of the plane stress state of triangular plates with a high degree of automation. As an illustrative example, a numerical calculation of triangular plates with a grid density of N = 8 is performed using a computer; the result of the study is the analysis of the stress state in the calculated grid nodes with a wide variation in the values of the a lateral edges angles of inclination to the base of the triangle, the orthotropy coefficients. The results obtained demonstrate the feasibility and effectiveness of using a finite difference scheme on irregular triangular grids for analyzing the plane stress state of orthotropic plates. The developed approach provides a solid theoretical foundation for modeling stress distributions in structures with geometric and material anisotropy. The flexibility of the method allows adapting it to a wide range of boundary conditions and geometric configurations, laying the groundwork for further analytical development and integration into engineering software tools for structural analysis.

This paper explores the application of numerical methods to solve the two-dimensional heat equation, focusing on modelling the urban heat island (UHI) effect in Beijing, China. By utilizing the finite difference method (FDM) and solving finite difference equations iteratively using Python, our team simulated temperature variations between urban and rural zones and the heat transfer across different areas. Our results confirmed the presence of heat flow from urban to rural areas and the intensifying UHI effect over time.

Abstract This paper investigates the impact of accidental event at the gas station which is situated in Dnipro City. The processes of chemical and thermal air pollution were simulated on the basis of developed numerical models. To simulate chemical and thermal air pollution mass conversation equation and energy equation were used. For the numerical integration of governing equations finite difference schemes of splitting were used. Also, the process of fragments scattering which appears as the result of explosion at the gas station was modelled. To simulate fragments scattering Newton second Law was used. To solve the governing equation Euler’s method was used. Results of numerical experiments are presented.

The “pile-wall” structural system involves the use of conventional temporary retaining piles as part of the underground structure during its operational phase. While this approach has been implemented in engineering practice, there is a significant research gap regarding the theoretical and numerical understanding of the internal force distribution and load-sharing characteristics of the “pile-wall” system, especially in relation to the reinforcement bars used as connection nodes at the interface between the retaining piles and basement walls. This study addresses this gap by proposing a composite structural model that integrates an elastic foundation beam with a continuous beam to simulate the mechanical behavior of the “pile-wall” system. A finite difference equation and computational method are developed to analyze the internal forces and displacements in the system during two key stages: excavation and completion of the basement wall construction. The analytical solutions are validated through comparison with finite element simulations, which show a high degree of consistency between the results. The main contributions of this study include the development of a reliable calculation method for analyzing “pile-wall” internal forces and providing new insights into the load-sharing mechanisms of the system. These findings contribute to a deeper understanding of the mechanical behavior of the “pile-wall” structure and offer a solid theoretical foundation for its broader application in engineering design.

With the development of pulse technology, reinforced concrete buildings are exposed to increasingly complex high-power electromagnetic pulse (EMP) environments, posing risks of functional degradation or destruction of indoor electronic equipment and systems. Therefore, it is imperative to assess the internal fields of buildings under EMP irradiation. The challenge lies in the multi-scale characteristics of reinforced concrete buildings, where fine grids are required for the accurate modelling of rebar, thereby consuming substantial computing resources. To address this challenge, this paper proposes a fast calculation method of the time–domain coupling characteristics between buildings and EMPs based on the electromagnetic parameter equivalence of reinforced concrete walls. The method first calculates the equivalent electromagnetic parameters from the S-parameters of the walls, which are then fitted into polynomial rational functions. Then, the auxiliary differential equation finite-difference time–domain (ADE-FDTD) method is used to analyze the time–domain coupling characteristics of reinforced concrete walls and buildings under EMP irradiation. The results show that the proposed method significantly enhances computational efficiency while maintaining high accuracy. Specifically, for a large two-story reinforced concrete building, the method achieves a 3.2-fold increase in computational speed and a 4.3-fold reduction in memory usage compared to conventional commercial software (CST Studio Suite 2022). This approach provides an effective solution for simulating the coupling characteristics between large reinforced concrete buildings and external EMPs.

Consider a square free polynomial P ∈ Q [ x ] P\in \mathbb {Q}[x] of degree n n and whose coefficients have binary length at most h h . We present an algorithm computing the linear relations with coefficients in Q \mathbb {Q} between the roots of P P , in polynomial time in n , h n,h . We also present an algorithm for computing multiplicative relations between the roots of P P , also running in polynomial time in n , h n,h . Previous methods were running in polynomial time in n ! n! . Alongside, we present how to build polynomials having many relations between their roots (additive or multiplicative). Finally, we present several applications: the elementary computation of hyperelliptic integrals, the Galois group computation of differential and difference equation with constant coefficients, and the computation of torsion number for C 2 C^2 finite difference equations.

This research explores the effects of different plate thicknesses on the three-dimensional temperature distribution modeling in Butt-Joint welding using the finite difference approximation. Through computational simulations, the study investigates how varying plate thickness influences temperature profiles within welded structures. The findings offer insights into the thermal behavior of welds with varying thicknesses, aiding in the optimization of welding processes for enhanced efficiency and quality. Key aspects addressed include the formulation of governing heat transfer equations, treatment of boundary conditions, and integration of welding parameters and material properties. Validation of the model against experimental data ensures its reliability and applicability across various welding scenarios. By providing insights into the intricate temperature dynamics during welding processes, this research facilitates improved weld quality, reduced defects, and enhanced process efficiency. Introduction: Welding is a key process used in industries like automotive, aerospace, construction, and manufacturing. The quality of welded joints is essential for the reliability and performance of structures and components. Proper temperature control in the weld is vital for achieving high-quality welds. Temperature distribution is affected by factors like welding parameters, material properties, joint configuration, and heat source characteristics. Understanding these factors is critical for optimizing welding processes and ensuring defect-free welds. Numerical modeling, particularly the finite difference method (FDM), helps simulate heat transfer and predict temperature distribution accurately. Objectives: The objective of this research is to investigate the effects of plate thickness on the three-dimensional temperature distribution during Butt-Joint welding using the finite difference method. Specifically, the study aims to: Model the temperature profiles within welded structures for varying plate thicknesses. Analyze the influence of welding parameters, material properties, and joint configurations on temperature distribution. Optimize welding processes to enhance weld quality and minimize defects. Validate the developed model through comparison with experimental data to ensure its reliability and applicability across different welding scenarios. Provide insights into the thermal behavior of welds, contributing to improved process efficiency and quality control in welding operations. Methods: This study uses the heat conduction equation to model heat transfer in submerged arc welding. A fine mesh (1 mm resolution) is applied for x- and z-axes, while a coarser mesh is used along the welding axis. Only half of the plate is considered for computation. The arc efficiency (η) is 0.9, with heat distribution modeled as ηVI. Boundary Conditions: Heat source along the y-axis. Convective heat dissipation along edges and surfaces. Convection Coefficient Newtonian convection cooling is applied: h = 4.5 × 10⁻⁴ W/mm²K (within 50 mm of the weld line). h = 1.8 × 10⁻⁵ W/mm²K (remaining plate area). Material Properties: Carbon-manganese steel is used, with temperature-dependent thermal properties. Finite Difference Approximation (FDA) FDA replaces differential equations with finite difference equations (FDE): Discretize the domain into a grid. Convert the heat equation to FDE. Formulate and solve equations using Leibmann’s iteration method. A central difference implicit scheme ensures stability. The iteration scheme incorporates temperature-dependent properties. The computational code is implemented in C. Results: The model used grid systems for different plate thicknesses: 127×29×8 (6 mm), 132×36×10 (8 mm), 127×27×12 (10 mm), and 122×32×14 (12 mm). Temperature history was computed over 1500-time steps (0.2 s each) using heat flow equations and material properties. Heat input followed welding parameters. Figure shows the temperature distribution for a 6 mm plate, and presents contour plots for a 10 mm plate. Peak temperature decreases away from the weld line, with a fusion zone and heat-affected zone (HAZ). Thicker plates dissipate more heat, reducing fusion depth and increasing cooling rates. Experimental Verification: Bead-on-plate welding experiments validated the numerical model using C-Mn steel samples of varying thicknesses. Welding parameters were applied in submerged arc welding. Flux and filler metal compositions were analysed via scanning electron microscopy. Temperature was measured with K-type thermocouples on upper and lower surfaces and recorded at one-second intervals using an Agilent 34970A Data Acquisition system. Welding was performed with direct current electrode positive polarity and a 25 mm contact tube-to-workpiece distance. A total of 76 mild steel samples (6 mm, 8 mm, 10 mm, 12 mm) were tested, with peak temperature deviations under 3% at 15 mm from the weld line. Figures 8 and 9 confirm a close match between numerical and experimental results, validating the heat flow model. Conclusions: This study develops and validates a theoretical heat flow model for welding using experiments on C-Mn steel (6–12 mm thick). Key findings include: An implicit central difference approximation effectively simulates 3D thermal cycles, closely matching real conditions. The model accurately captures temperature flow, including convection effects and material property variations. Numerical results align with experiments (deviation <3%), ensuring reliable predictions for weld microstructures, distortion, and residual stress. Cooling rates increase with thickness, raising hydrogen embrittlement risks but stabilizing faster in thicker plates. The model proves effective for welding process optimization and material control.

This work presents a detailed theoretical analysis of the multiple-relaxation-time (MRT) lattice Boltzmann method (LBM), formulated on central moment (CM) space, for the numerical modeling of the one-dimensional advection-diffusion equation (ADE) with a constant velocity and diffusion coefficient, based on the D1Q3 lattice. Other LBM collision operators, such as single-relaxation-time Bhatnagar-Gross-Krook (BGK), regularized (REG) and MRT in raw moment (RM) space are also considered in this study. Without recurring to asymptotic analyses, such as the Chapman-Enskog expansion, we investigate the approximation of the MRT-CM with respect to the ADE by deriving its equivalent finite difference (EFD) scheme, which obeys an explicit four-level finite difference scheme at discrete level. Its steady-state limit follows a standard central differencing scheme for the steady ADE, yet with possible artefacts in the effective diffusion coefficient. Then, through the Taylor expansion of the EFD scheme, a detailed accuracy analysis, based on the equivalent partial differential (EPD) equation, reveals the leading order truncation errors associated with each collision model under study. Although MRT-CM and MRT-RM models have similar error structures, the former has a much reduced and simpler form, particularly in the dispersion error term, which might explain the improved Galilean invariance of the CM model. Through a suitable combination of the MRT free parameters (either in RM or CM bases), it is possible to improve its accuracy from second- to fourth-order. After that, we study the necessary and sufficient stability conditions of the MRT-CM, and its relation with other collision operators, based on the von Neumann stability analysis of the derived EFD schemes. Unexpectedly, the MRT-CM appears to support a narrower stability domain than the MRT-RM model, particularly at higher advection velocities, which can be tracked down to the inclusion of additional terms in the stability condition of the former that scale with higher order polynomials of the advection velocity. Finally, some numerical tests for the ADE on 1D unbounded domains are conducted, which confirm this work theoretical conclusions on the MRT-CM performance.

The shape control of a discretely represented curve (DRC) in the static-geometric method can be achieved not only by varying the functional external load but also through the coefficients in computational templates. These templates form the basis for constructing systems of finite-difference equations for DRC formation and indicate the proportional contribution of adjacent nodes to the desired formation. This article proposes a general approach to creating computational templates for modeling geometric objects (GOs) using superpositions of point sets. This aims to further study the influence of superposition coefficients, both arbitrary and of adjacent nodes of numerical sequences, on the formation of discrete analogs of elementary functional dependencies. One of the objectives of this study is to continue exploring the modeling of discrete geometric objects (DGOs) based on the classical finite difference method, the static-geometric method, and the geometric apparatus of superpositions. Since any polynomial of degree n is defined by n+1 points, determining the ordinate of any point given its abscissa requires substituting the coordinates of n+1 points into the polynomial function equation. This results in a system of algebraic equations containing n+1 equations and n+1 variables. Solving this system yields the polynomial coefficients a0 , a1 , a2 , a3 , … , an are found. In contrast to this approach, the recursive formula and the formula for determining the superposition coefficients proposed in this study allow for calculating the ordinate of any point of a polynomial of degree n given its abscissa without constructing and solving a system of n+1 equations. The ordinate of any curve point is determined as a superposition of the ordinates of n+1 points. In the proposed method of geometric curve modeling, the superposition coefficients are derived from systems of equations that contain one equation fewer than those used to calculate polynomial coefficients . Computational templates have been developed for the discrete formation of polynomial functional dependencies using superpositions of adjacent points’ coordinates. The approach presented in the article can be used to obtain expressions similar to formula (3) for calculating superposition coefficients for adjacent points of polynomials with two variables. Varying the superposition coefficients in the developed computational templates allows for studying the impact of these coefficients, both arbitrary and for adjacent nodes of numerical sequences, on the formation of discrete analogs of elementary functional dependencies. The shape control of a discretely represented curve (DRC) in the static-geometric method can be achieved not only by varying the functional external load but also through the coefficients in computational templates. These templates form the basis for constructing systems of finite-difference equations for DRC formation and indicate the proportional contribution of adjacent nodes to the desired formation. This article proposes a general approach to creating computational templates for modeling geometric objects (GOs) using superpositions of point sets. This aims to further study the influence of superposition coefficients, both arbitrary and of adjacent nodes of numerical sequences, on the formation of discrete analogs of elementary functional dependencies. One of the objectives of this study is to continue exploring the modeling of discrete geometric objects (DGOs) based on the classical finite difference method, the static-geometric method, and the geometric apparatus of superpositions. Since any polynomial of degree n is defined by n+1 points, determining the ordinate of any point given its abscissa requires substituting the coordinates of n+1 points into the polynomial function equation. This results in a system of algebraic equations containing n+1 equations and n+1 variables. Solving this system yields the polynomial coefficients a0 , a1 , a2 , a3 , … , an are found. In contrast to this approach, the recursive formula and the formula for determining the superposition coefficients proposed in this study allow for calculating the ordinate of any point of a polynomial of degree n given its abscissa without constructing and solving a system of n+1 equations. The ordinate of any curve point is determined as a superposition of the ordinates of n+1 points. In the proposed method of geometric curve modeling, the superposition coefficients are derived from systems of equations that contain one equation fewer than those used to calculate polynomial coefficients . Computational templates have been developed for the discrete formation of polynomial functional dependencies using superpositions of adjacent points’ coordinates. The approach presented in the article can be used to obtain expressions similar to formula (3) for calculating superposition coefficients for adjacent points of polynomials with two variables. Varying the superposition coefficients in the developed computational templates allows for studying the impact of these coefficients, both arbitrary and for adjacent nodes of numerical sequences, on the formation of discrete analogs of elementary functional dependencies.

In this paper the finite-difference time-domain general vector auxiliary differential equation method [Opt. Express14, 8305 (2006)10.1364/OE.14.008305], nonlinear polarization vector, the nonlinear electric dipole moment per unit volume, is extended to include anisotropy, in nonlinear isotropic media at optical frequencies. The theory is presented for extending the numerical method in 3D Cartesian coordinates, and then example simulation results are presented for two isotropic media. First, the simplified 2D transverse magnetic case is revisited for the fused silica example introduced in the 2006 paper, including the anisotropic part of the nonlinear polarization vector in the simulation; the simulation results including the anisotropic part of the polarization vector were compared with the purely isotropic polarization vector simulation results. Second, a simplified 2D transverse magnetic example was simulated in carbon disulfide, with its strong molecular re-orientation-induced polarization anisotropy.