Search

In this paper, we intend to utilize the generalized scalar auxiliary variable (GSAV) approach proposed in recent paper (Ju et al., SIAM J. Numer. Anal., 60 (2022), 1905–1931) for the nonlocal coupled sine-Gordon equation to construct a class of fully decoupled, linear, and second-order accurate energy-preserving scheme. The unconditional unique solvability and discrete energy conservation law of the proposed scheme are rigorously discussed, and the unconditional convergence is then proved by the mathematical induction. Particularly, the convergence analysis covers the proposed scheme in multiple dimensions due to the corresponding nonlinear terms satisfy the global Lipschitz continuity straightforwardly. Finally, time evolution of dynamical behavior of the governing equation with different nonlocal parameters are observed, and ample numerical comparisons demonstrate that the proposed scheme manifests high efficiency in long-time computations. [ABSTRACT FROM AUTHOR]

This-36pt paper considers a discontinuous Galerkin finite element method for the 2D transient and incompressible Navier-Stokes model. Following the analysis of Heywood and Rannacher (SIAM J. Numer. Anal. 19:275–311, 1982), we derive optimal velocity and pressure error estimates in L ∞ (L 2) and L ∞ (L 2) -norms, respectively, for the discontinuous Galerkin case. We use standard L 2 -projection and modified Stokes operator but on appropriate broken Sobolev spaces, and then standard duality arguments to achieve these results. For sufficiently small data, uniform in time estimates are proved. Based on the backward Euler method, time discretization is carried out and fully discrete error estimates are derived. Finally, numerical experiments are conducted to verify our theoretical findings. [ABSTRACT FROM AUTHOR]

We consider operator preconditioning B − 1 A , which is employed in the numerical solution of boundary value problems. Here, the self-adjoint operators A , B : H 0 1 (Ω) → H − 1 (Ω) are the standard integral/functional representations of the partial differential operators −∇⋅ (k(x)∇u) and −∇⋅ (g(x)∇u), respectively, and the scalar coefficient functions k(x) and g(x) are assumed to be continuous throughout the closure of the solution domain. The function g(x) is also assumed to be uniformly positive. When the discretized problem, with the preconditioned operator B n − 1 A n , is solved with Krylov subspace methods, the convergence behavior depends on the distribution of the eigenvalues. Therefore, it is crucial to understand how the eigenvalues of B n − 1 A n are related to the spectrum of B − 1 A . Following the path started in the two recent papers published in SIAM J. Numer. Anal. [57 (2019), pp. 1369-1394 and 58 (2020), pp. 2193-2211], the first part of this paper addresses the open question concerning the distribution of the eigenvalues of B n − 1 A n formulated at the end of the second paper. The approximation of the spectrum studied in the present paper differs from the eigenvalue problem studied in the classical PDE literature which addresses individual eigenvalues of compact (solution) operators. In the second part of this paper, we generalize some of our results to general bounded and self-adjoint operators A , B : V → V # , where V# denotes the dual of V. More specifically, provided that B is coercive and that the standard Galerkin discretization approximation properties hold, we prove that the whole spectrum of B − 1 A : V → V is approximated to an arbitrary accuracy by the eigenvalues of its finite dimensional discretization B n − 1 A n . [ABSTRACT FROM AUTHOR]

In a very interesting paper (SIAM J. Control Optim. 37(3): 765–776, 1999), Solodov and Svaiter introduced an effective projection algorithm with linesearch for finding a solution of a variational inequality problem (VIP) in Euclidean space. They showed that the iterative sequence generated by their algorithm converges to a solution of (VIP) under the main assumption that the cost mapping is pseudomonotone and continuous. In this paper, we propose to modify this algorithm for solving variational inequality problems in which the cost mapping is not required to be satisfied any pseudomonotonicity. Moreover, we do not use the embedded projection methods as in methods used in literature and the linesearch procedure is not necessary when the cost mapping is Lipschitz. Several numerical examples are also provided to illustrate the efficient of the proposed algorithms. [ABSTRACT FROM AUTHOR]

