Numerical Solutions of Nonlinear Schrödinger Equation: An Application Example of Nonlinear Analysis

Abstract

The nonlinear Schrödinger equation is used to show how numerical methods can be used to solve mathematical problems present in nonlinear analysis. The Lanzos-Chevbychev Pseudospectral method is shown to be effective, flexible, and economical to meet various demands in practical applications of mathematical simulations using nonlinear differential equations. The electromagnetic wave propagation through an inhomogeneous, anisotropic, and complex space is used as an example to show how successful mathematical modeling could be used to explain the complex phenomenon of astronomical redshift that is the central issue in the widely debated Hubble tension.

1. Introduction

In recent times, nonlinear analysis has been increasingly used in science and technology. Many advanced and innovative applications in those fields include nonlinear effects in their design and development. To be useful to real-world problems, those mathematical models need to be solved by methods developed in nonlinear analysis. Out of many possible mathematical methods, some are developed specifically for nonlinear differential equations (NDEs). This chapter will concentrate on a specific method for the solutions of second order NDEs. As a specific example, nonlinear Schrödinger equation (NLSE) is being chosen. The emphasis is to show how such numerical methods can be used to investigate how electromagnetic waves propagate under various realistic physical conditions in space.

For many years researchers have had extensive interest in how to solve NDEs analytically. But, because of the nonlinear nature, little success has been achieved in solving them directly. That is, starting from NDE itself and finding the solutions analytically in a forward direction. However, for methods starting from some assumed solutions and working out how to satisfy the NDE analytically as an inverse problem, there are many successes. Some of these inverse examples include the inverse differential and integral methods such as, for wave propagation [1], the G’/G expansion method [2] and its various variants [3, 4, 5], and inverse scattering methods for antenna design [6]. Generally, as an inverse problem, there is no limit to how many solutions can be found because there are an infinite number of choices for the set of system parameters that define the chosen base functions. However, the need to have a matching background medium is a notable limitation [7]. From the nature of those solutions, it could be concluded that this inverse approach is more suitable for qualitative analysis that the performance of a design, or the characteristics of a system could be assessed qualitatively. For quantitative assessment, the direct approach is a better choice, because the solutions are obtained by satisfying not just the NDE but also the initial/boundary conditions.

Numerical solutions of NDEs start with a scheme to discretize the problem into a set of simultaneous nonlinear algebraic equations. Linear algebraic algorithms are then used to solve those equations with an iterative scheme to cater for the nonlinear terms. For transient problems, especially when a long history involving a large set of equations is needed, the computational efficiency of the chosen method becomes important. As many different models under different prescribed conditions may be encountered, the flexibility of the method is also a factor for consideration.

In Section 2 of this chapter, we describe the Lanczos-Chevbychev Pseudospectral (LCPS) method [8, 9] that we have used to solve many different NDEs. The LCPS method uses an economized power series and has been shown to perform as well as similar orthogonal eigenfunction series expansion methods such as the Chebyshev Pseudospectral method. However, the advantage of using LCPS is that an ordinary power series is involved that would be the simplest and the most economical computing method. The details of the LCPS method are given in this section together with some application examples.

In Section 3, we apply NLSE to electromagnetic wave propagation through space [10, 11], together with two simple examples. We show in Section 4 that long distance, and other characteristic nature of space such as anisotropy, inhomogeneity, and gravitational effect, can be effectively included. How to calibrate our findings with empirical data is also described there. In Section 5, we discuss the usefulness and limitations of mathematical simulations based on examples we have solved. The check list includes items such as the appropriateness of the model, the variable ranges within which the model is applicable, and the implication of any assumptions made. To be realistic, we make use of our findings on astronomical redshift and compare them to popularly accepted theories in that the debate on Hubble tension is receiving considerable attention [12, 13]. We present our conclusions in Section 6.

2. Numerical solution methods for second order nonlinear partial differential equations

As higher order can be reduced to second order by introducing additional second order equations, we can restrict ourselves to consider only a second order differential equation in a dispersive field.

2.1 Nonlinear partial differential equations

Consider a time-dependent two-dimensional boundary problem

$$i{\mathbf{u}}_t=\mathbf{D}\left(\mathbf{x},t\right){\nabla }^2\mathrm{u}+\mathbf{F}\left(\mathbf{x},t,\mathbf{u}\right)\mathbf{u},$$E1

$$\alpha {\mathit{u}}^{\prime }\left(\mathit{x},t\right)+\mathrm{\beta }\mathit{u}\left(\mathit{x},\mathrm{t}\right)+\mathrm{\gamma }=0,\mathrm{at}\mathit{x}={\mathit{x}}_{\mathit{b}},$$E2

and,

$$\mathit{u}\left(\mathit{x},0\right)={\mathit{u}}_{\mathit{o}}\left(\mathit{x}\right),$$E3

where D is the dispersion coefficient, and the nonlinear F is a spatial and time-dependent potential. α, β, and γ are coefficients associated with the boundary conditions.

As an example, we use a single mode solution u(X,Y,t) in a two-dimensional Cartesian system with spatial variables, X and Y, coefficients D₁ and D_2, and a nonlinear potential F, such that

