High-precision numerical analysis method and system for wave propagation in nonlinear media

Abstract

The present invention discloses a high-precision numerical analysis method and system for wave propagation in nonlinear media, which relates to the field of wave propagation application technology. The method includes: constructing a model equation for wave propagation in a nonlinear medium; discretizing the model equation using a finite difference method to obtain a high-precision difference format that approximates the solution of the model equation; designing a numerical algorithm for solving the model equation based on the high-precision difference format; and numerically simulating and analyzing the propagation of waves in nonlinear media. The present invention uses a computer method to simulate and analyze the propagation phenomenon of waves in nonlinear media. Because the numerical method proposed by the present invention has higher precision in both time and space directions, it can more accurately and effectively simulate the wave propagation process. The time extrapolation algorithm and the fast solution algorithm are effectively utilized, which will improve the overall computational efficiency of the numerical algorithm and reduce the cost of computer numerical simulation.

Description

Technical Field

[0001] The present invention relates to the technical field of wave propagation applications, and in particular to a high-precision numerical analysis method and system for wave propagation in nonlinear media. Background Art

[0002] Since the 1950s, the phenomenon of single-wave propagation or nonlinear interactions between multiple waves has been widely observed in various fields of physics and applied physics, particularly fluid dynamics, solid-state physics, nonlinear optics, and plasma physics. For example, when multiple waves meet specific resonance conditions in terms of wavenumber and frequency, they form resonant interactions, a phenomenon widely found in engineering applications in the field of water waves. In addition, several common wave propagation phenomena exist in nonlinear optics, including light beam transmission and the propagation of light in photonic crystals. These phenomena primarily reflect the influence of nonlinear effects on light propagation in nonlinear media. For example, in optical fibers, the propagation of light beams is affected by nonlinear effects, resulting in self-focusing of the beam and the generation of optical solitons.

[0003] To better understand wave propagation, effective methods are needed to describe its behavior in nonlinear media. Currently, there are two main approaches to studying wave propagation: experimental and numerical simulation. While experimental methods can reveal the actual propagation behavior of waves under specific conditions, any change in conditions necessitates repeated experiments, resulting in high costs and long processing times. In contrast, numerical simulations, by establishing equations modeling wave propagation in nonlinear media and applying corresponding numerical methods to simulate different propagation scenarios, are more cost-effective and time-efficient. Currently, a variety of numerical methods have been proposed for wave propagation in nonlinear media, such as the finite difference method, the finite volume method, the finite element method, and the differential integral method. However, most existing numerical calculation methods suffer from low precision in both time and space, resulting in inaccurate results, large errors, and low resolution when using coarse grids. To obtain more accurate and high-resolution results, finer grids are required, but this significantly increases the computational workload and reduces computational efficiency. Summary of the Invention

[0004] The present invention aims to propose an equation model describing the temporal evolution of the wave function, based on the physical properties of wave propagation in nonlinear media. Subsequently, a corresponding high-precision difference scheme is established using the finite difference method in numerical computation. To improve the efficiency of numerical simulations, the present invention designs a time extrapolation algorithm and a fast solution algorithm compatible with the high-precision difference scheme, enabling accurate and efficient numerical simulations of wave propagation in nonlinear media.

[0005] According to a first aspect of an embodiment of the present disclosure, a high-precision numerical analysis method for wave propagation in a nonlinear medium is provided, comprising the following steps:

[0006] Construct model equations for wave propagation in nonlinear media;

[0007] The model equation is discretized using the finite difference method to obtain a high-precision difference format for the approximate solution of the model equation;

[0008] Designing a numerical algorithm for solving the model equations based on a high-precision difference scheme;

[0009] Numerical simulation and analysis of wave propagation in nonlinear media.

[0010] According to a second aspect of an embodiment of the present disclosure, a high-precision numerical analysis system for wave propagation in a nonlinear medium is provided, comprising:

[0011] Model acquisition module, which constructs the model equations of wave propagation in nonlinear media;

[0012] The format design module uses the finite difference method to discretize the model equation and obtain a high-precision difference format that approximates the solution of the model equation;

[0013] A solution algorithm design module, which designs a numerical algorithm for solving the model equation based on a high-precision difference format;

[0014] The analysis module performs numerical simulation and analysis on wave propagation in nonlinear media.