In this paper, four accelerated subspace minimization conjugate gradient methods based on 2-regularization or 3-regularization models with different matrices are proposed. The search directions are generated by minimizing a quadratic approximation and different regularization models of the objective function, and the new directions satisfy the sufficient descent condition. In addition, an acceleration technique is used to improve the stepsize. With a modified nonmonotone line search, we establish the global convergence of the proposed methods under mild assumptions. By using the Kurdyka-Łojasiewicz property, we analyze the convergence rates of function sequence and iterative sequence without the convexity of objective function. Numerical experiments show that the proposed methods are very effective, especially for CUTEr test problems with no less than 50 dimensions, the best of the new methods is competitive to the latest limited memory conjugate gradient software package CG_DESCENT(6.8) developed by Hager and Zhang (SIAM J. Optim. 23:2150–2168, 24). [ABSTRACT FROM AUTHOR]

In this paper, we use a fully discrete Crank-Nicolson scheme to solve the Cauchy problem for one-dimensional linearized Green-Naghdi system with fast convolution boundary condition which is derived through the Padé approximation for the square root function. We also introduce a constant damping term to guarantee the stability. While the damping term meets certain conditions, the stability is proved for the numerical solutions with the fast convolution boundary condition. Therefore, the difficulty of numerical instability which rises in Kazakova and Noble (SIAM J. Num. Anal. 1, 657–683, 2020) is overcome. The computational cost of the convolution integral is reduced from O (N 2) to O (N ln (N)) for the total number of time steps N. Numerical examples verify the theoretical results and demonstrate the performance for the fast numerical method. [ABSTRACT FROM AUTHOR]

In this paper, we propose a divergence-free nonconforming virtual element discrete scheme for (Navier-) Stokes problem based on Raviart-Thomas-like virtual element space. By choosing different (compared to Zhao et al. (SIAM J. Numer. Anal. 57, 2730–2759, 2019)) degrees of freedom and global virtual element space, we realize H1- and L2-error estimates of the velocity that are independent of the continuous pressure. Moreover, the presented scheme can also deal with polygonal meshes (including non-convex and degenerate elements), arbitrary approximation orders k, and large Reynolds numbers. Finally, we investigate several classical numerical experiments of the incompressible flow to present the performance of this numerical scheme. [ABSTRACT FROM AUTHOR]

Recently, the condition numbers of the total least squares (TLS) problems having a unique solution have been studied at length in Zheng et al. (SIAM J. Matrix Anal. Appl. 38: 924–948, 2017). However, it is known that the TLS problem may have no solution, and even if an existing solution, it may not be unique. As a continuation of their work, in this paper, we investigate the condition numbers of the minimum Frobenius norm solution of the (multidimensional) TLS problem when having more than one solution. The tight and computable upper bound estimates of the normwise, mixed, and componentwise condition numbers are respectively derived. Some numerical experiments are performed to illustrate the tightness of these upper bounds. [ABSTRACT FROM AUTHOR]

In this paper, we present an adaptive variable step size numerical scheme for the integration of linear stochastic oscillator equations driven by additive Brownian white noise. We first show that traditional adaptive schemes based on local error estimation destroy the long-time behavior of the underlying method. As a remedy, we extend the idea presented in Hairer and Söderlind (SIAM J. Sci. Comput. 26(6), 1838-1851 2005) to the stochastic setting and show that using step density control mechanisms based on time regularization and local error tracking, we are able to obtain numerical schemes which preserve the important qualitative features of the solution process such as symmetry, time reversibility, symplecticity, linear growth rate of the second moment, and infinite oscillation. Numerical experiments confirm the theoretical findings of the paper. [ABSTRACT FROM AUTHOR]

In this paper, two new subspace minimization conjugate gradient methods based on p-regularization models are proposed, where a special scaled norm in p-regularization model is analyzed. Different choices of special scaled norm lead to different solutions to the p-regularized subproblem. Based on the analyses of the solutions in a two-dimensional subspace, we derive new directions satisfying the sufficient descent condition. With a modified nonmonotone line search, we establish the global convergence of the proposed methods under mild assumptions. R-linear convergence of the proposed methods is also analyzed. Numerical results show that, for the CUTEr library, the proposed methods are superior to four conjugate gradient methods, which were proposed by Hager and Zhang (SIAM J. Optim. 16(1):170–192, 2005), Dai and Kou (SIAM J. Optim. 23(1):296–320, 2013), Liu and Liu (J. Optim. Theory. Appl. 180(3):879–906, 2019) and Li et al. (Comput. Appl. Math. 38(1):2019), respectively. [ABSTRACT FROM AUTHOR]