$$\mathit{x}=\left\{\begin{array}{l}X\\ Y\end{array}\right\},\mathit{D}\left(\mathit{x},t\right)=\left[\begin{array}{ll}{D}_1\left(X,Y,t\right)& 0\\ 0& D_2\left(X,Y,t\right)\end{array}\right],\mathit{F}\left(x,t,u\right)=\left\{\begin{array}{l}{F}_u\left(X,Y,t,u,v\right)\\ F_v\left(X,Y,t,u.v\right)\end{array}\right\}.$$E4

For numerical reasons, if a function is not smooth, such as in the case solitons, it is desirable to adopt a multidomain approach. The given rectangular two-dimensional domain of interest is divided into M x N subdomains. The affine transformation is used to scale each subdomain to [−1, 1]², in the new coordinates {x,y},

$$\Omega ≔{\Omega }^{i,j},i=1,2,\dots ,M;j=1,2,\dots ,N.$$E5

In the subdomains, the associated surfaces are

$$\begin{array}{lll}{S}_x^{0,j}& =& {\Omega }^{1,j}\left[-1,y\right],S_x^{i,j}={\Omega }^{i,j}\left[1,y\right]\cap {\Omega }^{i+1,j}\left[-1,y\right],i=1,2,\dots ,M-1,\\ S_x^{M,j}& =& {\Omega }^{M,j}\left[1,y\right],j=1,2,\dots ,N,\\ S_y^{i,0}& =& {\Omega }^{i,1}\left[x,-1\right],S_y^{i,j}={\Omega }^{i,j}\left[x,1\right]\cap {\Omega }^{i,j+1}\left[x,-1\right],j=1,2,\dots ,N-1,\\ S_y^{i,N}& =& {\Omega }^{i,N}\left[x,1\right],i=1,2,\dots ,M.\end{array}$$E6

Boundary conditions specified in Eqs. (2) and (3) apply only on those surfaces that form part of the boundary. For inter-subdomain surfaces, the specified conditions are continuities of both the function and its derivative. These also apply to the four corner points.

Based on the Lanczos-Chebushev Pseudospectral (LCS) method [8, 9], the function u^i,j in each subdomain Ω^i,j, is be represented by the tensor product of two truncated power series,

$$u^{i,j}\left(x,y,t\right)=\sum _{l=0}^L\sum _{k=0}^K{u}_{l,k}^{i,j}\left(t\right)x^k{y}^l$$E7

and, using term-by-term differentiation, the derivatives

$${\left(u_x\right)}^{i,j}\left(x,y,t\right)=\sum _{l=0}^L\sum _{k=1}^{K}k\underset{l,k}{\overset{i,j}{u}}\left(t\right)x^{k-1}{y}^l,$$

$${\left(u_y\right)}^{i,j}\left(x,y,t\right)=\sum _{l=1}^L\sum _{k=0}^{K}l\underset{l,k}{\overset{i,j}{u}}\left(t\right)x^k{y}^{l-1},$$

$${\left(u_{\mathit{xx}}\right)}^{i,j}\left(x,y,t\right)=\sum _{l=0}^L\sum _{k=2}^{K}k\left(k-1\right)\underset{l,k}{\overset{i,j}{u}}\left(t\right)x^{k-2}{y}^l,$$

$${\left(u_{\mathit{yy}}\right)}^{i,j}\left(x,y,t\right)=\sum _{l=2}^L\sum _{k=0}^{K}l\left(l-1\right)\underset{l,k}{\overset{i,j}{u}}\left(t\right)x^k{y}^{l-2}.$$E8

Based on the approach proposed by Lanczos [14, 15], the discretization of the problem is done by collocation at specially chosen gride points. For example, the grid points for the x variable over the interval [−1,1] in each subdomain are the K - 1 roots of a Chebyshev function, where K is the highest order of the power series used,

$$x_k=-\mathit{cos}\left\{\frac{\left(2k+1\right)\pi }{2\left(K-2\right)}\right\},k=0,1,\dots ,K-2,$$E9

and for the y variable,

$$y_l=-\mathit{cos}\left\{\frac{\left(2l+1\right)\pi }{2\left(L-2\right)}\right\},l=0,1,\dots ,L-2.$$E10

For each subdomain, the function u^i,j(x,y,t) as well as its derivatives are substituted into the governing differential equation at the grid points, (x_l, y_k), l = 0, 1,…, L-2 and k = 0, 1,…, K-2 to give (L - 1) x (K – 1) ODEs. On the four surfaces of each subdomain, boundary or interfacial continuity condition is specified on grid points: [x = ±1, y = ±1], [(x_l, l = 0,1 … L-2), (y = ±1)], and [(x = ±1), (y_k, k = 0, 1 … K-2)] to give a further (2 L + 2 k + 4) ODEs. The assemble of ODEs for the system is in the form,

$$\mathit{iA}{U}_t-L_{1}U-H_1\left(U,t\right)=0.$$E11

In Eq. (11), the unknown coefficients U is a [M x N x (K + 1) x (L + 1)]. A and L₁ are linear matrices, H₁ is a nonlinear vector. But their row dimensions are larger than the length of U. We use a discrete least square method to rectify this problem by multiplying Eq. (11) with A^T, the matrix transpose of A. The resultant matrix equations are well-posed to be solve by a linear equation solver with an iterative procedure to carter for the nonlinear term.

