Time Step Size Research Articles

Challenges in Numerical Solutions of Higher-Dimensional Differential Equations

Differential equations constitute a fundamental tool in modeling various natural phenomena across scientific disciplines such as physics, engineering, and finance. We provide an overview of fractional differential equations, focusing on the computational requirements associated with their numerical solutions from a computer science perspective. We analyze the computational intricacies concerning First-Order Linear ODE, First-Order Nonlinear ODE, Second-Order Linear ODE, Second-Order Nonlinear ODE, Heat Equation (PDE), and Wave Equation (PDE). This comparative assessment delves into the computational demands of solving these equations using differential equation methodologies. While analytical solutions provide deep insights, obtaining numerical solutions, particularly in higher dimensions, remains a persistent challenge. Finite difference methods commonly employed for numerical solutions, In higher-dimensional problems, traditional numerical methods face challenges stemming from an exponential surge in grid points and the consequent demand for substantially decreased time step sizes. This paper explores the challenges posed by higher-dimensional differential equations in numerical solutions. It highlights the infeasibility of finite difference methods in such scenarios and emphasizes the need for innovative numerical techniques capable of efficiently handling the complexities of higher-dimensional differential equations. Overcoming these challenges is crucial for advancing our understanding and modeling capabilities in complex real-world systems governed by differential equations. Continued research efforts strive to develop novel numerical methodologies capable of addressing these challenges, aiming to broaden the scope of solvable higher-dimensional differential equations and expand their application across diverse scientific domains.

Hybrid methods for solving structural geometric nonlinear dynamic equations: Implementation of fifth-order iterative procedures within composite time integration methods

This paper studies a class of hybrid methods that implement multi-point iterative procedures as the nonlinear solver within optimized composite implicit methods. These multi-point methods are employed to enhance convergence order and accuracy without requiring higher-order derivatives for solving structural geometric nonlinear dynamic equations. The promising approach provides an opportunity to delve deeper into its implications and explore the potential applications of multi-point techniques and composite methods across the analysis of dynamical systems. This study involves the utilization of two-point and five-point iterative methods with fifth-order convergence within three-step and four-step optimized composite implicit time integration methods. Moreover, a comprehensive hybrid scheme is introduced by integrating a family of optimized composite implicit methods with single-point and multi-point iterative methods. Then, the performance of the hybrid scheme is analyzed through structural examples, considering factors such as time step size, tolerance threshold, and the degrees of freedom of structures. The findings and numerical results illustrate that multi-point methods outperform the single-point method in terms of the number of iterations. Simultaneously, the four-step time integration method exhibits superior performance compared to the three-step time integration method, albeit with a marginal difference. It is observed that the hybrid method, consisting of the four-step time integration method and the five-point iterative method, performs better than other methods in terms of the number of iterations and errors for solving geometric nonlinear dynamic equations, especially in large structures. Thus, it stands out as a favorable choice for practical analyses and a possible candidate for generalization to other types of nonlinear problems in computational mechanics.

An efficient two-grid high-order compact difference scheme with variable-step BDF2 method for the semilinear parabolic equation

Due to the lack of corresponding analysis on appropriate mapping operator between two grids, high-order two-grid difference algorithms are rarely studied. In this paper, we firstly discuss the boundedness of a local bi-cubic Lagrange interpolation operator. And then, taking the semilinear parabolic equation as an example, we first construct a variable-step high-order nonlinear difference algorithm using compact difference technique in space and the second-order backward differentiation formula with variable temporal stepsize in time. With the help of discrete orthogonal convolution kernels, temporal-spatial error splitting idea and a cut-off numerical technique, the unique solvability, maximum-norm stability and corresponding error estimate of the high-order nonlinear difference scheme are established under assumption that the temporal stepsize ratio satisfies rk := τk/τk−1 < 4.8645. Then, an efficient two-grid high-order difference algorithm is developed by combining a small-scale variable-step high-order nonlinear difference algorithm on the coarse grid and a large-scale variable-step high-order linearized difference algorithm on the fine grid, in which the constructed piecewise bi-cubic Lagrange interpolation mapping operator is adopted to project the coarse-grid solution to the fine grid. Under the same temporal stepsize ratio restriction rk < 4.8645 on the variable temporal stepsize, unconditional and optimal fourth-order in space and second-order in time maximum-norm error estimates of the two-grid difference scheme is established. Finally, several numerical experiments are carried out to demonstrate the effectiveness and efficiency of the proposed scheme.