In this paper, we combine the S-iteration process introduced by Agarwal et al. (J. Nonlinear Convex Anal., 8(1), 61–79 2007) with the proximal point algorithm introduced by Rockafellar (SIAM J. Control Optim., 14, 877–898 1976) to propose a new modified proximal point algorithm based on the S-type iteration process for approximating a common element of the set of solutions of convex minimization problems and the set of fixed points of nearly asymptotically quasi-nonexpansive mappings in the framework of CAT(0) spaces and prove the △-convergence of the proposed algorithm for solving common minimization problem and common fixed point problem. Our result generalizes, extends and unifies the corresponding results of Dhompongsa and Panyanak (Comput. Math. Appl., 56, 2572–2579 2008), Khan and Abbas (Comput. Math. Appl., 61, 109–116 2011), Abbas et al. (Math. Comput. Modelling, 55, 1418–1427 2012) and many more. [ABSTRACT FROM AUTHOR]

B-method is a novel method developed by Beck et al. (SIAM J. Sci. Comput. 37(5), A2998–A3029, 2015), and has been shown theoretically to be very advantageous in time discretization of the second-order parabolic equations with blow-up solutions. In this paper, we extend the B-method to approximate the blow-up solution of a class of fourth-order parabolic equations, which plays very important role in many engineering applications. First, by following the systematic means of constructing numerical schemes based on the technique of variation of constants proposed by Beck et al., we give some B-method schemes for the fourth-order semilinear parabolic equations. Second, we perform a truncation error analysis to show when and why the B-method scheme is advantageous over its classical counterpart. Third, we take one of the constructed numerical schemes as an example to show the well-posedness using the technique of upper and lower solutions. Last, we carry out numerical experiments to approximate the blow-up solutions and illustrate the efficiency of our numerical schemes. [ABSTRACT FROM AUTHOR]

With the recent emergence of mixed precision hardware, there has been a renewed interest in its use for solving numerical linear algebra problems fast and accurately. The solution of total least squares problems, i.e., solving min E , r ‖ [ E , r ] ‖ F subject to (A + E) x = b + r , arises in numerous applications. Solving this problem requires finding the smallest singular value and corresponding right singular vector of [ A , b ] , which is challenging when A is large and sparse. An efficient algorithm for this case due to Björck et al. (SIAM J. Matrix Anal. Appl. 22(2), 413–429 2000), called RQI-PCGTLS, is based on Rayleigh quotient iteration coupled with the preconditioned conjugate gradient method. We develop a mixed precision variant of this algorithm, RQI-PCGTLS-MP, in which up to three different precisions can be used. We assume that the lowest precision is used in the computation of the preconditioner and give theoretical constraints on how this precision must be chosen to ensure stability. In contrast to standard least squares, for total least squares, the resulting constraint depends not only on the matrix A , but also on the right-hand side b . We perform a number of numerical experiments on model total least squares problems used in the literature, which demonstrate that our algorithm can attain the same accuracy as RQI-PCGTLS albeit with a potential convergence delay due to the use of low precision. Performance modeling shows that the mixed precision approach can achieve up to a 4 × speedup depending on the size of the matrix and the number of Rayleigh quotient iterations performed. [ABSTRACT FROM AUTHOR]

A randomized quaternion singular value decomposition algorithm based on block Krylov iteration (RQSVD-BKI) is presented to solve the low-rank quaternion matrix approximation problem. The upper bounds of deterministic approximation error and expected approximation error for the RQSVD-BKI algorithm are also given. It is shown by numerical experiments that the running time of the RQSVD-BKI algorithm is smaller than that of the quaternion singular value decomposition, and the relative errors of the RQSVD-BKI algorithm are smaller than those of the randomized quaternion singular value decomposition algorithm in Liu et al. (SIAM J. Sci. Comput., 44(2): A870-A900 (2022)) in some cases. In order to further illustrate the feasibility and effectiveness of the RQSVD-BKI algorithm, we use it to deal with the problem of color image inpainting. [ABSTRACT FROM AUTHOR]