For a one-dimensional problem, N = 1 and L = 0, and in the m^th subdomain,

$${\mathbf{U}}^m={\left\{u_0^1,u_1^1,\dots ,u_k^1,u_0^2,u_1^2,\dots ,u_K^2,\dots ,u_0^M,u_1^M,\dots ,u_K^M\right\}}^{\prime }.$$E12

The matrix A in Eq. (11) consists of

$$\mathit{A}=\left[\begin{array}{llllll}{A}_1& & & & & \\ & A_2& & & & \\ & & ·& & & \\ & & & ·& & \\ & & & & ·& \\ & & & & & A_M\\ S_{1,1}& S_{1,2}& & & & \\ S_{2,1}& S_{2,2}& S_{2,3}& & & \\ & & ·& & & \\ & & & ·& & \\ & & & S_{M-1,M-2}& S_{M-1,M-1}& S_{M-1,M}\\ & & & & S_{M,M-1}& S_{M,M}\end{array}\right],$$E13

where $A_m=\left[\begin{array}{llllll}{x}_1^0& x_1^1& ·& ·& ·& x_1^M\\ x_2^0& x_2^1& ·& ·& ·& x_2^M\\ ·& ·& & & & ·\\ ·& ·& & & & ·\\ x_K^0& x_K^1& ·& ·& ·& x_K^M\end{array}\right],m=1,2,\dots M$,

and A_m is independent of m. There are M rows of S’s and each S is a 2 x (K + 1) matrix:

$$S_{1,1}=\left[\begin{array}{llllll}{\left(-1\right)}^0& {\left(-1\right)}^1& ·& ·& ·& {\left(-1\right)}^{K-1}{\left(-1\right)}^K\\ 0& 1& 2& ·& ·& \begin{array}{ll}\left(K-1\right)& K\end{array}\end{array}\right],$$

$$S_{1,2}=\left[\begin{array}{lllllll}0& & ·& & & ·& ·\\ 0& -1& -2{\left(-1\right)}^1& ·& ·& -\left(K-1\right){\left(-1\right)}^{K-2}& -K{\left(-1\right)}^{K-1}\end{array}\right],$$

and

$$S_{M,M-1}=\left[\begin{array}{lllllll}1& 1& ·& ·& ·& 1& 1\\ 0& & ·& & ·& & ·\end{array}\right]$$

$$S_{M,M}=\left[\begin{array}{llllll}{\left(-1\right)}^0& {\left(-1\right)}^1& ·& ·& ·& \begin{array}{ll}{\left(-1\right)}^{K-1}& {\left(-1\right)}^K\end{array}\\ 1& 1& ·& ·& ·& \begin{array}{lll}1& & 1\end{array}\end{array}\right].$$E14

Depending on the actual field equation involved, L₁ and H₁ can be constructed accordingly. Details of how the elements of each matrix could be determined are given in Ref. [8].

2.2 Solution by the real time evolution (RTE) method

This well-known method has been used to solve an initial-boundary problem [16] that evolves into a stationary and periodic solution. This method could be modified to cover cases where the solutions are periodical in time. In this special application, U is the same at the beginning and at the end of a period. To march from the beginning of a period to the end, we have chosen the implicit Crank-Nicholson stepwise formulation that is unconditionally stable. For Eq. (11),

$$\mathit{iA}\left(U^{m+1}-U^m\right)-\frac{\mathit{\Delta t}}{2}\left[L_1\left(U^{m+1}+U^m\right)+H_1\left(U^{m+1},t^{m+1}\right)+H_1\left(U^m,t^m\right)\right],$$E15

where the superscript m refers to the time step number. With the superscript r refers to the iteration number and the symbol ‘→’ means an integration, step, the iterative approach would be

$$U^{m+1,0}=U^m;\text{then}{U}^{m+1,r-1}\to U^{m+1,r}.$$E16

It should be noted that the iteration approach is needed due to the nonlinear H₁ terms. Since A and L₁ are both linear, only one inversion is required for all the iterative and time steps:

$$U^{m+1,r+1}={\left[\mathit{iA}-\frac{\mathit{\Delta t}}{2}{L}_1\right]}^{-1}\left(\left[\mathit{iA}+\frac{\mathit{\Delta t}}{2}{L}_1\right]U^m+\frac{\mathit{\Delta t}}{2}[H_1\left(U^{m+1,r},t^{m+1}\right)+H_1(U^m,t^m\left)\right]\right).$$E17

But starting with any pulse energy, a periodic solution may not exist. For this reason, we have developed a version of RTE method that, as shown later, the iteration will converge to an exactly periodic (EP) solution.

For a single mode problem, a term exp(iμ t₀) could be factored out from the solution of Eq. (17) at the position of the pulse peak. At a given time, t = t_o,

$$\mathrm{u}\left(\mathrm{x},\mathrm{y},{\mathrm{t}}_0\right)=exp\left({\mathrm{i\mu t}}_0\right)\widehat{\mathrm{u}}\left(\mathrm{x},\mathrm{y},{\mathrm{t}}_0\right).$$E18

Generally, if T is the period,

$$\mathrm{u}\left(x,y,t_0+T\right)=exp\left[\mathit{i\mu }\left(t_0+T\right)\right]\widehat{\mathrm{u}}\left(x,y,t_0+T\right),$$E19

