JP3394989B2 - Numerical calculation method and numerical calculation device, and recording medium storing numerical calculation program - Google Patents

Description

Detailed Description of the Invention

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a numerical calculation technique for physical quantities such as pressure and temperature.
In particular, it belongs to a numerical calculation technique using an iterative solution algorithm .

[0002]

In the past, iterative approximation algorithm has been used in control techniques for physical quantities such as pressure and temperature. The iterative approximation algorithm is an algorithm for obtaining an optimal solution by numerical calculation while repeatedly performing iterative approximation based on various optimal control theories.

However, when the iterative approximation algorithm is used, the number of iterations of the calculation is increased when obtaining the solution of each parameter to be controlled, and therefore the control is performed at high speed and with high accuracy in terms of calculation time. Has become difficult. Therefore, it is desired to improve the convergence of the iterative solution algorithm so that it can be applied to control technologies in various fields.

In view of the above-mentioned problems, the present invention provides a numerical counter using an iterative solution algorithm with higher convergence than before.
It is an object to provide a calculation method, a numerical calculation device, and a recording medium recording a numerical calculation program.

[0005]

In order to solve the above-mentioned problems, the solving means taken by the invention of claim 1 is that the coefficient matrix when the partial differential equation to be satisfied by the physical quantity U is discretized as A, Assuming that the homogeneous term (source term) is f, a computer is used to solve AU = f, and a numerical value of the physical quantity U is calculated. An initial value U ⁰ of the physical quantity U is set and the number of repetitions m is set. 0 is given as the initial value, 0 is given as the initial value of the perturbation amount φ, and as the initial value r ⁰ of the residual r, (f−A · U ⁰ ).
And incrementing the number of repetitions m,
A · φ = r ^m predicted approximation Ψ ^m with internal solver
The step of iterative calculation by the calculation unit of
An optimized approximation value φ ^m that minimizes the L 2 norm of the residual r ^m is obtained from the predicted approximation value Ψ ^m by the optimization routine by the second computing unit, and as an approximate solution U ^{m + 1} (U ^m + φ ^m ). And a step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} until the approximate solution U ^m converges, and the step of obtaining the predicted approximate value Ψ ^m is as follows. A variable κ ^m, which represents the residual resection rate for each number of repetitions ^m, is a predetermined value K
When it exceeds, the iterative calculation is finished.

In the invention of claim 2 , the above-mentioned claim
The variable κ ^m in the numerical calculation method of 1 is

[Equation 1] Shall be given by the formula.

Further, the solving means taken by the invention of claim 3 is that, when the coefficient matrix when the partial differential equation to be satisfied by the physical quantity U is discretized is A and the inhomogeneous term (source term) is f, As an apparatus for solving U = f and numerically calculating the physical quantity U, an initial value U ⁰ of the physical quantity U is set, 0 is given as an initial value for the number of repetitions m, 0 is given as an initial value of the perturbation quantity φ, and the remaining (F−A · U ⁰ ) is set as the initial value r ⁰ of the difference r, and A · φ
= R ^m predicted approximate value Ψ ^m by iterative calculation by the first arithmetic unit having an internal solver, and residual r ^m
Of the optimized approximated value φ ^m that minimizes the L2 norm of
obtained from ^m, the approximate solution U ^{m + 1} gives (U ^m + phi ^m), until the residual r ^{m + 1} as (r ^m -A · φ ^m) approximate solution U ^m and providing the converge, This step is repeatedly executed, and the step of obtaining the predicted approximate value Ψ ^m ends the iterative calculation when the variable κ ^m representing the residual cutoff rate for each number of repetitions m exceeds a predetermined value K.

In the invention of claim 4 , the above-mentioned claim
The variable κ ^m in the numerical computer of 3 is

[Equation 1] Shall be given by the formula.

