The key aim of this paper is to provide a new nth generalized order perturbed isogeometric fast multistage technique of boundary elements to compute the propagation of time harmonics in an infinite region. Structural geometry and boundary integral equations are constructed by using non-uniform rational B-splines. The source of system uncertainty is believed to be the incident plane wave number’s unpredictability. The actual field, depending on the input random variables, is simulated using the extended nth-order perturbation method. The field and kernel values for boundary integral formulas are generated via the nth-order generalized series of Taylor expansions using perturbation parameters. The fast multipole method (FMM) is utilized to speed up the process. The effectiveness and correctness of the proposed algorithm are verified by Monte Carlo simulations (MCs) with numerical examples.