and, for u(x,y,t₀) to be periodic,

$$\mathit{\mu T}=2\pi$$E20

and

$$\widehat{\mathrm{u}}\left(x,y,t_0\right)=\widehat{\mathrm{u}}\left(x,y,t_0+T\right).$$E21

As Eq. (17) can only be used to solve for exactly the number of coefficients in the series expansion, the pulse energy needs to be specified so that μ will be unique. For this reason, we have designed a set of iterative algorithms based on the pulse energy being of a specific value. To achieve this purpose, the pulse energy at the end of each iterative step is adjusted to the specific energy, eventually, the procedures lead to the correct μ and the converged pulse shape and pulse energy. The rate of convergence could be improved, if we also use the well-used averaging method [17] for u:

1. Start with

   $${\widehat{\mathrm{u}}}^0\left(x,y,0\right)=u\left(x,y,0\right),\to u^0\left(x,y,T\right).$$E22

2. Find μ_uT from u^m(x_o,y_o,T), then

$$\begin{array}{c}{\widehat{\mathrm{u}}}^m\left(x,y,T\right)=exp\left[-i{\mu }_{u}T\right]u^m\left(x,y,T\right),\\ w\left(x,y\right)=\frac{1}{2}\left[{\widehat{u}}^m\left(x,y,T\right)+u^m\left(x,y,0\right)\right],\\ u^{m+1}\left(x,y,0\right)=w^{m+1}\left(x,y\right)\sqrt{E/⟨|w|·|w|⟩}\to u^{m+1}\left(x,y,T\right),\end{array}$$E23

where the superscript m is the iteration number and E is the specified energy. The symbol ‘→’ indicates that, in each iteration, u^m(x,y,T) is obtained from u^m(x,y,0) using Eq. (17).

An obvious pre-condition is that the initial input used must be close to the converged solution. There are no set rules, but a Gaussian pulse is a good start. If error in the input has a negative imaginary component in its eigenvalue, the iteration will not converge due to modulation instability.

2.3 Numerical example of bimodal wave propagation

The governing equation for the spatiotemporal evolution of complex wave u(z, x, τ) and v(z, x, τ) in a planar waveguide is known [18] as

$$\begin{array}{c}{\mathrm{iu}}_z+\frac{1}{2}{\mathrm{u}}_{\mathit{xx}}+\frac{1}{2}{\mathrm{D}}_1{\mathrm{u}}_{\mathit{\tau \tau }}+\mathrm{vu}\ast =0,\\ 2\mathrm{i}\left({\mathrm{v}}_z+{\mathrm{cv}}_{\tau }\right)+\frac{1}{2}{\mathrm{v}}_{\mathit{xx}}+\frac{1}{2}{\mathrm{D}}_2{\mathrm{v}}_{\mathit{\tau \tau }}-\mathrm{qv}+\frac{1}{2}{\mathrm{u}}^2=0.\end{array}$$E24

where c is the group velocity mismatch parameter and q the phase mismatch constant. These equations are the same in form to those dealt with previously but with t and y replaced by z and τ respectively.

The system where v has twice the frequency of u is known as second harmonic generation. Traveling waves would split into a fundamental and a second harmonic modes. For such a system, the solution principles used in the RTE method remain the same with both u and v involved [19]:

$$\begin{array}{c}{\widehat{\mathrm{u}}}^m\left(x,,z\right)=exp\left(-\mathit{i\mu z}\right)u^m\left(x,,z\right),\\ {\widehat{v}}^m\left(x,\mathrm{\tau },z\right)=exp\left(-2\mathit{i\mu z}\right)v^m\left(x,\mathrm{\tau },z\right).\end{array}$$E25

As there are two pulse energies E_u and E_v now involved, we have three choices for assigning energy: The total energy E = E_u + E_v, or E_u and E_v by itself.

A system that supports the copropagating of two pulses of arbitrary frequencies may support also a continuously varying spectrum. For this reason, errors in the initial guess could grow with distance traveled. It is important to design an algorithm to ensure that the iterative procedures will lead only to the ground state solutions for u. It is noted that, at convergence, v will also assume equilibrium state. To implement these ideas into our algorithms, we set μ_u = 1 for u. At the end of each iterative cycle, we re-scale u so that μ could be forced to converge to 1. For v we do not preset μ_v and just let it assume its own value at convergence. The constraint we use for v is a specified energy ratio R = E_v/E_u. Again, E_v is not given a value at the beginning, but it will assume a value once a converged E_u is found. The algorithms are as follows:

1. Start with

   $$\begin{array}{c}{\widehat{\mathrm{u}}}^0\left(x,y,0\right)=u\left(x,y,0\right),\to u^0\left(x,y,T\right),\\ {\widehat{v}}^0\left(x,y,0\right)=v\left(x,y,0\right),\to v^0\left(x,y,T\right).\end{array}$$E26

2. Find μ_uT from u^m(x_o, y_o, T) and μ_vT from v^m(x_o, y_o, T), then