[0015] According to a third aspect of an embodiment of the present disclosure, an electronic device is provided, comprising a memory, a processor, and a computer program stored and running on the memory, wherein when the processor executes the program, the high-precision numerical analysis method for wave propagation in nonlinear media is implemented.

[0016] According to a fourth aspect of an embodiment of the present disclosure, a computer-readable storage medium is provided, on which a computer program is stored. When the program is executed by a processor, the high-precision numerical analysis method of wave propagation in nonlinear media is implemented.

[0017] The above technical solution adopted by the present invention has the following advantages compared with the prior art:

[0018] 1. The present invention uses computer methods to simulate and analyze wave propagation in nonlinear media. Because the numerical method proposed in the present invention has higher precision in both time and space, it can more accurately and effectively simulate the wave propagation process.

[0019] 2. The present invention effectively utilizes the time extrapolation algorithm and the fast solution algorithm, which will improve the overall computational efficiency of the numerical algorithm and reduce the cost of computer numerical simulation.

[0020] 3. This invention can deepen people's understanding of wave propagation phenomena and provide strong theoretical support, verifiable high-precision numerical calculation and analysis tools for practical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] The drawings in the specification, which constitute a part of this application, are used to provide further understanding of this application. The illustrative embodiments of this application and their descriptions are used to explain this application and do not constitute improper limitations on this application.

[0022] Figure 1 This is a comparison chart of the computation time of the two algorithms under different grids;

[0023] Figure 2 Surface diagram of soliton formation and propagation in nonlinear media;

[0024] Figure 3 Projection diagram of the formation and propagation of solitons in nonlinear media;

[0025] Figure 4 This is the projection diagram of two waves propagating in opposite directions in a nonlinear medium. Specific implementation methods

[0026] The present disclosure will be further described below with reference to the accompanying drawings and embodiments.

[0027] It should be noted that the following detailed descriptions are illustrative and intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present application belongs.

[0028] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or combinations thereof.

[0029] It should be noted that the flowcharts and block diagrams in the accompanying drawings illustrate the possible implementation architecture, functions and operations of the methods and systems according to the various embodiments of the present disclosure. It should be noted that each box in the flowchart or block diagram can represent a module, program segment, or a part of code, and the module, program segment, or a part of code can include one or more executable instructions for implementing the logical functions specified in the various embodiments. It should also be noted that in some alternative implementations, the functions marked in the box can also occur in an order different from that marked in the accompanying drawings. For example, two boxes represented in succession can actually be executed substantially in parallel, or they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the flowchart and / or block diagram, and the combination of the boxes in the flowchart and / or block diagram, can be implemented using a dedicated hardware-based system that performs the specified function or operation, or can be implemented using a combination of dedicated hardware and computer instructions.

[0030] Example 1:

[0031] This embodiment provides a high-precision numerical analysis method for wave propagation in nonlinear media, including the following steps:

[0032] Step 1. Construct the model equations for wave propagation in nonlinear media;

[0033] Specifically, assume that the wave propagates in the two-dimensional region Ω = [a1, b1] × [a2, b2], where the wave amplitude is the function φ(x, y, t), the potential energy in the nonlinear medium is represented by f(x, y), and the nonlinear coefficient during propagation is β. In this way, the model equation satisfied by the wave amplitude during propagation in the nonlinear medium is

[0034]

[0035] Where Δ is the Laplace operator, α is a constant, i is an imaginary unit, x and y are spatial direction variables, t is a time variable, and T is the total time required for the calculation to arrive. The displacement and velocity of the wave propagation at the initial moment are

[0036]

[0037] where ψ(x,y) and is a smooth function.

[0038] In order to simulate the propagation of waves in nonlinear media, it is necessary to discretize the model equation (1) to establish a high-precision difference format. Based on the above goal, a new variable function is introduced Formula (1)-(2) can be rewritten as

[0039]

[0040] Step 2. Use the finite difference method to discretize the model equation and obtain a high-precision difference format for the approximate solution of the model equation;

[0041] Specifically, the time step is defined as The calculation time is evenly divided into K intervals. At the same time, the calculation domain is evenly meshed. Given two positive integers M x and M y , define the x-direction spatial step size as and the spatial step length in the y direction is That is, [a1, b1] × [a2, b2] is evenly divided into M x ×M y intervals. Each calculation grid point is recorded as (x j ,y l ,t k ), where x j =a1+jh x ,y l =a2+lh y ,t k =kτ. In addition, the grid function is defined