This work is devoted to show the efficiency of a new numerical approach in solving geometrical shape optimization problems constrained to partial differential equations, on a family of convex domains. More precisely, we are interested to an improved numerical optimization process based on the new shape derivative formula, using the Minkowski deformation of convex domains, recently established in Boulkhemair and Chakib (J. Convex Anal. 21( n ∘ 1 ), 67–87 2014), Boulkhemair (SIAM J. Control Optim. 55( n ∘ 1 ), 156–171 2017). This last formula allows to express the shape derivative by means of the support function, in contrast to the classical one expressed in term of vector fields Henrot and Pierre 2005, Delfour and Zolésio 2011, Sokolowski and Zolesio 1992. This avoids some of the disadvantages related to the classical shape derivative approach, when one use the finite elements discretization for approximating the auxiliary boundary value problems in shape optimization processes Allaire 2007. So, we investigate here the performance of the proposed shape optimization approach through the numerical resolution of some shape optimization problems constrained to boundary value problems governed by Laplace or Stokes operator. Notably, we carry out a comparative numerical study between its resulting numerical optimization process and the classical one. Finally, we give some numerical results showing the efficiency of the proposed approach and its ability in producing good quality solutions and in providing better accuracy for the optimal solution in less CPU time compared to the classical approach. [ABSTRACT FROM AUTHOR]

In this paper, we are concerned with the numerical treatment of a recent diffuse interface model for two-phase flow of electrolyte solutions (Campillo-Funollet et al., SIAM J. Appl. Math. 72(6), 1899-1925, 2012). This model consists of a Nernst-Planck-system describing the evolution of the ion densities and the electrostatic potential which is coupled to a Cahn-Hilliard-Navier-Stokes-system describing the evolution of phase-field, velocity field, and pressure. In the first part, we present a stable, fully discrete splitting scheme, which allows to split the governing equations into different blocks, which may be treated sequentially and thereby reduces the computational costs significantly. This scheme comprises different mechanisms to reduce the induced numerical dissipation. In the second part, we investigate the impact of these mechanisms on the scheme's sensitivity to the size of the time increment using the example of a falling droplet. Finally, we shall present simulations showing ion induced changes in the topology of charged droplets serving as a qualitative validation for our discretization and the underlying model. [ABSTRACT FROM AUTHOR]

Recently, an augmented IIM (Li et al. SIAM J. Numer. Anal. 55, 570-597 2017) has been developed and analyzed for some interface problems. The augmented IIM can provide not only a second-order accurate solution in the entire domain but also second-order accurate gradients from each side of the interface. In the augmented IIM, the PDE is reformulated and an augmented variable of co-dimension one is introduced and solved through a Schur complement system. Nevertheless, the augmented IIM is somewhat difficult to implement and understand for non-experts in the area. In this paper, a direct method (thus simpler than augmented IIM) without using the augmented variable is proposed for elliptic interface problems with a variable coefficient that can have a finite jump across an interface. The resulting finite difference scheme is the standard five-point central scheme at regular grid points, while it is a compact nine-point scheme at irregular grid points. The computed solution is second-order accurate and will be used to recover the gradient from each side of the interface to second-order accuracy. The discrete Green's function is constructed and used to prove second-order convergence of both the solution and gradient for the one-dimensional algorithm with a piecewise constant coefficient. For the two-dimensional algorithm with a piecewise constant coefficient, we numerically prove second-order convergence of the solution by applying the discrete elliptic maximum principle using an optimization solver, while the second-order convergence of the gradient and the same conclusion for more general problems will be only demonstrated in the numerical examples. [ABSTRACT FROM AUTHOR]

Distributed algorithms are receiving renewed attention across multiple disciplines due to the dramatically increasing demand of big data processing. We consider a class of consensus optimization problems over a static network system of multiple agents, where each of local cost functions is a sum of a smooth function and a non-Lipschitz regularization term. This kind of problems is widely found in scientific and engineering areas such as machine learning and data analysis. Inspired by Bian and Chen (SIAM J. Opt. 23(3), 1718–1741 2017), we propose a decentralized smoothing quadratic regularization algorithm (abbreviated as D-SQRA) for solving the composite consensus problem with non-Lipschitz singularities. To some extent, D-SQRA can be seen as an extension in the decentralized setup of the smoothing quadratic regularization algorithm (SQRA) proposed in (SIAM J. Opt. 23(3), 1718–1741 2017). Our main contribution is to show that D-SQRA can inherit the theoretical properties of its centralized counterpart, i.e., SQRA, in both sides of convergence and worst-case iteration complexity to achieve an ϵ scaled stationary point. We also present some numerical examples on sparse sensing problems based on synthetic and real datasets to corroborate the effectiveness of the proposed decentralized algorithm. [ABSTRACT FROM AUTHOR]