$$\begin{array}{lll}{\widehat{\mathrm{u}}}^m\left(x,y,T\right)& =& exp\left[-i{\mu }_{u}T\right]u^m\left(x,y,T\right),\\ {\widehat{v}}^m\left(x,y,T\right)& =& exp\left[-i{\mu }_{v}T\right]v^m\left(x,y,T\right),\\ \mathit{uu}\left(x,y\right)& =& \frac{1}{2}\left[{\widehat{u}}^m\left(x,y,T\right)+u^m\left(x,y,0\right)\right],\\ \mathit{vv}\left(x,y\right)& =& \frac{1}{2}\left[{\widehat{v}}^m\left(x,y,T\right)+v^m\left(x,y,0\right)\right],\\ r_u& =& \sqrt{\frac{1}{{\mu }_u}},\\ u^{m+1}\left(x,y,0\right)& =& r_u\mathit{uu}\left(x,y\right)\to u^{m+1}\left(x,y,T\right),\\ \\ r_v& =& \sqrt{R\frac{⟨|u^{m+1}\left(x,y,0\right)|·|u^{m+1}(x,y,0)|⟩}{⟨|\mathit{vv}\left(x,y\right)|·|\mathit{vv}(x,y)|⟩}},\\ v^{m+1}\left(x,y,0\right)& =& r_v\mathit{vv}\left(x,y\right)\to v^{m+1}\left(x,y,T\right).\end{array}$$E27

We have applied the LCPS method to Eq. (24) and obtain a set of ODEs that was solved with RTE method to give a set of stationary solutions. The complicated waveforms found could be seen from Figure 1. We then propagate this set of solutions over a distance z = 10. The propagation histories show that there is no change in the pulses amplitude and energy as can be seen from the plots in Figure 2.

Figure 1.

Stationary solutions found for Eq. (24) with E = 400, D2 = − 0.2 and q = 2. (symmetry in the central x-plane was used with 4 x 2 subdomains and K = L = 8. The initial guesses for u and v were Gaussian pulses in both directions).

Figure 2.

Stable propagation of the stationary solutions shown in Figure 1.

3. Electromagnetic wave propagation through space

Numerical procedures described in the previous sections have been modified and used to study the propagation characteristics of electromagnetic waves in the form of bright, dark, and anti-dark solitons [10, 11, 20]. The steps needed are to be described below.

3.1 Stable periodic (SP) soliton solutions of NLSE

For a plane wave the governing NLSE and boundary conditions are

$$u_x-\frac{i}{2}D\left(x\right)u_{\mathit{tt}}-\mathit{i\gamma }{\left|u\right|}^{2}u=0,u\left(0,x\right)=u\left(L,x\right)=0,$$E28

where u is the slow varying envelope of the axial electric field, $D\left(x\right)$and $\gamma$ represent the dispersion coefficient and self-phase modulation parameters, respectively. x and t are the spatial propagation distance and temporal local time, respectively. L is the width of the numerical window used for t. For application to space, x is a very large number while D and $\gamma$ are very small. To eliminate possible numerical complications associated with those numbers, scaling factors, x_o and t_o, are introduced so that

$$x^{\ast }=\frac{x}{x_o},t^{\ast }=\frac{t}{t_o},$$E29

then, together with

$$D^{\ast }=\frac{D{x}_o^{0.5}}{{\gamma }^{0.5}},u^{\ast }={\left(\gamma x_o\right)}^{0.5}u,$$E30

Eq. (28) becomes dimensionless,

$$u_x-\frac{i}{2}D\left(x\right)u_{\mathit{tt}}-i{\left|u\right|}^{2}u=0,$$E31

where the superscript * has been omitted for simplicity.

To solve Eq. (31) numerically, consider pulse propagation as a transient problem along the spatial distance, x, the discretization is one dimensional and only at the temporal local time domain t. Using M subdivisions

$$\Omega ≔{\Omega }^i,i=1,2,\dots ,M.$$E32

For each subdivision, the numerical window of length L is mapped into an interval varying from – 1 to +1, and an economized power series is used:

$$u\left(t,x\right)=\sum _{k=0}^K{u}_k\left(x\right)t^k.$$E33

Applying the LCPS method to Eq. (31) at the collocation points, t₀, t_1… t_K, to each subdomain leads to a set of ODEs. The assembly of all subdomains involves the series expansion coefficient as a vector of length $\left(K+1\right)\times M$:

$$\mathit{u}={\left\{u_0^1,u_1^1,\dots ,u_K^1,u_0^2,u_1^2,\dots ,u_K^2,\dots ,u_0^M,u_1^M,\dots ,u_K^M\right\}}^{\prime }.$$E34

Between two adjacent subdomains, i and i + 1, the continuity conditions are:

$$u^i{\left|{}_1=u^{i+1}\right|}_{-1};\frac{d}{\mathit{dx}}{u}^i{\left|{}_1=\frac{d}{\mathit{dx}}{u}^{i+1}\right|}_{-1}.$$E35

The set of transient ODEs obtained is in the form,

$$\mathit{A}{\mathit{u}}_x\left(x\right)-i\mathit{Lu}\left(x\right)=i\mathit{Q}\left(x,\mathit{u}\right).$$E36

Applying the RTE method described in Section 2.2 to the above equation,

$$\mathit{A}\left({\mathit{u}}^{m+1}+{\mathit{u}}^m\right)-\frac{i∆x}{2}\left[\mathit{L}\left({\mathit{u}}^{m+1}+{\mathit{u}}^m\right)\right]=\frac{i∆x}{2}\left[\mathit{Q}\left(x,{\mathit{u}}^{m+1}\right)+\mathit{Q}\left(x,{\mathit{u}}^m\right)\right],$$E37