[0042] Define the difference operator

[0043]

[0044] remember and is the point (x j ,y l ,t k ), using the Crank-Nicolson method for time and the sixth-order difference operator for space, we can obtain the following unconditionally stable high-precision difference format:

[0045]

[0046] The periodic boundary condition is expressed as

[0047] The above high-precision differential format can be expressed in matrix form

[0048]

[0049]

[0050] where Φ, P, and f are Φ j,l ,P j,l ,f j,l The matrix formed by all discrete points in the calculation area, and

[0051]

[0052] Here circ indicates that the matrix is ​​a circulant matrix. Represents the matrix A x The inverse matrix of A y The inverse matrix of this high-order format is

[0053] Step 3. Design a numerical algorithm to solve the model equation based on a high-precision difference format;

[0054] Specifically, the numerical algorithm for solving the model equations includes a time extrapolation algorithm. The accuracy of equations (9)-(11) in the time direction is only second order. In order to achieve the theoretical sixth order accuracy in the spatial direction during calculation, τ must be taken to be sufficiently small, but this is very time-consuming. Therefore, the established high-precision difference format is upgraded to sixth order through the following Richardson extrapolation technique:

[0055]

[0056] The numerical algorithm for solving the model equations includes a fast solution algorithm. In order to improve the computational efficiency in the subsequent numerical calculation, the following corresponding fast solution algorithm is designed, which equates the matrix form (9)-(11) of the high-precision difference format to:

[0057]

[0058] in

[0059]

[0060]

[0061] Circulant matrix B x , B y , G x and G y Can be diagonalized

[0062]

[0063] in They are the circulant matrices B x ,B y ,G x ,G y The diagonal matrix formed by the eigenvalues, T x ,T y is the reversible matrix composed of all eigenvectors of the circulant matrix. Further sorting out equations (13)-(14), we get:

[0064]

[0065] According to the diagonalization of the circulant matrix, Equation (17) can be rewritten as

[0066]

[0067] Multiply both sides of equation (18) by T x , right multiplication set up get

[0068]

[0069] in and 1 is a matrix whose elements are all 1. In the specific numerical calculation, the following iterative format is used

[0070]

[0071] Where . / and .* represent element division and multiplication of the matrix, when ||Φ k+1(s+1) -Φ k+1(s) || ∞ ≤10 -12 When Φ k+1 =Φ k+1(s+1) , we get Φ k+1 Then, calculate P by formula (16) k+1 .

[0072] The numerical algorithm for solving the model equations includes a time-advancing algorithm, which is based on a time extrapolation algorithm and a fast solution algorithm, combined with a high-precision difference format. The entire time-advancing process for solving equations (1)-(2) is given as follows:

[0073] Step 1: Take the time step τ=λ·max(h x ,h y ),λ∈R + , get the initial value and

[0074] Step 2: Iterate formula (20) to get the function value of the first time layer

[0075] Step 3: Obtain the function value of the first time layer through formula (16)

[0076] Step 4: Repeat Step 2-3 until the final moment, and get and

[0077] Step 5: Re-take the time step τ=λ·max(hx ,h y ) / 2, repeat Step 2-Step 4 to get and

[0078] Step 6: Re-take the time step τ = λ·max(h x ,h y ) / 4, repeat Step 2-Step 4 to get and

[0079] Step 7: Finally, we can get the new and

[0080] Step 4. Perform numerical simulation and analysis on wave propagation in nonlinear media.

[0081] To verify the theoretical results of the proposed format, three nonlinear problems were solved using the format. For simplicity, equations (13)-(15) proposed in this invention are denoted as SP-HOD, and equations (15)-(16) and (20) using the fast solution algorithm are denoted as SP-HOD-Fast. Unless otherwise specified, h = h is used to solve the following problems. x =h y .