Analysis and estimation of time step sizes and stiffness parameters for efficient simulations of macromolecules with hydrodynamic interactions

Recent studies have shown the importance of using highly resolved models for Brownian Dynamics (BD) simulations of long macromolecules. For computational efficiency, such models use stiff springs to mimic a single Kuhn step and use a single-step semi-implicit (SS) scheme. Somewhat unexpectedly, time step sizes for such a single-step method need to be reduced with increasing chain size and level of hydrodynamic interactions (HIs), for good convergence. The conventional predictor–corrector (PC) method works reasonably well but remains computationally slow, owing to multiple iterations per time step to convergence. In this study, we reveal how the time step size for the much faster SS method is tied to the physics of the problem. Using simple physical principles, we derive an analytical estimate of the upper limit on the time step size for given levels of HI, chain size, and stiffness of connecting springs. Detailed BD simulations at equilibrium and in flow fields highlight the success of our analytical estimate. We also provide a lower limit on spring stiffness parameter such that it remains effectively rigid and successfully mimics a Kuhn step. Our investigations show that the resulting BD simulations using our estimated time step size in the SS scheme are significantly faster than the conventional PC technique. The analysis presented here is expected to be useful in general for any type of simulations of macromolecules, with or without flow fields, owing to deep connections with the underlying physics.

The generalized Kelvin chain-based model for an orthotropic viscoelastic material

We propose a constitutive material model to describe the rheological (viscoelastic) mechanical response of timber. The viscoelastic model is based on the generalized Kelvin chain applied to the orthotropic material and is compared to the simple approach given by standards. The contribution of this study consists of the algorithmization of the viscoelastic material model of the material applied to the orthotropic constitutive law and implementation into the FEM solver. In the next step, the fitting of the input parameters of the Kelvin chain is described, and at least a material model benchmark and comparison to the approach given by standards were done. The standardized approach is based on the reduction of the material rigidity at the end of the loading period using a creep coefficient, whereas the loading history state variables are not considered when establishing the result for a specific time step. The paper presents the benefits of the rheological model. It also demonstrates the fitting algorithm based on particle swarm optimization and the least squares method, which are essential for the use of the generalized Kelvin chain model. The material model based on the orthotropic generalized Kelvin chain was implemented into the FEM solver for the shell elements. This material model was validated on the presented benchmark tasks, and the influence of the time step size on the accuracy of model results was analyzed.

Energy diminishing implicit-explicit Runge–Kutta methods for gradient flows

This study focuses on the development and analysis of a group of high-order implicit-explicit (IMEX) Runge–Kutta (RK) methods that are suitable for discretizing gradient flows with nonlinearity that is Lipschitz continuous. We demonstrate that these IMEX-RK methods can preserve the original energy dissipation property without any restrictions on the time-step size, thanks to a stabilization technique. The stabilization constants are solely dependent on the minimal eigenvalues that result from the Butcher tables of the IMEX-RKs. Furthermore, we establish a simple framework that can determine whether an IMEX-RK scheme is capable of preserving the original energy dissipation property or not. We also present a heuristic convergence analysis based on the truncation errors. This is the first research to prove that a linear high-order single-step scheme can ensure the original energy stability unconditionally for general gradient flows. Additionally, we provide several high-order IMEX-RK schemes that satisfy the established framework. Notably, we discovered a new four-stage third-order IMEX-RK scheme that reduces energy. Finally, we provide numerical examples to demonstrate the stability and accuracy properties of the proposed methods.

The stabilized exponential-SAV approach for the Allen–Cahn equation with a general mobility

