Abstract: To improve the efficiency of the aerodynamic shape optimization based on the discrete adjoint method in a wide speed range, a hybrid time integration method is proposed. The hybrid method consists of the lower-upper symmetric Gauss-Seidel (LU-SGS) scheme, approximate Newton-Krylov (ANK) scheme, and Newton-Krylov (NK) scheme. The LU-SGS scheme is first employed to start the iteration, and then the ANK and NK schemes are sequentially used to accelerate the convergence rate. The relative convergence of the residual norm is employed to switch the three methods above. The ANK method is selected until the residual norm decreases by 2–4 orders of magnitude. The NK method is specified when the residual norm is 10² to 10³ times the target residual. By this switching strategy, the hybrid scheme is able to accelerate the convergence rate and improve the efficiency of aerodynamic optimization. The results of the test cases demonstrate that the solver employing the LU-SGS+ANK+NK scheme is 75% faster than that using the LU-SGS scheme. Moreover, compared to the Runge-Kutta (RK)+ANK+NK and diagonalized diagonally dominant alternating direction implicit (D3ADI)+ANK+NK schemes, the LU-SGS+ANK+NK scheme exhibits greater robustness. Furthermore, in a wide-speed-range aerodynamic optimization case, the optimizer based on the hybrid integration scheme is about 70% faster than that based on the LU-SGS method by accelerating each flow field simulation.