where ∆x is the step size, and the superscript m refers to the time step number. To carter for the nonlinear nature of Eq. (37) that is associating with Q, an iterative algorithm [9, 10] is used.

The initial input pulse for a bright soliton could be

$$u\left(t,0\right)=\beta exp[-\alpha \left(t-0.5L){}^2\right],$$E38

where L is the numerical window used for t, $\alpha$ a chosen constant to give an input pulse as close to the SP soliton as possible, and $\beta$ an adjusting parameter to give a specified pulse energy, E,

$$E\left(x\right)=\underset{-\frac{L}{2}}{\overset{\frac{L}{2}}{\int }}(\left|{\left.u\left(t,x\right)\right|}^2\right)\mathit{dt}.$$E39

As u(t,x) is a truncated soliton pulse, it has been found [9] that to eliminate residual reflection, the following boundary conditions could be used:

$$u\left(t,x\right)=1000\frac{\partial u}{\partial t}-u\left(t,x\right)\mathit{at}x=\pm 0.5L.$$E40

To find the stable periodical (SP) solution, Eq. (37) is integrated to a selected distance Z, with the first half using a specified dispersion coefficient of – D, and for the second half D. For an SP solution, the input pulse must be the same as the output pulse. We use this fact to design an iterative scheme based on successive halves,

$${\mathit{u}}_{\mathit{\text{in}}}^{m+1}=0.5\left({\mathit{u}}_{\mathit{\text{in}}}^m+{\mathit{u}}_{\mathit{\text{out}}}^m\right),$$E41

where u_in, and u_out are the input and output pulse to the dispersion map respectively and the superscript m denote the iteration number. It should be noted that stable periodic solitons are special cases of cyclic solitons in that no phase matching is needed. For the exact periodic (SP) solutions described in Section 2.2. the input and output pulses have the same amplitude and phase.

3.2 A numerical example of a SP bright soliton

Using the procedure described in the previous section and a dispersion map with length Z = 6, Figure 3 shows how the solutions converged to stable and periodic pulses [10]. The distance, x, shown is the cumulated distance. As the step size is 0.0005, each iteration generates 12,000 pulses. In the last few iterative cycles, the pulse width is changing linearly with the distance traveled in both halves of the dispersion map. The input and output pulses in a dispersion map have the same shape and energy but not the same phase.

Figure 3.

Iteration convergence for the numerical example.

Figure 4 shows changes of the pulse shape as an SP soliton is traveling through a medium with negative D.

Figure 4.

Pulse changes when propagating through a medium with -D.

4. Electromagnetic wave propagation in a complex space system

Some of the transmission characteristics of electromagnetic waves through space have been investigated previously [11, 20]. We shall deal with a space that has other complex features in the following sections.

4.1 SP solitons with different pulse energies and in a space with random dispersion coefficient

We have solved for a segment consisting of piecewise continuous dispersion coefficients as given in Table 1 below. The pulse width histories for cases with pulse energy E = 0.4 and 0.8 found for the above cases are shown in Figure 5. The plots show that the overall pulse width change for the random dispersion case is the same as that based on the averaged dispersion. Also in the plots are equations of the trendlines showing the linear relationship between pulse width and distance traveled when the average coefficient D is used. Also, for two times increase pulse energy, the deviation from the average gradient is only +/− 4%.

	Case 1			Case 2
Section	Propagation distance, x	Dispersion coefficient, D	External source, uo	Propagation distance, x	Dispersion coefficient, D	External source, uo
1	2	−0.1	0	2	−0.2	0
2	0.5	−0.1	0	0.5	0.2	−0.2
3	1	−0.1	0	1	0.1	0
4	0.5	−0.1	0	0.5	−0.2	0.2
5	2	−0.1	0	2	−0.15	0

Table 1.

Two cases of propagation through different D and u_o.

Figure 5.

Propagation of SP soliton with different pulse energy.

4.2 A system of multiple SP solitons

The NLQE for multiple solitons is

$$u_x^j-\frac{i}{2}D\left(x\right)u_{\mathit{tt}}^j-\mathit{i\gamma }\left\{\sum _{k=1}^J{\left|u^k\right|}^2\right\}{u}^j=0,\mathrm{j}=1,2,..\mathrm{J}.$$E42

Although a set of equations, equal in number to the number of solitons, are involved, the same numerical procedures described previously can be used to obtain a set of SP solitons. For our purpose, however, we only use three well-spaced and identical pulses. Figure 6 shows how the central pulse has converged to an SP soliton. The same applies to the other two solitons.

Figure 6.

Iteration convergence for the central pulse.

How the pulse shape is changing can be seen in Figure 7. The gradient of pulse width changes against distance traveled is not sensitive to pulse energy as can be seen in Figure 8.

Figure 7.

Changing shape as pulse propagating through a region where D is negative.

Figure 8.

Half pulse width changes versus distance traveled at different pulse energy.

4.3 Propagating through space with a CW background

To include a constant CW background, u_o, into NLSE, let

$$u=v+u_o.$$

Substituting the above into Eq. (31) to give