In this paper, we construct a second-order accurate, energy stable and maximum bound principle-preserving scheme for the Allen–Cahn equation with a general mobility based on the stabilized exponential scalar auxiliary variable (SESAV) approach. Some extra stabilizing terms are added to the discretized scheme for the purpose of improving numerical stability. We first proved the maximum bound principle (MBP) under reasonable constraints on time step size and the stabilization parameter. Then, we found that the proposed scheme is unconditionally energy-stable. Finally, a numerical example is carried out to verify the theoretical results.

A novel approach to solve hyperbolic Buckley-Leverett equation by using a transformer based physics informed neural network

Solving hyperbolic partial differential equation is a challenging task due to the non-linear feature that requires to capture the shock wave. Numerical solution relies on discretization of both spatial and temporal domain, and iterative approach like Newton's method is involved and time step size is crucial for stability and convergence in the presence of non-linearity. Physics-informed neural networks (PINNs) offer a new and versatile approach for solving partial different equations by minimizing the residual of governing equations and approaching to the initial and boundary conditions. Currently, most PINNs are built based on a simple fully connected neural network which exhibits some limitations to model complex non-linear partial differential equations. In this paper, a novel method is developed to combine Transformer model and PINNs approach (Tr-PINN) to solving a hyperbolic partial differential equation directly without any prior knowledge. Tr-PINN method is based on a series of Transformer blocks where self-attention mechanism is used to capture the non-linearity features of the solution. Unlike most PINNs models generate inputs with spatial and temporal vectors only, Tr-PINN introduces the non-linearity term mobility ratio as additional input vector. The method is tested on a classical hyperbolic problem, called Buckley-Leverett equation with non-convex flux function. We found that the Tr-PINN method can capture the water shock front effectively and provide a general solution for Buckley-Leverett equation under various mobility ratio conditions.

Numerical coupling of aerosol emissions, dry removal, and turbulent mixing in the E3SM Atmosphere Model version 1 (EAMv1) – Part 2: A semi-discrete error analysis framework for assessing coupling schemes

Abstract. Part 1 (Wan et al., 2024) of this study discusses the motivation and empirical evaluation of a revision to the aerosol-related numerical process coupling in the atmosphere component of the Energy Exascale Earth System Model version 1 (EAMv1) to address the previously reported issue of strong sensitivity of the simulated dust aerosol lifetime and dry removal rate to the model's vertical resolution. This paper complements that empirical justification of the revised scheme with a mathematical justification leveraging a semi-discrete analysis framework for assessing the splitting error of process coupling methods. The framework distinguishes the error due to numerical splitting from the error due to the time integration method(s) used for each individual process. Such a distinction results in a framework that provides an intuitive understanding of the causes of the splitting error. The application of this framework to the dust life cycle in EAMv1 confirms (i) that the original EAMv1 scheme artificially strengthens the effect of dry removal processes and (ii) that the revised splitting reduces that artificial strengthening. While the error analysis framework is presented in the context of the dust life cycle in EAMv1, the framework can be broadly leveraged to evaluate process coupling schemes, both in other physical problems and for any number of processes. This framework will be particularly powerful when the various process implementations support a variety of time integration approaches. Whereas traditional local truncation error approaches require separate consideration of each combination of time integration methods, this framework enables evaluation of coupling schemes independent of particular time integration approaches for each process while still allowing for the incorporation of these specific time integration errors if so desired. The framework also explains how the splitting error terms result from (i) the integration of individual processes in isolation from other processes and (ii) the choices of input state and time step size for the isolated integration of processes. Such a perspective has the potential for the rapid development of alternative coupling approaches that utilize knowledge both about the desired accuracy and about the computational costs of individual processes.

Higher‐order inverse mass matrices for the explicit transient analysis of heterogeneous solids