[0082] First, we verify the high-precision numerical format through a numerical example 1 with an exact solution. We select τ = 0.01, t = 2 to compare the computational efficiency of the algorithms SP-HOD and SP-HOD-Fast. Figure 1 It can be seen that SP-HOD-Fast, which uses a fast algorithm, requires less calculation time under the same number of grids, and the calculation efficiency becomes higher and higher as the number of grids increases. This also reflects the superiority of the fast algorithm designed by the present invention. At the same time, Table 1 gives the numerical errors of the two algorithms under different spatial step sizes. It can be seen that the calculation results of the two algorithms are consistent and there is no loss of numerical accuracy. This is consistent with the theoretical derivation and analysis. In addition, the results in Table 2 verify that the numerical algorithm proposed in the present invention is still applicable when the step sizes in the two spatial directions are not equal. Table 3 verifies the stability of the format through numerical means, showing that the format is unconditionally stable. This means that a larger time step can be taken during calculation, which will undoubtedly improve the computational efficiency of the problem. Example 2 depicts the phenomenon that a soliton is formed in the center of the nonlinear medium region and diffuses outward, and the step size h = 1 / 10 and τ = 1 / 10 are used for simulation. From Figure 2 and Figure 3It can be seen that the soliton in the center of the region gradually spreads to the surrounding areas as time evolves, while maintaining the soliton's characteristics and propagation speed. This result is consistent with the physical problem described. Example 3 simulates the interaction between two waves propagating in opposite directions in a nonlinear medium, taking the step size h = 1 / 10 and τ = 1 / 10 for the simulation. Figure 4 It can be seen that as time evolves, a wave in the region splits into two waves that propagate in opposite directions. This is consistent with the wave propagation law characterized by the physical problem described by the equation model.

[0083] Example 1 considers the nonlinear medium region [-1,1] 2 The propagation of internal waves, taking the parameter α=β=1, the potential energy of nonlinear medium Initial displacement Initial velocity The exact solution to this problem is

[0084] Table 1 L obtained by two algorithms for different h when T=1 in Example 1 ∞ and L2 error and convergence order (Rate)

[0085]

[0086] Table 2 L calculated in Example 1 using the SP-HOD-Fast format when τ = 0.01, T = 1 ∞ and L2 error and convergence order (Rate)

[0087]

[0088] Table 3 L calculated in Example 1 using the SP-HOD-Fast format when h = 1 / 50 and T = 5 ∞ and L2 error

[0089]

[0090]

[0091] Example 2: Consider the nonlinear medium region Ω=[-20,20] 2 A soliton is formed at the center and propagates outward. The corresponding parameters of this problem are α=β=1, the potential energy of the nonlinear medium is zero, and the initial displacement is Initial velocity ψ(x,y)=sin(x+y)exp(-2(x 2 +y 2 )).

[0092] Example 3: Consider the nonlinear medium region Ω=[-40,40] 2The two waves propagating in opposite directions interact with each other. The corresponding parameters of this problem are α=β=1, the potential energy of the nonlinear medium is zero, the initial velocity ψ(x,y)=0, and the initial displacement

[0093] Example 2:

[0094] This embodiment provides a high-precision numerical analysis system for wave propagation in nonlinear media, including:

[0095] Model acquisition module, which constructs the model equations of wave propagation in nonlinear media;

[0096] The format design module uses the finite difference method to discretize the model equation and obtain a high-precision difference format that approximates the solution of the model equation;

[0097] A solution algorithm design module, which designs a numerical algorithm for solving the model equation based on a high-precision difference format;

[0098] The analysis module performs numerical simulation and analysis on wave propagation in nonlinear media.

[0099] Example 3:

[0100] An electronic device includes a memory, a processor, and a computer program stored and executed on the memory, wherein the processor, when executing the program, implements the above-mentioned high-precision numerical analysis method for wave propagation in nonlinear media, including:

[0101] Construct model equations for wave propagation in nonlinear media;

[0102] The model equation is discretized using the finite difference method to obtain a high-precision difference format for the approximate solution of the model equation;

[0103] Designing a numerical algorithm for solving the model equations based on a high-precision difference scheme;

[0104] Numerical simulation and analysis of wave propagation in nonlinear media.

[0105] Example 4:

[0106] A computer-readable storage medium having a computer program stored thereon, wherein when the program is executed by a processor, the program implements the above-mentioned high-precision numerical analysis method of wave propagation in nonlinear media, comprising:

[0107] Construct model equations for wave propagation in nonlinear media;

[0108] The model equation is discretized using the finite difference method to obtain a high-precision difference format for the approximate solution of the model equation;

[0109] Designing a numerical algorithm for solving the model equations based on a high-precision difference scheme;

[0110] Numerical simulation and analysis of wave propagation in nonlinear media.