We analyze two algorithms for computing the symplectic factorization A = LLT of a given symmetric positive definite symplectic matrix A. The first algorithm W1 is an implementation of the HHT factorization from Dopico and Johnson (SIAM J. Matrix Anal. Appl. 31(2):650–673, 2009), see Theorem 5.2. The second one is a new algorithm W2 that uses both Cholesky and Reverse Cholesky decompositions of symmetric positive definite matrices. We present a comparison of these algorithms and illustrate their properties by numerical experiments in MATLAB. A particular emphasis is given on symplecticity properties of the computed matrices in floating-point arithmetic. [ABSTRACT FROM AUTHOR]

The asymptotic stable region and long-time decay rate of solutions to linear homogeneous Caputo time fractional ordinary differential equations (F-ODEs) are known to be completely determined by the eigenvalues of the coefficient matrix. Very different from the exponential decay of solutions to classical ODEs, solutions of F-ODEs decay only polynomially, leading to the so-called Mittag-Leffler stability, which was already extended to semi-linear F-ODEs with small perturbations. This work is mainly devoted to the qualitative analysis of the long-time behavior of numerical solutions. By applying the singularity analysis of generating functions developed by Flajolet and Odlyzko (SIAM J. Disc. Math. 3(2), 216–240, 1990), we are able to prove that both L 1 scheme and strong A-stable fractional linear multistep methods (F-LMMs) can preserve the numerical Mittag-Leffler stability for linear homogeneous F-ODEs exactly as in the continuous case. Through an improved estimate of the discrete fractional resolvent operator, we show that strong A-stable F-LMMs are also Mittag-Leffler stable for semi-linear F-ODEs under small perturbations. For the numerical schemes based on α -difference approximation to Caputo derivative, we establish the Mittag-Leffler stability for semi-linear problems by making use of properties of the Poisson transformation. The new results and analyses provide not only the rigorous justifications and explanations of the Mittag-Leffler stability of numerical solutions with exact decay rate, but also establish some close connection between the continuous and discrete F-ODEs. Numerical experiments are presented for several typical time fractional evolutional equations, including time fractional sub-diffusion equations and semi-linear F-ODEs. All the numerical results exhibit the typical long-time polynomial decay rate, which is fully consistent with our theoretical predictions. [ABSTRACT FROM AUTHOR]

We develop a formal construction of a pointwise divergence-free basis in the nonconforming virtual element method of arbitrary order for the Stokes problem introduced in Zhao et al. (SIAM J. Numer. Anal. 57(6):2730–2759, 2019). The proposed construction can be seen as a generalization of the divergence-free basis in Crouzeix-Raviart finite element space (Brenner, Math. Comp. 55(192):411–437, 1990; Thomasset, 1981) to the virtual element space. Using the divergence-free basis obtained from our construction, we can eliminate the pressure variable from the mixed system and obtain a symmetric positive definite system. Several numerical tests are presented to confirm the efficiency and the accuracy of our construction. [ABSTRACT FROM AUTHOR]