AbstractNew methods are presented for the direct computation of higher‐order inverse mass matrices (also called reciprocal mass matrices) that are used for explicit transient finite element analysis. The motivation of this work lies in the need of having appropriate sparse inverse mass matrices, which present the same structure as the consistent mass matrix, preserve the total mass, predict suitable frequency spectrum and dictate sufficiently large critical time step sizes. For an efficient evaluation of the reciprocal mass matrix, the projection matrix should be diagonal. This condition can be satisfied by adopting dual shape functions for the momentum field, generated from the same shape functions used for the displacement field. A theoretically consistent derivation of the inverse mass matrix is based on the three‐field Hamilton principle and requires the projection matrix to be evaluated from the integral of these shape functions. Unfortunately, for higher‐order FE shape functions and serendipity FE elements, the projection matrix is not positive definitive and can not be employed. Therefore, we study several lumping procedures for higher order reciprocal mass matrices considering their effect on total‐mass preserving, frequency spectra and accuracy in explicit transient simulations. The article closes with several numerical examples showing suitability of the direct inverse mass matrix in dynamics.

On the non-global linear stability and spurious fixed points of MPRK schemes with negative RK parameters

Recently, a stability theory has been developed to study the linear stability of modified Patankar–Runge–Kutta (MPRK) schemes. This stability theory provides sufficient conditions for a fixed point of an MPRK scheme to be stable as well as for the convergence of an MPRK scheme towards the steady state of the corresponding initial value problem, whereas the main assumption is that the initial value is sufficiently close to the steady state. Initially, numerical experiments in several publications indicated that these linear stability properties are not only local but even global, as is the case for general linear methods. Recently, however, it was discovered that the linear stability of the MPDeC(8) scheme is indeed only local in nature. Our conjecture is that this is a result of negative Runge–Kutta (RK) parameters of MPDeC(8) and that linear stability is indeed global if the RK parameters are nonnegative. To support this conjecture, we examine the family of MPRK22(α\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\varvec{\\alpha }$$\\end{document}) methods with negative RK parameters and show that even among these methods there are methods for which the stability properties are only local. However, this local linear stability is not observed for MPRK22(α\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\alpha $$\\end{document}) schemes with nonnegative Runge–Kutta parameters. In particular, it is shown that MPRK22(α\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\alpha $$\\end{document}) schemes with 0<α<0.5\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0<\\alpha <0.5$$\\end{document} or -0.5<α<0\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$-0.5<\\alpha <0$$\\end{document} are only stable if the time step size is sufficiently small. But schemes with α≤-0.5\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\alpha \\le -0.5$$\\end{document} are stable and converge towards the steady state of the initial value problems for all time step sizes, at least for the test problem under consideration. Furthermore, it is shown that for some of the latter systems, the initial values must actually be close enough to the steady state to guarantee this result.

Several Numerical Simulation Methods for the Time-Dependent Schrödinger Equation

The time-dependent Schrödinger equation can be solved using the Houbolt difference scheme or the space-time polynomial particular solutions method, with the former performing well in dissipative problems and the latter being suitable for wave-like problems. To handle variations in time-step sizes while overcoming numerical oscillations in non-smooth regions, these two methods are coupled. The time-domain is handled using the Houbolt difference format, while the spatial domain is handled using the method of polynomial particular solutions. The numerical results show that the advantages of each method are preserved. Additionally, this hybrid approach can be replaced with other difference methods and meshless methods to meet various numerical solution needs.

A Hu-Washizu variational approach to self-stabilized quadrilateral Virtual Elements: 2D linear elastodynamics

A recent mixed formulation of the Virtual Element Method in 2D elastostatics, based on the Hu-Washizu variational principle, is here extended to 2D elastodynamics. The independent modeling of the strain field, allowed by the mixed formulation, is exploited to derive first order quadrilateral Virtual Elements (VEs) not requiring a stabilization (namely, self-stabilized VEs), in contrast to the standard VEs, where an artificial stabilization is always required for first order quads. Lumped mass matrices are derived using a novel approach, based on an integration scheme that makes use of nodal values only, preserving the correct mass in the case of rigid-body modes. In the case of implicit time integration, it is shown how the combination of a self-stabilized stiffness matrix with a self-stabilized lumped mass matrix can produce excellent performances both in the compressible and quasi-incompressible regimes with almost negligible sensitivity to element distortion. Finally, in the case of explicit dynamics, the performances of the different types of derived VEs are analyzed in terms of their critical time-step size.

CFD Analysis for Valve-Holding Camber Permanent Inhaler Spacer (AerospaAcer) with Different Valves