[0111] Those skilled in the art will appreciate that the modules or steps of the present disclosure described above can be implemented using a general-purpose computer device. Alternatively, they can be implemented using program code executable by a computing device, which can then be stored in a storage device and executed by the computing device. Alternatively, they can be fabricated into separate integrated circuit modules, or multiple modules or steps can be fabricated into a single integrated circuit module for implementation. The present disclosure is not limited to any specific combination of hardware and software.

[0112] The foregoing description is merely a preferred embodiment of the present application and is not intended to limit the present application. Persons skilled in the art will readily appreciate that various modifications and variations are possible. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present application shall be included within the scope of protection of the present application.

[0113] Although the above describes the specific implementation methods of the present disclosure in conjunction with the accompanying drawings, it is not intended to limit the scope of protection of the present disclosure. Those skilled in the art should understand that on the basis of the technical solution of the present disclosure, various modifications or variations that can be made by those skilled in the art without creative work are still within the scope of protection of the present disclosure.

Claims

1. A high-precision numerical analysis method for wave propagation in nonlinear media, characterized by: The following steps are involved: Construct model equations for wave propagation in nonlinear media; The model equation is discretized using the finite difference method to obtain a high-precision difference format for the approximate solution of the model equation; Designing a numerical algorithm for solving the model equations based on a high-precision difference scheme; Numerical simulation and analysis of wave propagation in nonlinear media; The model equation for wave propagation in nonlinear media is constructed as follows: Assume that the wave propagates in the two-dimensional region Ω = [a1, b1] × [a2, b2], where the wave amplitude is the function φ(x, y, t), the potential energy in the nonlinear medium is represented by f(x, y), and the nonlinear coefficient during propagation is β; then the model equation satisfied by the wave amplitude during propagation in the nonlinear medium is: Where Δ is the Laplace operator, α is a constant, i is an imaginary unit, x and y are spatial direction variables, t is a time variable, and T is the total time required for the calculation to arrive; the displacement and velocity of the wave propagation at the initial moment are: where ψ(x,y) and is a smooth function; In order to simulate the propagation of waves in nonlinear media, the model equation needs to be discretized, so a new variable function is introduced: Formula (1)-(2) can be rewritten as The finite difference method is used to discretize the model equation and obtain a high-precision difference format for the approximate solution of the model equation. Specifically, the time step is defined as The calculation time is evenly divided into N intervals. At the same time, the calculation domain is evenly meshed. Given two positive integers M x and M y , define the spatial step size as and That is, [a1, b1] × [a2, b2] is evenly divided into M x ×M y intervals, where a1, b1, a2, and b2 are the endpoints of the intervals; each calculation grid point is (x j ,y l ,t n ), where x j =a1+jh x ,y l =a2+lh y ,t n =nτ; In addition, the grid function is defined The difference operator is defined as: remember and is the point (x j ,y l ,t n ), using the Crank-Nicolson method for time and the sixth-order difference operator for space, we can obtain the following high-precision difference format: The periodic boundary condition is expressed as k is the kth time step, δ t is the time-centered difference operator; Express the above high-precision difference format in matrix form: where Φ, P, and f are Φ j,l ,P j,l ,f j,l The matrix formed by all discrete points in the calculation area, and Here circ indicates that the matrix is ​​a circulant matrix. Represents the matrix A x The inverse matrix of A y The inverse matrix of ; the accuracy of this high-order format is The numerical algorithm for solving the model equations includes a time extrapolation algorithm that raises the high-precision difference scheme to sixth order using the Richardson extrapolation technique as follows: The numerical algorithm for solving the model equations includes a fast solution algorithm that equates the matrix form of the high-precision difference format to in Circulant matrix B x , B y , G x and G y Diagonalizable in They are the circulant matrices B x ,B y ,G x ,G y The diagonal matrix formed by the eigenvalues, T x ,T y is a reversible matrix composed of all eigenvectors of the circulant matrix; further arranging equations (13)-(14), we get: According to the diagonalization of the circulant matrix, Equation (17) can be rewritten as Multiply both sides of equation (18) by T x , right multiplication set up get in and 1 is a matrix whose elements are all 1; in specific numerical calculations, the following iterative format is used Where . / and .* represent element-wise division and multiplication of matrices. When ||Φ k+1(s+1) -Φ k+1(s) || ∞ ≤10 -12 When Φ k+1 =Φ k+1(s+1) , we get Φ k+1 Then, calculate P by formula (16) k+1 .