$$v_x-\frac{i}{2}D\left(x\right)v_{\mathit{tt}}-i{\left|v+u_o\right|}^2\left(v+u_o\right)=0.$$E43

Using the same numerical procedures as described previously, Eq. (43) could be solved to give an SP solution. By propagating this SP soliton along the distance x, the transmission characteristics could be determined. We have solved for two cases with different system parameters as shown in Table 1. In fact, Case 1 is using a constant D that is the same as the distance weighted average of D in Case 2. The pulse width and energy histories are plotted out in Figure 9. It could be seen that u_o has little influence on the overall pulse width change.

Figure 9.

Two cases of propagation through space with CW background.

4.4 Propagation through an amplifying or attenuating space

The NLS equation for electromagnetic waves (solitons) propagation in dimensionless form and in an attenuating space is

$$u_x-\frac{i}{2}D\left(x\right)u_{\mathit{tt}}-i{\left|u\right|}^{2}u=S\left(x\right),$$E44

where $S\left(x\right)=\mathit{su}$ for amplification, and $S\left(x\right)=-\mathit{su}\left(x\right)$ for attenuation and s is a constant.

To solve the above equation numerical, S must be added to Q in Eq. (36), As a numerical example, a system consists of three sections was used. In all sections, D = −0.2 but s = 0.5, − 0.5 and 0.5, respectively. Solutions are shown in Figure 10. Features of the solution histories fond are: (a) based on the sign of s, the pulse energy increases or decreases steadily, and (b) practically, there is no change in the gradient of the wavelength half width versus x curve.

Figure 10.

Propagation histories in the present of an external source.

4.5 Propagation with gravitational deflection

Based on the general relativity theory, the approximate light path deflection angle, ∆ϴ is found to be [21],

$$∆ϴ=\frac{4\mathit{GM}}{∆},$$E45

where d the distance between the light path and the center of the mass, G is the gravitation constant, M is the mass, and ∆ is the distance between the wave front and the center of the mass.

Eq. (45) could be used in its dimensionless form,

$$∆ϴ=\frac{C}{\mathcal{E}},$$E46

where $\mathcal{E}=\frac{C∆}{4\mathit{GM}}$. Since the event is taking place in space, we have no way of knowing M and ∆. But we can still track the gravitational deflection history using a single arbitrarily chosen parameter C. Using a new rectangular coordinate system (x1, x2), at a particular step, let ((x1)_m, (x2)_m) be the position of the mass center and ((x1)₁, (x2)₁) the wave front; the straight line connecting the wavefront to the center of the mass is

$$\mathcal{E}=\sqrt{{\left({\left(x1\right)}_m-{\left(x1\right)}_1\right)}^2+{\left({\left(x2\right)}_m-{\left(x2\right)}_1\right)}^2}.$$E47

Then, by specifying a C, the deflection angle $∆ϴ$ can be found from Eq. (46). If a wavefront is propagating along a light path making an angle ϴ with the x1-axis. Integrating along x1 with step ∆x, the new wave front position would be [((x1)₁ + ∆x) cos (ϴ + $∆ϴ\left)\right)$, ((x2)₁ + ∆x) sin (ϴ + $∆ϴ$)]. Knowing the new position, Eq. (46) and Eq. (47) could be used to find $\mathcal{E}$ and $∆ϴ,$and the next wavefront position in the next integration step.

To implement this deflection scheme, let Z be the length of the propagation distance to be investigated. Let t(0, 0) be the starting position with the mass M located along the straight line x1 = 0.5 Z. If φ is the angel between the line joining the wavefront to the center of the mass and the x1-axis, the mass is located at (0.5 Z, 0.5 Z tan(φ)). For this example, let the initial φ = 20, ϴ = 0, and Z = 3, and D = − 0.2. Considering two cases, Case 3, C = 0.00005, and Case 4, C = 0.0001, respectively. After solving for the wavefront histories, the solutions are plotted out in Figure 11. It can be seen that, at x1 = 0.5 Z where the wavefront is closest to the mass, the deflection rate is the largest. As expected, a larger C will give a larger deflection. Without deflection, the wavefront will move along the x1-axis. With deflection, the wavefront has traveled the same distance along its path, but shorter in term of x1-coodinate. It should be pointed out x1 and x2 are scaled down to dimensionless quantities chosen according to local conditions. They would be many orders smaller than the entire propagation distance.

Figure 11.

Propagation in the present of gravitational deflection.

4.6 Calibrations with physical systems

The coefficients associated with dispersion and self-phase modulation for space cannot be measured directly. To apply our numerical findings to cosmological redshift, calibration with measured data must be used as described in Ref. [10, 19]. If the starting and ending wavelength is λ₁ and λ₂, and the half pulse width, W₁ and W₂, and the dimensionless redshift, z, as defined in astronomy, is given by,

$$z=\frac{{\lambda }_2-{\lambda }_1}{{\lambda }_1}$$

$$=\frac{W_2-W_1}{W_1}.$$E48

Also, in astronomy, the measured redshift-distance relation is given by Hubble’s law

$$z=H_{o}d/c,$$E49

where H_o is the Hubble constant given in units of km/s/Mpc. The distance between the light source and an observer, d is in Mpc unit, and c is the speed of light in km/s. Since Eq. (49) is linear with specific physical dimensions, and our numerical results are also linear but dimensionless, we could choose any two given points in Eq. (49) for calibration, so that our results agree with Hubble’s law.