During the COVID-19 pandemic, data statistics showed that patients with respiratory problems become infected. One of the therapy techniques for the respiratory condition was the use of a metered-dose inhaler and spacer. The objective of this paper is to determine the flow parameters of three types of valves which is duckbill valve, cross slit valve and umbrella valve in inhaler spacer to compare fluid flow between valve Previous researchers chose the duckbill valve to control fluid flow in inhaler spacer. The flow characteristics are unaffected by the materials used in the new disposable inhaler spacers, such as paper and polylactic acid (PLA). Several design valves reduced the skewness below 0.94 by suppressing the fillet and chamfer. ANSYS Workbench Fluent 19.2 is used to calculate flow parameters such as turbulence kinetic energy (TKE), turbulence eddy dissipation (TED), velocity, particle velocity magnitude, streamline, and vector velocity. The setup input data is based on the previous researcher's specified parameters such as viscosity model, drug characteristics (salbutamol and propellant), discrete phase model (DPM) equal to 80, boundary condition model, and SIMPLE technique. For the three types of valves, the nozzle injection used a 0.50-millimetre dimension. The simulation work is cross-checked against the results of prior simulations. Within each iteration, the transient flow employed a time step size of 0.01 for 200 steps. The results show that computational analysis can distinguish between models of varying complexity. The TED, TKE, and velocity graphs showed the approximate value between the model geometries. Overall, the study was successful in achieving the desired velocity magnitude in terms of visual and graph representations of the various valve.

A collision-based hybrid method for the BGK equation

We apply the collision-based hybrid method introduced by Hauck and McClarren [1] to the Boltzmann equation with the BGK operator and a hyperbolic scaling. An implicit treatment of the source term is used to handle stiffness associated with the BGK operator. Although it helps the numerical scheme become stable with a large time step size, it is still not obvious to achieve the desired order of accuracy due to the relationship between the size of the spatial cell and the mean free path. Without asymptotic preserving property, a very restricted grid size is required to resolve the mean free path, which is not practical. Our approaches are based on the noncollision-collision decomposition of the BGK equation. We introduce the arbitrary order of nodal discontinuous Galerkin (DG) discretization in space with a semi-implicit time-stepping method; we employ the backward Euler time integration for the uncollided equation and the 2nd order predictor-corrector scheme for the collided equation, i.e., both source terms in uncollided and collided equations are treated implicitly and only streaming term in the collided equation is solved explicitly. This improves the computational efficiency without the complexity of the numerical implementation. Numerical results are presented for various Knudsen numbers to present the effectiveness and accuracy of our hybrid method. Also, we compare the solutions of the hybrid and non-hybrid schemes.

Energy-Preserving/Group-Preserving Schemes for Depicting Nonlinear Vibrations of Multi-Coupled Duffing Oscillators

In the paper, we first develop a novel automatically energy-preserving scheme (AEPS) for the undamped and unforced single and multi-coupled Duffing equations by recasting them to the Lie-type systems of ordinary differential equations. The AEPS can automatically preserve the energy to be a constant value in a long-term free vibration behavior. The analytical solution of a special Duffing–van der Pol equation is compared with that computed by the novel group-preserving scheme (GPS) which has fourth-order accuracy. The main novelty is that we constructed the quadratic forms of the energy equations, the Lie-algebras and Lie-groups for the multi-coupled Duffing oscillator system. Then, we extend the GPS to the damped and forced Duffing equations. The corresponding algorithms are developed, which are effective to depict the long term nonlinear vibration behaviors of the multi-coupled Duffing oscillators with an accuracy of O(h4) for a small time stepsize h.

A positivity-preserving, linear, energy stable and convergent numerical scheme for the Poisson–Nernst–Planck (PNP) system

