In this paper, we are concerned with the numerical solutions of the coupled fractional Klein-Gordon-Schrödinger equation. The numerical schemes are constructed by combining the Crank-Nicolson/leap-frog difference methods for the temporal discretization and the Galerkin finite element methods for the spatial discretization. We give a detailed analysis of the conservation properties in the senses of discrete mass and energy. Then the numerical solutions are shown to be unconditionally bounded in L2 −norm, Hα2−\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$H^{\frac {\alpha }{2}}-$\end{document}semi-norm and L∞−\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$L^{\infty }-$\end{document}norm, respectively. Based on the well-known Brouwer fixed-point theorem and the mathematical induction, the unique solvability of the discrete solutions is proved. Moreover, the schemes are proved to be unconditionally convergent with the optimal order Oτ2+hr+1\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$O\left (\tau ^{2}+h^{r+1}\right )$\end{document}, where τ is the temporal step, h is the spatial grid size, and r is the order of the selected finite element space. Furthermore, by using the proposed decoupling and iterative algorithms, several numerical examples are included to support theoretical results and show the effectiveness of the schemes. Finally, the fast Krylov subspace solver with suitable circulant preconditioner is designed to effectively solve the Toeplitz-like linear systems. In each iterative step, this method can effectively reduce the memory requirement of above each finite element scheme from OM2\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}${{O}\left (M^{2}\right )}$\end{document} to O(M), and the computational complexity from OM3\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}${O\left (M^{3}\right )}$\end{document} to O(MlogM)\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}${O(M \log M)}$\end{document}, where M is the number of grid nodes. Numerical tests are carried out to show that this fast algorithm is more practical than the traditional backslash and LU factorization/Cholesky decomposition methods, in terms of memory requirement and computational cost.