To use real physical data for calibration, we choose, as an example, H_o = 70 km/sec/Mpc for Eq. (49). We also choose the pulse representing the spectral line due to Lyman-alpha hydrogen that has a wavelength of 121.6 nm., and the corresponding period is 405 ps. For this example, we choose the linear relationship that we have found previously for a random medium and shown in Figure 5,

$$W=0.9539x+2.5834,$$E50

and the half pulse width at the two calibration points, A at x_A = 2, and B at x_B = 6, from Eq. (50), are W_A = 4.9286, and W_B = 8.3068. Then, the temporal time multiplying scaling factor to convert W into period in picosecond,

$$t_o=\frac{405}{W_A}=82.1734\mathrm{ps}.$$E51

Now, the redshift between A and B,

$$z_{\mathit{AB}}=\frac{W_B-W_A}{W_A}=0.6854.$$E52

Using Eq. (49), the distance variable between the source and the observer,

$$d/c=\frac{z_{\mathit{AB}}}{H_o}=\frac{0.6854}{70}0009792\text{Mpcs}/\mathrm{km}$$E53

$$d=0.009792c\times \frac{3.08567758\times {10}^{19}}{9.461\times {10}^{12}}=9.581{10}^9\mathit{\text{light year}}$$E54

Using Eq. (53), the distance scaling factor used to convert x to d/c is

$$f_x=\frac{d/c}{x_{A-x_B}}0002448\text{Mpcs}/\mathrm{km}$$E55

5. Discussion

We have demonstrated that numerical methods can be used to solve many nonlinear field problems. The general mathematical procedures involve the reduction of the governing differential equations to a set of nonlinear algebraic equations that is solved by an iterative approach. Depending on the expected characteristics of the solutions, numerical algorithms may need to be modified. For an exactly stationary and periodic traveling wave, the iterative procedures as described in Section 2.2 are needed. But for the studies of redshift, only a stable and periodic solution in the sense, that the traveling pulse can reproduce itself anywhere in the whole propagation history, the iterative procedures needed are less elaborative as described in Section 3.1. Since the system is a nonlocal problem, the finding of the correct initial pulse is vital for the investigation of redshift in starlight.

The soliton solutions’ spike-like characteristic needs a very high order of series approximation. We choose to sub-divide the computational domain into zones so that a lower order series could be used for each zone. In theory, the tail of a soliton extends to zero value at infinity. Because we can only use a finite computational window, we are solving for a truncated soliton. Therefore, the use of boundary condition, Eq. (40), is important to limit reflection at the boundary. As the pulse width is changing with distance traveled, the choices for an initial pulse width and the size of the computational window need careful design. Providing a numerical method can satisfactorily reproduce a stable and periodic soliton, there is no reason why such a method cannot also be used.

For mathematical simulation to be an effective tool employed to design, operate, or control a system, a set of well-proven relationships between the dependent and independent variables must be used, together with their system parameters. If not all the parametric data are known, calibration can be carried out to find the missing ones. However, it is important that the simulation must be used within a valid range of all parameters. What has been done in Section 4.6 is special because the linear relationship found in the numerical simulation is the same as the empirical Hubble law.

In Section 4, we have used numerical examples to show that electromagnetic wave can travel through sections of space that have different dispersion coefficient, but the overall wavelength change is equal to that based on the averaged coefficient (See Section 4.1). This characteristic is because soliton behave both as a wave and a particle. Similarly, we have shown that the present of other solitons (See Section 4.2), the present of a CW background (See Section 4.3), or the present of a source, do not significantly change the rate of wavelength changes against distance traveled. The last case, Section 4, is about the effect of gravity; the estimate found is approximate. But, since the space is known to be a virtual empty void even with the present of more than 200 billion galaxies, a journey through the universe would experience gravitational effect only over a miniscule portion of the total length. We can conclude that gravity will contribute to redshift but not significantly.

Hubble tension is a typical case where prediction from the standard model of cosmology is at variant to the experimentally observed empirical Hubble law [12, 13]. As some of the parameters used in the model were based on small redshift data available on the time, it not surprising to see that the model predicts differently when data for larger redshift are available. Based on the cosmological principle, the model uses a homogeneous and isotropic universe. But the real space has large clusters and voids even measured in the cosmological scale. Redshift in starlight traveling through a cluster will be quite different to the one traveling through a void.

The propagation theory used in our investigation is built on NLSE that ensures a balance between linear dispersion and nonlinear self-phase modulation. The space is never assumed as isotropic or homogeneous. By solving the NLSE, the traveling histories are constructed according to local conditions encountered. In this approach, we can accommodate the complex nature of the universe.

6. Conclusion

1. The LCPS method is an effective and economical method to convert a differential equation into a set of nonlinear algebraic equations that can then be solved numerically.

2. The propagation theory based on the NLSE predicts redshift in electromagnetic wave through space can be calibrated to agree with Hubble law.

3. The complex nature of the universe does alter the particle nature of solitons, allowing changes during their transmission histories to be accumulated.

4. In mathematical simulations, the limits on their ranges of applicability should be recognized.

Conflict of interest

The author declares no conflict of interest.