A new computational framework for evaluation of the gamma function Γ(z) over the complex plane is developed. The algorithm is based on interpolation by rational functions, and generalizes the classical methods of Lanczos (SIAM J. Numer. Anal. B 1:86–96, (1964) and Spouge (SIAM J. Numer. Anal., 31(3):931–944, (1994) (which we show are also interpolatory). This framework utilizes the exact poles of the gamma function. By relaxing this condition and allowing the poles to vary, a near-optimal rational approximation is possible, which is demonstrated using the adaptive Antoulous Anderson (AAA) algorithm, developed in Nakatsukasa et al. (Appl. Math. Comp., 40:1494-1522, (2016)) and Nakatsukasa and Trefethen (SIAM J. Sci. Comp., 42(5):A3157–A3179, (2020)). The resulting approximations are competitive with Stirling's formula in terms of overall efficiency. [ABSTRACT FROM AUTHOR]

The k th Fréchet derivative of a matrix function f is a multilinear operator from a cartesian product of k subsets of the space ℂ n × n into itself. We show that the k th Fréchet derivative of a real-valued matrix function f at a real matrix A in real direction matrices E1, E2, ... , Ek can be computed using the complex step approximation. We exploit the algorithm of Higham and Relton (SIAM J. Matrix Anal. Appl.35(3):1019–1037, 2014) with the complex step approximation and mixed derivative of complex step and central finite difference scheme. Comparing with their approach, our cost analysis and numerical experiment reveal that half and seven-eighths of the computational cost can be saved for the complex step and mixed derivative, respectively. When f has an algorithm that computes its action on a vector, the computational cost drops down significantly as the dimension of the problem and k increase. [ABSTRACT FROM AUTHOR]

We study Leja sequences on the unit disc and resulting mapped sequences to conventional domain, ellipses and real intervals, for the problem of relaxation of Richardson iteration. Using simple considerations, we establish upper and lower estimates for the growth of associated Newton polynomials and the so-called residual polynomials. These results broaden the understanding of such sequences and add to results established in Calvetti et al. (Numer. Math. 67(1), 21–40, 1994), Calvetti and Reichel (Numerical Algorithms 11(1-4), 79–98, 1996; J. Comput. Appl. Math. 71(2), 267–286 1996), Fischer and Reichel (Numer. Math. 54(2), 225–242, 1989), Nachtigal et al. (SIAM Journal on Matrix Analysis and Applications 13(3), 796–825, 1992), Reichel (SIAM J. Numer. Anal. 25(6), 1359–1368, 1988), Tal-Ezer (J. Sci. Comput. 4(1), 25–60, 1989). We also propose adaptive strategies for the selection of relaxation parameters that build on fast strategies (Balgama et al., Electron. Trans. Numer. Anal. 7, 124–140, 1998; Reichel, Linear Algebra and its Applications 154-156(C), 389–414, 1991) which are known to be suitable for the problem. In particular a Leja ordering is enforced on the mapped sequences. New strategies for conveniently confining the spectrums of certain matrices arising in PDEs discretization is described and demonstrate the flexibility of the Richardson iteration. [ABSTRACT FROM AUTHOR]

Computing the greatest common divisor of a set of polynomials is a problem which plays an important role in different fields, such as linear system, control, and network theory. In practice, the polynomials are obtained through measurements and computations, so that their coefficients are inexact. This poses the problem of computing an approximate common factor. We propose an improvement and a generalization of the method recently proposed in Guglielmi et al. (SIAM J. Numer. Anal. 55, 1456–1482, 2017), which restates the problem as a (structured) distance to singularity of the Sylvester matrix. We generalize the algorithm in order to work with more than 2 polynomials and to compute an Approximate GCD (Greatest Common Divisor) of degree k ≥ 1; moreover, we show that the algorithm becomes faster by replacing the eigenvalues by the singular values. [ABSTRACT FROM AUTHOR]

In this work, we develop geometric multigrid algorithms for the immersed finite element methods for elliptic problems with interface (Chou et al. Adv. Comput. Math. 33, 149–168 2010; Kwak and Lee, Int. J. Pure Appl. Math. 104, 471–494 2015; Li et al. Numer. Math. 96, 61–98 2003, 2004; Lin et al. SIAM J. Numer. Anal. 53, 1121–1144 2015). We need to design the transfer operators between levels carefully, since the residuals of finer grid problems do not satisfy the flux condition once projected onto coarser grids. Hence, we have to modify the projected residuals so that the flux conditions are satisfied. Similarly, the correction has to be modified after prolongation. Two algorithms are suggested: one for finite element spaces having vertex degrees of freedom and the other for edge average degrees of freedom. For the second case, we use the idea of conforming subspace correction used for P1 nonconforming case (Lee 1993). Numerical experiments show the optimal scalability in terms of number of arithmetic operations, i.e., O (N) for V -cycle and CG algorithms preconditioned with V -cycle. In V -cycle, we used only one Gauss-Seidel smoothing. The CPU times are also reported. [ABSTRACT FROM AUTHOR]

In the present work, orthogonal spline collocation (OSC) method with convergence order O(τ3−α + hr+ 1) is proposed for the two-dimensional (2D) fourth-order fractional reaction-diffusion equation, where τ, h, r, and α are the time-step size, space size, polynomial degree of space, and the order of the time-fractional derivative (0 < α < 1), respectively. The method is based on applying a high-order finite difference method (FDM) to approximate the time Caputo fractional derivative and employing OSC method to approximate the spatial fourth-order derivative. Using the argument developed recently by Lv and Xu (SIAM J. Sci. Comput. 38, A2699-A2724, 2016) and mathematical induction method, the optimal error estimates of proposed fully discrete OSC method are proved in detail. Then, the theoretical analysis is validated by a number of numerical experiments. To the best of our knowledge, this is the first proof on the error estimates of high-order numerical method with convergence order O(τ3−α + hr+ 1) for the 2D fourth-order fractional equation. [ABSTRACT FROM AUTHOR]