Further, in the solution means taken by the invention of claim 5 , the coefficient matrix when the partial differential equation to be satisfied by the physical quantity U is discretized in the computer is A, and the inhomogeneous term (source term) is f. At this time, the initial value U ⁰ of the physical quantity U is set as a recording medium in which a numerical calculation program for numerically calculating the physical quantity U is solved, and A · U = f is set. 0 is given as the initial value of the perturbation amount φ, (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r, and A · φ is increased while incrementing the number of repetitions m.
= R ^m predicted approximate value Ψ ^m by iterative calculation by the first arithmetic unit having an internal solver, and residual r ^m
Of the optimized approximated value φ ^m that minimizes the L2 norm of
obtained from ^m, the approximate solution U ^{m + 1} gives (U ^m + φ ^m), repeated until the residual r ^{m + 1} as (r ^m -A · φ ^m) approximate solution U ^m and providing the converges In the step of obtaining the predicted approximate value Ψ ^m , the computer performs a process of ending the iterative calculation when the variable κ ^m representing the residual cutoff rate for each number of iterations m exceeds a predetermined value K. It is a record of a numerical calculation program to be executed.

According to the invention of claim 6 , in the recording medium on which the program for numerical calculation of claim 5 is recorded, the variable κ ^m used in the program for numerical calculation is:

[Equation 1] Shall be given by the formula.

[0011]

BEST MODE FOR CARRYING OUT THE INVENTION An embodiment of the present invention will be described below with reference to the drawings.

FIG. 1 is a block diagram showing a configuration of a control device according to an embodiment of the present invention, showing an example in which a numerical calculation technique according to the present invention is applied to automatic vehicle traveling control. .
In FIG. 1, it is assumed that the vehicle includes a circuit that automatically travels to a destination within an allowable error range on the track set by the target value / route setting unit 11. During traveling, the position information detector 12 detects traveling information such as position / direction information and speed. The display device 18 displays information such as set values and running conditions on a monitor in real time.

The deviation φ between the vehicle position and the track is constantly monitored by the information management device 13. Then, when this deviation φ exceeds the allowable value δmax, the travel information at that time is used as a starting value, and the optimization parameter control of the direction angle, the angular acceleration, the speed, the acceleration, etc. A simulation is performed using the algorithm according to the present invention.

FIG. 2 is a flow chart showing a numerical calculation method according to an embodiment of the present invention.

(Formulation of residual equation) Let U be the physical quantity to be obtained, A is the coefficient matrix when the partial differential equation to be satisfied by the physical quantity is discretized, and f is the inhomogeneous term (source term). The equation to be solved is represented by AU = f (1). This formula is generally used in the SOR method, ADI
It can be solved by a successive approximation method such as a method.

In the present invention, the following solution method is used.

The solution U∞ of the equation (1) is expressed as follows by the approximate solution U and the perturbation amount (difference between the true solution) and φ. U∞ = U + φ (2) In the solution according to the present invention, the perturbation amount φ is obtained and the approximate solution U is modified to obtain the true solution U∞ of the equation (1).

Here, the residual r for the approximate solution U is defined as follows.

R = f−A · U (3) From the equations (1) to (3), A · (U + φ) = f∴A · φ = f−A · U = r Therefore, the perturbation amount φ is obtained. In order to do so, solve the following equation. A · φ = r (4)

(Algorithm of Residual Exclusion Method) In the algorithm according to the present invention, the predicted approximate value Ψ is obtained by repeating the smallest unit with the steepest convergence gradient by the ADI method or the like, instead of obtaining the convergent solution of the equation (4). Is calculated, and the calculated predicted approximate value Ψ is supplied to the optimization control routine. And
The calculation of the predicted approximate value Ψ is repeatedly executed until the execution result of the optimization control routine satisfies the predetermined condition.

<Optimization Control Routine> Now, assuming that the number of iterations is m, a synthetic perturbation φ ^m that minimizes the L2 norm of the residual and a new approximate solution U ^{m + 1} are defined as follows.

[Equation 2] U ^{m + 1} = U ^m + φ ^m (6) α _l (l = 1, 2, 3, ..., L) is a residual minimization coefficient, which is a constant obtained by a calculation method described later. The residual is minimized as follows so that the L2 norm of the residual becomes smaller.

When the new approximate solution U ^{m + 1} is expressed by the equation (6), the residual error r ^{m + 1} for it is expressed by the following equation (3).

R ^{m + 1} = f−A · U ^{m + 1} (7) Substituting equation (6) into equation (7) and using equations (3) and (5),

[Equation 3]

The L2 norm ‖r ^{m + 1} ‖ of the residual r ^{m + 1} with respect to the approximate solution U ^{m + 1} is, in the case of three dimensions, expressed as the square root of the sum of all interior points of (r ^{m + 1} ) ^2. It is given in (9).

[Equation 4]

In order to minimize the L2 norm ‖r ^{m + 1} ‖ of this residual r ^{m + 1} , the residual minimization coefficient α _l ( l =
1, 2, 3, ..., L ) are determined by the least squares method. That is,

[Equation 5] Equation (10) becomes an L-element simultaneous equation of α, and α _l is determined by numerically solving it. When α _l is determined, φ ^m is obtained from the equation (5) and U ^{m + 1} is obtained from the equation (6).
A method of calculating the residual minimization coefficient α _l will be described later.

By repeating the same routine while incrementing m, the solution converges to the solution U ∞ that makes the L2 norm of the residual of equation (9) zero or minimum.

The numerical calculation method according to this embodiment will be described with reference to the flowchart shown in FIG.

First, in step S1, an initial value U ⁰ of the target physical quantity U is set. Then, in step S2, a 0 as an initial value in number of repetitions m, give 0 as an initial value of the perturbation quantity phi, as the initial value r ⁰ of the residual r (f
-A · U ⁰ ) is set.

Thereafter, steps S3 to S are performed while incrementing the number of iterations m until the approximate solution U ^m converges.
9 is repeatedly executed.

In steps S3 to S8, the ADI method, the SO
R method, by internal solver to perform successive approximation method such as CG method, the approximate solution of A · φ = r ^m (predicted approximate value) is obtained [psi ^m. However, if this equation is to be solved by the existing iterative approximation method, it will take many iterations for practical problems to converge, and inevitably the computer time will be enormous.

Therefore, in the present embodiment, the upper limit N of the number of iterations of the sequential calculation is set (steps S7 and S8), and at the same time, the variable κ representing the reduction rate of the residual, that is, the residual cutoff rate.
^{When m} is introduced (step S5) and the variable κ ^m exceeds a predetermined value K, the iterative calculation is finished (step S6).

When the number of iterations of the sequential calculation is fixed,
The rate of reduction of residuals (residual cut-off rate) fluctuates greatly with each iteration m of the loop until the convergence determination. This variation is due to the degree of approximation of the predicted approximation value Ψ ^m at each iteration m. Therefore, the number of loop iterations until the solution converges may increase more than necessary.

On the other hand, in this embodiment, the iterative calculation is terminated when the variable κ ^m representing the residual cutoff rate exceeds the predetermined value K in the step of obtaining the predicted approximate value Ψ ^m . The residual excision rate for each iteration m of the loop can be maintained above a certain value, and therefore the number of loop iterations required until the solution converges can be reduced.

Here, the variable κ ^m representing the residual cutoff rate is defined as follows.

[Equation 6]

This variable κ ^m becomes 1.0 because the accuracy of the predicted approximation value Ψ ^m becomes higher as the iterative calculation is repeated.
Approach endlessly. Therefore, by performing the iterative calculation until the value of the variable κ ^m exceeds the arbitrarily set predetermined value K, the residual cutoff rate for each loop iteration number m becomes substantially constant.

In steps S3 to S8, the predicted approximation value Ψ
⁰ has been obtained, in step S9, by using the predicted approximate value [psi ^0, the optimization routine, as the next time of the residual r ¹ is minimized, approximated phi ⁰ optimized. With this optimized approximation value φ ⁰ , the approximation value U ¹ of the first-order physical quantity and the residual error r ¹ of the first order are obtained from the equations (2) and (3) by the following equation. U ¹ = U ⁰ + φ ⁰ (12) r ¹ = f−A · U ¹ (13) Substituting the formula (12) into the formula (13), r ¹ = f−A · (U ⁰ + φ ⁰ ) = F−A · U ⁰ −A · φ ⁰ = r ⁰ −A · φ ⁰ (14) Therefore, the first-order residual r ¹ is obtained by the equation (14).

In step S11, m is incremented, the process returns to step S3, and the same calculation is performed. By repeating such processing, (U ² , r ² ), (U ³ , r ³ ) ... (U ^m , r ^m ) ... are sequentially obtained, and the residual norm ‖r ^m ‖ decreases. .

[0038] When the number of repetitions of m, the equation for determining the predicted approximate value [psi ^m in step S4, A · φ = r ^m, and the step S9, the optimized An approximation phi ^m by the optimization control routine , U ^{m + 1} = U ^m + φ ^m r ^{m + 1} = r ^m −A · φ ^m .

The solution is when ‖r ^m ‖ = 0, but the final residual rfinal may remain depending on the boundary conditions. That is, when all the boundary conditions are given by the derivative of the physical quantity (Neumann problem), there is a relation (Neumann constraint condition) that should be mathematically satisfied between the source term f and the boundary condition. The difference norm becomes ‖r ^m ‖ = ‖r ^m-1 ‖ = ‖r ^m-2 ‖ ...

In the other solution methods, in the case of such a Neumann problem, the final residual rfinal does not converge as it is.
, The source term f or the boundary condition is modified to obtain a convergent solution. This is changing the problem. However, the solution method according to the present embodiment has a characteristic that the final residual rfinal is automatically determined and a convergent solution is obtained even without such correction.

(How to Obtain Residual Minimization Coefficient α _l ) Now, in order to simplify the notation, δ is introduced as follows. Letting δ ^m = Ψ ^m δ ^{m-l + 1} = φ ^{m-l + 1} (l = 2, 3, ..., L), the optimized approximation φ ^m is

[Equation 7] Is expressed as Therefore,

[Equation 8] Expressed without the Σ symbol, r ^{m + 1} = r ^m − (α ₁ A · δ ^m + α ₂ A · δ ^m-1 + α ₃ A · δ
^m-2 + ... + α _L A · δ ^{m-L + 1} ) .

The residual minimization coefficient α _l is obtained by giving the condition that the (m + 1) th order residual norm ‖r ^{m + 1} ‖ (square root of sum of squares) is minimized by using, for example, the least squares method. ,Desired.

[Equation 9] Differentiating this by α _l ,

[Equation 10]

By solving this equation for α _l ,
The residual minimization coefficient α _l can be obtained. If you perform differentiation,

[Equation 11] Becomes Here, l ′ = 1, 2, ..., L. By solving this simultaneous linear equation, the residual minimization coefficient α _l
Can be asked.

For example, when L = 3, the above equation becomes

[Equation 12] It becomes a three-dimensional simultaneous linear equation like. _{^{Then, φ m = α 1 δ m}} + α 2 δ m-1 + α 3 δ m-2 r m + 1 = r m - (α 1 A · δ m + α 2 A · δ m-1 + α 3 A · δ m ^-2 ) = r ^m −A · φ ^m U ^{m + 1} = U ^m + φ ^m .

In the traveling control of the vehicle shown in FIG. 1, first, the predicted correction amount calculation unit 14 repeatedly calculates the predicted correction amount Ψ ^m based on the current traveling information. Here, iterative calculation is performed until the residual excision rate κ ^{m for} each number of iterations exceeds a predetermined value K. Next, the optimum addition correction amount calculation unit 15 calculates the m-th order addition correction amount φ ^m . The solution obtained by the predicted correction amount calculation unit 14 and the optimum addition correction amount calculation unit 15 is finally determined to be converged by the correction amount convergence determination unit 16, and the optimum control method for returning to a predetermined trajectory is determined. It The determined control information is the direction control unit 17
And the traveling route of the vehicle is corrected.

Although the present invention has been described here by taking the traveling control of the vehicle as an example, it goes without saying that the present invention can be applied to various controls in other fields.

The above numerical calculation method can be easily realized by an apparatus equipped with a computer that executes a program for realizing the method. Also,
It can be realized by recording a program for realizing the method on a computer-readable recording medium and causing a computer to execute the program recorded on the recording medium.

Here, a supplementary explanation will be given regarding the comparison with the prior art. Conventionally, the technologies used to control equipment include
There are proportional control and fuzzy control. Taking position control as an example, in the case of proportional control, basically, the comparison between the current position and the set position is binarized and judged, and the deviation is tried to approach zero as much as possible. Actions are frequent. Therefore, there are drawbacks that the operation is not stable and the disturbance cannot be dealt with flexibly.

As a technique for solving such a drawback, a control technique employing a fuzzy theory has been developed. this is,
The state quantity is "largely deviated" or "a little deviated"
Judgment is "fuzzy" such as "subtle deviation", and control is performed based on "fuzzy" rules such as "correct largely", "correct a little" or "just" for these states. Is. With this technology, "smooth control" similar to AI control was realized, but the number of rules increases in proportion to the amount of information and control parameters, which makes setting and tuning work troublesome or subjective setting. However, there is also a drawback that the control results vary because of the possibility of

The present invention is a technique used to determine the optimum means for controlling a controlled object to a target state, and its basis is a nonlinear equation corresponding to a control parameter. For this reason, it is possible to adapt to the case where the target state is not clear numerically depending on the target, and it also has the function of automatically tuning the optimum state.

(Other Application Example 1) FIG. 3 is a diagram showing another application example of the algorithm according to the present invention. FIG. 3 is a flow showing the application of the two-dimensional or three-dimensional elliptic equation method to the grid generation method as an application example to the numerical calculation.

The symbols used in FIG. 3 have the following meanings. U: Coordinate value (x, y, z) of the grid point on the physical surface P: Control function of grid generation A: Linearized coefficient matrix Ψ ^(m) : Predicted correction value φ ^(m) at m-th step: Optimized modified value αl in m-th step: undetermined coefficient matrix r ^(m) : residual L in m-th step L: number of undetermined coefficient matrix

In step S21, the coefficient matrix of the discretized elliptic equation system is linearized based on the values in the m-th step. In step S22, the number of iterations is controlled by the κ technique to obtain the predicted correction amount Ψ ^(m) in the (m + 1) th step. In step S23,
An undetermined coefficient matrix αi for minimizing the residual of the (m + 1) th step is obtained. In step S24, the predicted correction value Ψ ^(m) , the optimized correction before the (m−1) th order
Based on the values (φ ^(m-1) , ..., φ ^{(m-L + 1)} ) and the undetermined coefficient matrix, the optimized correction value φ ^{(m + 1)} at the ^{(m + 1) th step} is obtained. In step S25, the optimized
The coordinate value of the grid point on the physical surface is corrected based on the correction value φ ^{(m + 1)} .

FIG. 4 shows a lattice generated by this application example,
FIG. 5 shows a lattice generated by the SOR method as a comparative example. In FIGS. 4 and 5, a two-dimensional grid is generated around an elliptic cylinder placed in the air flow. The number of grid points is 90 ×
90. In this application example, the ADI method is used as the internal solver of the residual excision method.

FIG. 6 is a graph for comparing convergence speeds. In FIG. 6, a1 and a2 are X in this application example.
Direction and Y direction residuals, and b1 and b2 are X direction and Y direction residuals in the comparative example. As can be seen from FIG. 6, when the algorithm according to the present invention is used, S
It converges more than 6 times faster than the OR method.

(Other Application Example 2) FIG. 7 is a diagram showing another application example of the algorithm according to the present invention. FIG. 7 is a flow chart showing an example in which the residual ablation method is applied to the pressure calculation unit of the MAC method of fluid analysis, as an application example to numerical calculation.

The symbols used in FIG. 7 have the following meanings. t: time U: flow velocity P: pressure Ψ ^m : predicted correction amount in m-th step φ ^m : optimized correction amount in m-th step α: optimization coefficient r ^m : residual in m-th step

S31 is an initialization part including calculation of the initial residual r ⁰ , S32 is a calculation part of the estimated approximate value Ψ ^m , and S33 is κ ^m.
Calculation unit, the loop S34 includes S32 and S33 which are controlled repetition number by kappa ^m calculated in S33, the optimization coefficient S35 alpha, optimized correction quantity phi ^m and residual r in the next step It is a calculation unit of ^{m + 1} .

In order to evaluate the effect of this application example, M
Two-dimensional analysis around a cylinder placed in a uniform flow was performed by the AC method. The SOR method was used as a comparative example. The number of lattice points is 60 × 60, and the main flow velocity is 5 cm / s. In this application example, the SOR method is used as the internal solver of the residual cutting method.

FIG. 8 is a graph showing the relationship between the number of time steps and the CPU time. In FIG. 8, a shows the result by this application example, and b shows the result by the SOR method as a comparative example. FIG. 9 is a graph showing the relationship between the number of steps and the relative residual / relative error. In FIG. 9, a1 and a2 are results of this application example, and b1 and b2 are SO as a comparative example.
The result by the R method is shown. From FIGS. 8 and 9, it can be seen that according to this application example, the speed is increased by a factor of 5 compared to the SOR method.

[0061]

As described above, according to the present invention, when the variable κ ^m representing the residual cut rate for each number of iterations m exceeds the predetermined value K in the step of obtaining the predicted approximate value Ψ ^m , the iterative calculation is performed. Since the process is completed, the residual cutoff rate for each m iterations of the loop can be maintained at a certain value or more, and thus the value of m required for convergence can be reduced. Therefore,
Since the convergence of the iterative solution algorithm is improved, an example
For example, when used for control, high-speed and high-precision control becomes easy.

[Brief description of drawings]

FIG. 1 is a block diagram showing a configuration of a control device according to an embodiment of the present invention.

FIG. 2 is a flowchart showing a numerical calculation method according to an embodiment of the present invention.

FIG. 3 is a flowchart of an application example of the present invention, in which the algorithm according to the present invention is applied to a grid generation method.

FIG. 4 is a diagram showing a lattice generated by the application example of FIG. 3;

FIG. 5 is a diagram showing a lattice generated by a SOR method as a comparative example.

FIG. 6 is a graph for comparing the convergence speeds of the application example of FIG. 3 and the comparative example.

FIG. 7 is an application example of the present invention, which is a flowchart when the algorithm according to the present invention is applied to fluid analysis.

8 is a graph for comparing the CPU times of the application example of FIG. 7 and the comparative example.

9 is a graph for comparing the convergence speeds of the application example of FIG. 7 and the comparative example.

─────────────────────────────────────────────────── ─── Continuation of the front page (72) Inventor Atsuhiro Tamura 7-44-1, Jindaiji East, Chofu City, Tokyo Inside the Aerospace Research Institute (72) Inventor Koji Taniguchi 1-1835 Edobori, Nishi-ku, Osaka City Higobashi IP Building Vinus Co., Ltd. (56) Reference JP-A-6-95707 (JP, A) JP-A-8-5450 (JP, A) (58) Fields investigated (Int.Cl. ⁷ , DB name) G05B 11/00-13/04 G06F 15 / 20,15 / 328

Claims (6)

(57) [Claims] 1. When a coefficient matrix for discretizing a partial differential equation to be satisfied by a physical quantity U is A and an inhomogeneous term (source term) is f, a computer is used to solve A · U = f, A method for numerically calculating a physical quantity U, wherein an initial value U ⁰ of the physical quantity U is set, 0 is given as an initial value for the number of iterations m, 0 is given as an initial value of the perturbation quantity φ, and an initial value of the residual r is given. As r ⁰ (f−A ·
Set U ^0), incrementing the number of repetitions m, the predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a step of determining by an iterative calculation, the residual r ^m The optimized approximation value φ ^m that minimizes the L2 norm of is calculated from the predicted approximation value Ψ ^m by the second calculation unit by the optimization routine, and is calculated as (U ^{m +1)} as the approximation solution U ^{m + 1.}
+ Giving phi ^m), and providing a residual r ^{m + 1} a ^{(r m -A · φ m)} , until the approximate solution U ^m converges, which repeatedly executed, the predicted approximate value [psi ^m The step of obtaining is a numerical calculation method characterized in that the iterative calculation is terminated when the variable κ ^m representing the residual cutting rate for each number of repetitions m exceeds a predetermined value K. 2. The numerical calculation method according to claim 1 , wherein the variable κ ^m is Numerical calculation method characterized by being given by the formula. 3. A numerical calculation of the physical quantity U by solving A · U = f where A is a coefficient matrix when the partial differential equation to be satisfied by the physical quantity U is discretized and f is an inhomogeneous term (source term). An apparatus for performing the following, wherein an initial value U ⁰ of the physical quantity U is set, 0 is given as an initial value to the number of iterations m, 0 is given as an initial value of the perturbation quantity φ, and an initial value r ⁰ of the residual r is given as (f -A
Set U ^0), incrementing the number of repetitions m, the predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a step of determining by an iterative calculation, the residual r ^m The optimized approximation value φ ^m that minimizes the L2 norm of is calculated from the predicted approximation value Ψ ^m by the second calculation unit by the optimization routine, and is calculated as (U ^{m +1)} as the approximation solution U ^{m + 1.}
+ Giving phi ^m), and providing a residual r ^{m + 1} a ^{(r m -A · φ m)} , until the approximate solution U ^m converges, which repeatedly executed, the predicted approximate value [psi ^m The numerical calculation device characterized in that the step of obtaining is to end the iterative calculation when the variable κ ^m representing the residual ablation rate for each number of repetitions m exceeds a predetermined value K. 4. The numerical calculation device according to claim 3 , wherein the variable κ ^m is Numerical calculation device characterized by being given by the formula. 5. When a coefficient matrix for discretizing a partial differential equation to be satisfied by a physical quantity U is A and a non-homogeneous term (source term) is f, the computer is caused to solve A · U = f A recording medium in which a numerical calculation program for numerical calculation of U is recorded, wherein an initial value U ⁰ of a physical quantity U is set, 0 is given as an initial value of the number of repetitions m, and 0 is given as an initial value of the perturbation quantity φ. , As the initial value r ⁰ of the residual r (f−A ·
Set U ^0), incrementing the number of repetitions m, the predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a step of determining by an iterative calculation, the residual r ^m The optimized approximation value φ ^m that minimizes the L2 norm of is calculated from the predicted approximation value Ψ ^m by the second calculation unit by the optimization routine, and is calculated as (U ^{m +1)} as the approximation solution U ^{m + 1.}
+ Giving phi ^m), and providing a residual r ^{m + 1} a ^{(r m -A · φ m)} , until the approximate solution U ^m converges, which repeatedly executed, the predicted approximate value [psi ^m The step of obtaining is a recording medium recording a numerical calculation program for causing a computer to execute the processing for ending the iterative calculation when the variable κ ^m representing the residual cutting rate for each number of repetitions m exceeds a predetermined value K. . 6. A recording medium recording the numerical calculation program according to claim 5 , wherein the variable κ ^m is A recording medium recording a numerical calculation program, which is given by