This article focuses on the convergence analysis for a fully discrete finite difference scheme for the time-dependent Poisson–Nernst–Planck system. The numerical scheme, a three-level linearized finite difference algorithm based on the reformulation of the Nernst–Planck equations, which was developed in [J. Sci. Comput. 81(2019),436–458]. The positivity-preserving property and energy stability were theoretically established. In this paper, we rigorously prove first-order convergence in time and second-order convergence in space for the numerical scheme in the ℓ∞(0,T;ℓ2)∩ℓ2(0,T;Hh1) norm under the condition that time step size is linearly proportional to spatial mesh size, in which the natural property of the exponential terms gives a lot of challenges. Moreover, the higher-order asymptotic expansion (up to third-order temporal accuracy and fourth-order spatial accuracy) has to be involved, due to the leading local truncation error will not be enough to recover an ℓ∞ bound for ion concentrations n and p. To our knowledge, this scheme will be the first linear decoupled algorithm to combine the following theoretical properties for the PNP system: ion concentration positivity preserving, unconditionally energy stability and optimal rate convergence. Numerical results are shown to be consistent with theoretical analysis.

A higher order stable numerical approximation for time‐fractional non‐linear Kuramoto–Sivashinsky equation based on quintic B$$ \mathfrak{B} $$‐spline

This article deals with designing and analyzing a higher order stable numerical analysis for the time‐fractional Kuramoto–Sivashinsky (K‐S) equation, which is a fourth‐order non‐linear equation. The fractional derivative of order present in the considered problem is taken into Caputo sense and approximated using the scheme. In space direction, the discretization process uses quintic ‐spline functions to approximate the derivatives and the solution of the problem. The present approach is unconditionally stable and is convergent with rate of accuracy , where and denote the space and time step sizes, respectively. We have also noted that the linearized version of the K‐S equation leads the rate of accuracy to . The present approach is also highly effective for the time‐fractional Burgers' equation. We have shown that the present approach provides better accuracy than the scheme with the same computational cost for several linear/non‐linear problems, with classical as well as fractional time derivatives.

Accuracy Analysis for Explicit-Implicit Finite Volume Schemes on Cut Cell Meshes

The solution of time-dependent hyperbolic conservation laws on cut cell meshes causes the small cell problem: standard schemes are not stable on the arbitrarily small cut cells if an explicit time stepping scheme is used and the time step size is chosen based on the size of the background cells. In May and Berger (J Sci Comput 71: 919–943, 2017), the mixed explicit-implicit approach in general and MUSCL-Trap (monotonic upwind scheme for conservation laws and trapezoidal scheme) in particular have been introduced to solve this problem by using implicit time stepping on the cut cells. Theoretical and numerical results have indicated that this might lead to a loss in accuracy when switching between the explicit and implicit time stepping. In this contribution, we examine this in more detail and will prove in one dimension that the specific combination MUSCL-Trap of an explicit second-order and an implicit second-order scheme results in a fully second-order mixed scheme. As this result is unlikely to hold in two dimensions, we also introduce two new versions of mixed explicit-implicit schemes based on exchanging the explicit scheme. We present numerical tests in two dimensions where we compare the new versions with the original MUSCL-Trap scheme.

Two-grid finite element method with an H2N2 interpolation for two-dimensional nonlinear fractional multi-term mixed sub-diffusion and diffusion wave equation

&lt;abstract&gt;&lt;p&gt;In this paper, we studied the two-grid method (TGM) for two-dimensional nonlinear time fractional multi-term mixed sub-diffusion and diffusion wave equation. A fully discrete scheme with the quadratic Hermite and Newton interpolation (H2N2) method was considered in the temporal direction and the expanded finite element method is used to approximate the spatial direction. In order to reduce computational time, a dual grid method based on Newton iteration was constructed with order $ \alpha\in(0, 1) $ and $ \beta\in(1, 2) $. The global convergence order of the two-grid scheme reaches $ O(\tau^{3-\beta}+h^{r+1}+H^{2r+2}) $, where $ \tau $, $ H $ and $ h $ are the time step size, coarse grid mesh size and fine grid mesh size, respectively. The error estimation and stability of the fully discrete scheme were derived. Theoretical analysis shows that the two grid algorithms maintain asymptotic optimal accuracy while saving computational costs. In addition, numerical experiments further confirmed the theoretical results.&lt;/p&gt;&lt;/abstract&gt;