2. The high-precision numerical analysis method for wave propagation in nonlinear media according to claim 1, characterized in that: The numerical algorithm for solving the model equations includes a time-advancing algorithm, which is based on a time extrapolation algorithm and a fast solution algorithm, combined with a high-precision difference format. The entire time-advancing process for solving equations (1)-(2) is given as follows: Step 1: Take the time step τ=λ·max(h x ,h y ),λ∈R + , get the initial value and Step 2: Iterate formula (20) to get the function value of the first time layer Step 3: Obtain the function value of the first time layer through formula (16) Step 4: Repeat Step 2-3 until the final moment, and get and Step 5: Re-take the time step τ=λ·max(h x ,h y ) / 2, repeat Step 2-Step 4 to get and Step 6: Re-take the time step τ = λ·max(h x ,h y ) / 4, repeat Step 2-Step 4 to get and Step 7: Finally, we can get the new and 3. A high-precision numerical analysis system for wave propagation in nonlinear media, characterized by: include: Model acquisition module, which constructs the model equations of wave propagation in nonlinear media; The format design module uses the finite difference method to discretize the model equation and obtain a high-precision difference format that approximates the solution of the model equation; A solution algorithm design module, which designs a numerical algorithm for solving the model equation based on a high-precision difference format; Analysis module, which performs numerical simulation and analysis on wave propagation in nonlinear media; In the model acquisition module, assume that a wave propagates in a two-dimensional region Ω = [a1, b1] × [a2, b2], where the wave amplitude is a function φ(x, y, t), the potential energy in the nonlinear medium is represented by f(x, y), and the nonlinear coefficient during propagation is β. Then, the model equation satisfied by the wave amplitude during propagation in the nonlinear medium is: Where Δ is the Laplace operator, α is a constant, i is an imaginary unit, x and y are spatial direction variables, t is a time variable, and T is the total time required for the calculation to arrive; the displacement and velocity of the wave propagation at the initial moment are: where ψ(x,y) and is a smooth function; In order to simulate the propagation of waves in nonlinear media, the model equation needs to be discretized, so a new variable function is introduced: Formula (1)-(2) can be rewritten as The finite difference method is used to discretize the model equation and obtain a high-precision difference format for the approximate solution of the model equation. Specifically, the time step is defined as The calculation time is evenly divided into N intervals. At the same time, the calculation domain is evenly meshed. Given two positive integers M x and M y , define the spatial step size as and That is, [a1, b1] × [a2, b2] is evenly divided into M x ×M y intervals, where a1, b1, a2, and b2 are the endpoints of the intervals; each calculation grid point is (x j ,y l ,t n ), where x j =a1+jh x ,y l =a2+lh y ,t n =nτ; In addition, the grid function is defined The difference operator is defined as: remember and is the point (x j ,y l ,t n ), using the Crank-Nicolson method for time and the sixth-order difference operator for space, we can obtain the following high-precision difference format: The periodic boundary condition is expressed as k is the kth time step, δ t is the time-centered difference operator; Express the above high-precision difference format in matrix form: where Φ, P, and f are Φ j,l ,P j,l ,f j,l The matrix formed by all discrete points in the calculation area, and Here circ indicates that the matrix is ​​a circulant matrix. Represents the matrix A x The inverse matrix of A y The inverse matrix of ; the accuracy of this high-order format is The numerical algorithm for solving the model equations includes a time extrapolation algorithm that raises the high-precision difference scheme to sixth order using the Richardson extrapolation technique as follows: The numerical algorithm for solving the model equations includes a fast solution algorithm that equates the matrix form of the high-precision difference format to in Circulant matrix B x , B y , G x and G y Diagonalizable in They are the circulant matrices B x ,B y ,G x ,G y The diagonal matrix formed by the eigenvalues, T x ,T y is a reversible matrix composed of all eigenvectors of the circulant matrix; further arranging equations (13)-(14), we get: According to the diagonalization of the circulant matrix, Equation (17) can be rewritten as Multiply both sides of equation (18) by T x , right multiplication set up get in and 1 is a matrix whose elements are all 1; in specific numerical calculations, the following iterative format is used Where . / and .* represent element-wise division and multiplication of matrices. When ||Φ k+1(s+1) -Φ k+1(s) ||≤10 -12 When Φ k+1 =Φ k+1(s+1) , we get Φ k+1 Then, calculate P by formula (16) k+1 .

Images