US20090063597A1

US20090063597A1 - Numerical analysis device and numerical analysis program

Info

Publication number: US20090063597A1
Application number: US11/887,003
Authority: US
Inventors: Hiroaki Sono; Nobuhiro Yamada; Haruya Kitagawa; Tsuyoshi Nomura
Original assignee: Hokuriku Electric Power Co
Current assignee: Hokuriku Electric Power Co
Priority date: 2005-03-25
Filing date: 2006-03-22
Publication date: 2009-03-05
Also published as: US20120047192A1; JPWO2006103996A1; US8793293B2; WO2006103996A1

Abstract

A second derivative of a second-order differential equation is calculated at a reference variable value. The second derivative is multiplied by an analytical small variable value, the first derivative at the reference variable value is added, and a result is output as a first derivative after an increment of the analytical small variable value. The first derivative after an increment of the analytical small variable value is multiplied by the analytical small variable value, a physical value at the reference variable value is added, and a result is output as a physical value after an increment of the analytical small variable value.

Description

TECHNICAL FIELD

The present invention relates to a numerical analysis device and a numerical analysis program that can be applied to analyze various aspects of motion in natural phenomena and mechanical movements.

BACKGROUND ART

Traditionally Euler's method, in which the displacement and the velocity at the next time point are calculated by using an initial displacement, initial velocity, and initial acceleration in an interval, has been widely used for the numerical analysis using ordinary differential equation (see patent documents 1 and 2).
Patent document 1: Publication of unexamined patent No. 2004-062676
Patent document 2: Publication of unexamined patent No. 2002-258933

DISCLOSURE OF INVENTION

Problem to be Solved by the Invention

Euler's method, however, has a problem that a relatively large error occurs in an analysis result when processing steps are repeated with a small time increment Δt. Furthermore, there is another problem that when a very small time increment Δt is applied to avoid this error so much processing time to obtain an analysis result is needed that the method is eventually impractical to apply.
The present invention is made considering these problems and aims at providing numerical analysis devices and numerical analysis programs that can simply perform a numerical analysis process within a short time period without causing a large error in the analysis result.

Means for Solving Problem

To solve the above problems and to achieve the object, a numerical analysis device according to the present invention performs a numerical analysis of second-order differential equation for a target object to be analyzed. The numerical analysis device includes a setting unit that performs a setting of parameters including an initial condition of the second-order differential equation and a setting of an analytical small variable value for the numerical analysis; a second-derivative calculation unit that calculates a second derivative of the second-order differential equation at a reference variable value that is a variable at the start of the numerical analysis; a first-derivative calculation unit that performs a multiplication of the second derivative calculated by the second-derivative calculation unit by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputs a value obtained by the addition as a first derivative after an increment of the analytical small variable value; a physical-value calculation unit that performs a multiplication of the first derivative calculated by the first-derivative calculation unit by the analytical small variable value and an addition of a physical value of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputs a value obtained by the addition as a physical value after an increment of the analytical small variable value; and an operation control unit that repeats a control of setting a variable value obtained by adding the analytical small variable value to the reference variable value as a new reference variable value and causing the second-derivative calculation unit, the first-derivative calculation unit, and the physical-value calculation unit to calculate the second derivative, the first derivative, and the physical value, using the new reference variable value, respectively.
Furthermore, a numerical analysis device according to the present invention performs a numerical analysis of second-order differential equation for a target object to be analyzed. The numerical analysis device includes a setting unit that performs a setting of parameters including an initial condition of the second-order differential equation and a setting of an analytical small variable value for the numerical analysis; a second-derivative calculation unit that calculates a second derivative of the second-order differential equation at a reference variable value that is a variable at the start of the numerical analysis; a first-derivative calculation unit that performs a multiplication of the second derivative calculated by the second-derivative calculation unit by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputs a value obtained by the addition as a first derivative after an increment of the analytical small variable value; a physical-value calculation unit that performs a multiplication of the second derivative calculated by the second-derivative calculation unit by a square value of the analytical small variable value, a multiplication of the first derivative calculated by the first-derivative calculation unit by the analytical small variable value, and an addition of values obtained by the multiplications and a physical value of the second-order differential equation at the reference variable value, and outputs a value obtained by the addition as a physical value after an increment of the analytical small variable value; and an operation control unit that repeats a control of setting a variable value obtained by adding the analytical small variable value to the reference variable value as a new reference variable value and causing the second-derivative calculation unit, the first-derivative calculation unit, and the physical-value calculation unit to calculate the second derivative, the first derivative, and the physical value, using the new reference variable value, respectively.
Moreover, in the numerical analysis device according to the present invention, the analytical small variable value and the reference variable value are an analytical small time and reference time, respectively, and the first derivative and the second derivative are a first-order time derivative and a second-order time derivative, respectively.
Furthermore, in the numerical analysis device according to the present invention, the setting unit sets a termination condition of the numerical analysis, and when the termination condition of the numerical analysis is satisfied, the operation control unit performs a control of terminating the numerical analysis.
Moreover, the numerical analysis device according to the present invention further includes an output processing unit that outputs the physical value of each reference variable value calculated by the physical-value calculation unit, in association with the increment of the analytical small variable value.
Furthermore, in the numerical analysis device according to the present invention, the output processing unit generates a graph of the physical value for each reference variable value calculated by the physical-value calculation unit corresponding to the increment of the analytical small variable value.
Moreover, in the numerical analysis device according to the present invention, the second-order differential equation includes a variable value term.
Furthermore, in the numerical analysis device according to the present invention, the second-order differential equation includes a time term.
Moreover, in the numerical analysis device according to the present invention, the second-order differential equation is a simultaneous equation with multiple unknowns.
Furthermore, in the numerical analysis device according to the present invention, the numerical analysis is performed using a spreadsheet program.
Moreover, a numerical analysis program according to the present invention is for a numerical analysis of second-order differential equation for a target object to be analyzed. The numerical analysis program causes a computer to execute a setting procedure of performing a setting of parameters including an initial condition of the second-order differential equation and a setting of an analytical small variable value for the numerical analysis; a second-derivative calculating procedure of calculating a second derivative of the second-order differential equation at a reference variable value that is a variable at the start of the numerical analysis; a first-derivative calculating procedure of performing a multiplication of the second derivative calculated at the second-derivative calculating procedure by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputting a value obtained by the addition as a first derivative after an increment of the analytical small variable value; a physical-value calculating procedure of performing a multiplication of the first derivative calculated at the first-derivative calculating procedure by the analytical small variable value and an addition of a physical value of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputting a value obtained by the addition as a physical value after an increment of the analytical small variable value; and an operation control procedure of repeating a control of setting a variable value obtained by adding the analytical small variable value to the reference variable value as a new reference variable value and causing the second-derivative calculating procedure, the first-derivative calculating procedure, and the physical-value calculating-procedure to calculate the second derivative, the first derivative, and the physical value, using the new reference variable value, respectively.
Furthermore, a numerical analysis program according to the present invention is for a numerical analysis of second-order differential equation for a target object to be analyzed. The numerical analysis program causes a computer to execute a setting procedure of performing a setting of parameters including an initial condition of the second-order differential equation and a setting of an analytical small variable value for the numerical analysis; a second-derivative calculating procedure of calculating a second derivative of the second-order differential equation at a reference variable value that is a variable at the start of the numerical analysis; a first-derivative calculating procedure of performing a multiplication of the second derivative calculated at the second-derivative calculating procedure by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputting a value obtained by the addition as a first derivative after an increment of the analytical small variable value; a physical-value calculating procedure of performing a multiplication of the second derivative calculated at the second-derivative calculating procedure by a square value of the analytical small variable value, a multiplication of the first derivative calculated at the first-derivative calculating procedure by the analytical small variable value, and an addition of values obtained by the multiplications and a physical value of the second-order differential equation at the reference variable value, and outputting a value obtained by the addition as a physical value after an increment of the analytical small variable value; and an operation control procedure of repeating a control of setting a variable value obtained by adding the analytical small variable value to the reference variable value as a new reference variable value and causing the second-derivative calculating procedure, the first-derivative calculating procedure, and the physical-value calculating procedure to calculate the second derivative, the first derivative, and the physical value, using the new reference variable value, respectively.
Moreover, in the numerical analysis program according to the present invention, the analytical small variable value and the reference variable value are an analytical small time and reference time, respectively, and the first derivative and the second derivative are a first-order time derivative and a second-order time derivative, respectively.
Furthermore, in the numerical analysis program according to the present invention, the setting procedure includes setting a termination condition of the numerical analysis, and when the termination condition of the numerical analysis is satisfied, the operation control procedure includes performing a control of terminating the numerical analysis.
Moreover, the numerical analysis program according to the present invention further causes the computer to execute an output processing procedure of outputting the physical value of each reference variable value calculated at the physical-value calculating procedure, in association with the increment of the analytical small variable value.
Furthermore, in the numerical analysis program according to the present invention, the output processing procedure includes generating a graph of the physical value for each reference variable value calculated at the physical-value calculating procedure corresponding to the increment of the analytical small variable value.
Moreover, in the numerical analysis program according to the present invention, the second-order differential equation includes a variable value term.
Furthermore, in the numerical analysis program according to the present invention, the second-order differential equation includes a time term.
Moreover, in the numerical analysis program according to the present invention, the second-order differential equation is a simultaneous equation with multiple unknowns.
Furthermore, in the numerical analysis program according to the present invention, the numerical analysis is performed using a spreadsheet program.

EFFECT OF THE INVENTION

By using a numerical analysis device or a numerical analysis program according to the present invention, the effect is obtained that a numerical analysis process can simply be performed within a short time period without causing a large error in the analysis result, because a setting unit sets parameters, including initial conditions for a second-order differential equation, and analytical small variable value of the numerical analysis, a second-derivative calculation unit calculates a second derivative of a second-order differential equation at a reference variable value that is the variable value at the start of the analysis, a first-derivative calculation unit multiplies the second derivative calculated by the second-derivative calculation unit by the analytical small variable value, adds the first derivative of the second-order differential equation at the reference variable value to this multiplied value, and equates this added value to the first derivative after increment of the analytical small variable value, a physical-value calculation unit multiplies the first derivative calculated by the first-derivative calculation unit by the analytical small variable value, adds the physical value of the second-order differential equation at the reference variable value to this multiplied value, and equates this added value to the physical value after increment of the analytical small variable value, an operation control unit repeats control to set the variable value that is addition of the reference variable value with the analytical small variable value to the new reference variable value, and to calculate the second derivative, the first derivative, and the physical value by the second-derivative calculation unit, the first-derivative calculation unit, and the physical-value calculation unit, by using the new reference variable value, respectively.
Also by using a numerical analysis device or a numerical analysis program according to the present invention, the effect is obtained that a numerical analysis process can simply be performed within a short time period without causing a large error in the analysis result, because a setting unit sets parameters, including initial conditions for a second-order differential equation, and analytical small variable value of a numerical analysis, a second-derivative calculation unit calculates a second derivative of a second-order differential equation at a reference variable value that is the variable value at the start of the analysis, a first-derivative calculation unit multiplies the second derivative calculated by the second-derivative calculation unit by the analytical small variable value, adds the first derivative of the second-order differential equation at the reference variable value to this multiplied value, and equates this added value to the first derivative after increment of the analytical small variable value, a physical-value calculation unit multiplies the second derivative calculated by the second-derivative calculation unit by the square value of the analytical small variable value, multiplies the first derivative calculated by the first-derivative calculation unit by the analytical small variable value, adds the physical value of the second-order differential equation at the reference variable value to these multiplied values, and equates this added value to the physical value after increment of the analytical small variable value, the operation control unit repeats control to set the variable value that is addition of the reference variable value with the analytical small variable value to the new reference variable value, and to calculate the second derivative, the first derivative and the physical value by the second-derivative calculation unit, the first-derivative calculation unit, and the physical-value calculation unit, by using the new reference variable value, respectively.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a block diagram of an outline of configuration of a dynamic analysis device that is a numerical analysis device according to the first embodiment of the present invention.

FIG. 2 is a flow chart of a dynamic analysis procedure by the dynamic analysis device shown in FIG. 1.

FIG. 3 is a diagram of the concept of the dynamic analysis process according to the first embodiment of the present invention.

FIG. 4 is a diagram of the concept of the dynamic analysis process according to conventional Euler's method.

FIG. 5 is a diagram of the concept of the dynamic analysis process according to conventional modified Euler's method.

FIG. 6 is a schematic diagram illustrating the concept of configuration of a spring-mass system to which analysis equation is applied.

FIG. 7 is a diagram of a dynamic analysis result of the spring-mass system shown in FIG. 6.

FIG. 8 is a view of an example of a display screen of a spreadsheet program for dynamic analysis of the spring-mass system shown in FIG. 7.

FIG. 9 is a chart of an example of a dynamic analysis for the spring-mass system shown in FIG. 6 performed by a spreadsheet program.

FIG. 10 is a graph of the result shown in FIG. 9.

FIG. 11 is a chart of an example of a dynamic analysis result for the same spring-mass system by Euler's method.

FIG. 12 is a graph of the result shown in FIG. 11.

FIG. 13 is a graph of an example of a dynamic analysis result for the same spring-mass system by modified Euler's method.

FIG. 14 is a block diagram of an outline of configuration of a dynamic analysis device according to the second embodiment of the present invention.

FIG. 15 is a flow chart of a dynamic analysis procedure by the dynamic analysis device shown in FIG. 14.

FIG. 16 is a diagram of the concept of the dynamic analysis process according to the second embodiment of the present invention.

FIG. 17 is a schematic diagram of an outline of configuration of a spring-mass dashpot system to which dynamic analysis is applied.

FIG. 18 is a graph of a dynamic analysis result of the mass-dashpot system shown in FIG. 17.

FIG. 19 is a schematic diagram of an outline of configuration of a transporting machine to which dynamic analysis is applied, accompanied by a simple harmonic motion and a rotational motion.

FIG. 20 is a graph of a dynamic analysis result of the displacement of the transporting machine shown in FIG. 19.

FIG. 21 is a graph of a dynamic analysis result of the angular displacement of the transporting machine shown in FIG. 19.

FIG. 22 is a schematic diagram of an outline of configuration of multiple bodies of revolution, to which dynamic analysis is applied, connected by a power transmission such as a torsion axis, joint and gear.

FIG. 23 is a graph of a dynamic analysis result of the angular displacement of each body of revolution connected by a power transmission shown in FIG. 22.

FIG. 24 is a graph of a dynamic analysis result of the angular displacement of each body of revolution connected by a power transmission shown in FIG. 22.

FIG. 25 is a schematic diagram of an outline of configuration of a slider for a speed governor and a turbine, to which dynamic analysis is applied.

FIG. 26 is a graph of a dynamic analysis result of the displacement of the slider of the speed governor shown in FIG. 25.

FIG. 27 is a graph of dynamic analysis result of the angular displacement of the turbine shown in FIG. 25.

FIG. 28 is a schematic diagram of an outline of configuration of a cutting tool of a turning machine such as a lathe, to which dynamic analysis is applied.

FIG. 29 is a graph of dynamic analysis result of the displacement of a cutting tool of a turning machine such as a lathe shown in FIG. 28.

FIG. 30 is a schematic diagram of an outline of configuration of a reciprocating machine, to which dynamic analysis is applied, having a crank.

FIG. 31 is a graph of a dynamic analysis result of the reciprocating machine having a crank shown in FIG. 30.

FIG. 32 is a graph of a dynamic analysis result of the piston of the reciprocating machine having a crank shown in FIG. 30.

FIG. 33 is a schematic diagram of an outline of configuration of a brake disc to which dynamic analysis is applied.

FIG. 34 is a graph of dynamic analysis result of the brake disc shown in FIG. 33.

FIG. 35 is a schematic diagram of an outline of configuration of a manipulator-like device, to which dynamic analysis is applied, connected by a rigid body link.

FIG. 36 is a graph of dynamic analysis result of the manipulator-like device connected by a rigid body link shown in FIG. 35.

FIG. 37 is a graph of a dynamic analysis result of the manipulator-like device connected by a rigid body link shown in FIG. 35.

FIG. 38 is a schematic diagram of an outline of configuration of particles or celestial bodies, to which dynamic analysis is applied, that influence forces like gravity each other.

FIG. 39 is a graph of a dynamic analysis result of the displacement in x direction of the Earth shown in FIG. 38.

FIG. 40 is a diagram of a coordinate system and configuration relationship to analyze three-body problem of the Sun, the Earth and the Moon.

FIG. 41 is a graph of a dynamic analysis result of the orbit of the Moon in three-body problem placing the Earth in the center.

FIG. 42 is a graph of a dynamic analysis result of the orbit of the Moon by conventional Euler's method placing the Earth in the center.

FIG. 43 is a schematic diagram of an outline of configuration of a surge tank, to which dynamic analysis is applied, with a centrifugal pump.

FIG. 44 is a graph of a dynamic analysis result of the change of water level of the surge tank shown in FIG. 43.

FIG. 45 is a schematic diagram of an outline of configuration of an LCR circuit to which dynamic analysis is applied.

FIG. 46 is a graph of a dynamic analysis result of the change of electric charge in the LCR circuit shown in FIG. 45.

FIG. 47 is a system configuration schematic diagram of a configuration of a computer system that implements the dynamic analysis device.

FIG. 48 is a block diagram of a configuration of a main unit of the computer system shown in FIG. 47.

EXPLANATIONS OF LETTERS OR NUMERALS

- 1 DYNAMIC ANALYSIS DEVICE
- 10 INPUT UNIT
- 11 OUTPUT UNIT
- 12 OPERATION UNIT
- 13 OUTPUT PROCESSING UNIT
- 14 MEMORY UNIT
- 20 OPERATION CONTROL MODULE
- 21 INCLUSION MODULE
- 22 SECOND-ORDER-DERIVATIVE CALCULATION MODULE
- 23 FIRST-ORDER-DERIVATIVE CALCULATION MODULE
- 24,24A PHYSICAL-VALUE CALCULATION MODULE

BEST MODE(S) FOR CARRYING OUT THE INVENTION

Exemplary embodiments of a dynamic analysis device as an example of a numerical analysis device and a dynamic analysis program as an example of a numerical analysis program according to the present invention will be explained in detail.

First Embodiment

FIG. 1 is a block diagram of an outline of configuration of a dynamic analysis device that is an example of numerical analysis device according to the first embodiment of the present invention.
As shown in FIG. 1 dynamic analysis device 1 includes an input unit 10 that is implemented with a keyboard and a pointing device to input various information, an output unit 11 that is implemented with an LCD display, a printer and the like to output processing result, a memory unit 14 that is used to store programs and data as well as to store calculation results temporarily, and a control unit C that controls performing a dynamic analysis process.
The control unit C includes an operation unit 12 and an output processing unit 13, the operation unit 12 performs a dynamic analysis processing, and the output processing unit 13 performs a process to output the result of dynamic analysis to the output unit 11. The operation unit 12 includes an operation control module 20, an inclusion module 21, a second-order-derivative calculation module 22, a first-order-derivative calculation module 23, and a physical-value calculation module 24.
Referring the flowchart shown in FIG. 2, the procedure of control process of the control unit C is explained. Firstly, in FIGS. 1 and 2, analysis information is input by the inclusion module 21 by using the input unit 10 and the output unit 11 (step S101). The analysis information includes an analysis equation α(t)=f(x(t),v(t)), an analytical small time increment Δt, further, parameters, initial conditions and termination conditions and the like. The analysis equation here is an ordinary second-order differential equation and the α(t) is a second-order time derivative such as acceleration, x(t) is a physical value such as position, and v(t) is a first-order time derivative such as velocity.
The second-order-derivative calculation module 22 calculates an acceleration α(t) at the time of reference t using initial conditions (step S102). Then the first-order-derivative calculation module 23 multiplies the calculated acceleration α(t) by the analytical small time increment Δt, adds the velocity at the reference time to this multiplied value and equates this added value to the velocity v(t+Δt) after elapse of the analytical small time increment Δt (step S103). Namely, v(t+Δt)=v(t)+α(t)·Δt is calculated.
Then the physical-value calculation module 24 multiplies the calculated velocity v(t+Δt) by the analytical small time increment Δt, adds the displacement x(t) at the reference time t to this added value, and equates this added value to the displacement x(t+Δt) after elapse of the analytical small time increment Δt. Namely, x(t+Δt)=x(t)+v(t+Δt)·Δt is calculated.
The operation unit 12 stores the calculated value at step S102 to S104 into the memory unit 14 and also performs an output processing to enable output from the output unit 11 via the output processing unit 13 (step S105). Further the operation unit 12 decides whether t exceeds a specified value (step S106). Here the specified value is time elapsed from the reference time t. This elapsed time can be given as Δt×n, where n is the number of repetition thereafter. If t does not exceed the specified value (step S106, No), value (t+Δt) that is the addition of the reference time t with the analytical small time increment Δt is set as a new reference time, the step moves to step S102 (step S107), and processes from step S102 is repeated. If t exceeds the specified value (step S106, Yes), the process is terminated.
In other words, a displacement x(t+Δt) after lapse of the analytical small time increment Δt is calculated in this dynamic analysis process by the dynamic analysis device 1 as follows: firstly the acceleration α(t) that is an analysis equation at the reference time t is calculated by equation (1) using the displacement x(t) that is a physical value at the reference time t and the velocity v(t) that is a first-order time derivative at the reference time t. Then v(t+Δt) that is the velocity after elapse of the analytical small time increment Δt is calculated by equation (2), further x(t+Δt) that is the displacement after elapse of the analytical small time increment Δt is calculated by equation (3).
α(t)=f(x(t),v(t)) (1)
v(t+Δt)=v(t)=+α(t)·Δt (2)
x(t+Δt)=x(t)+v(t+Δt)·Δt (3)
Some conventional dynamic analysis devices utilize Euler's method where displacements are calculated by equation (4) below where v(t) is used instead of v(t+Δt) in the equation (3). Namely, displacements are calculated by
x(t+Δt)=x(t)+v(t)·Δt (4)
However, there are some cases with Euler's method where divergence occurs when a dynamic analysis is performed on a simple harmonic motion and the like. Then modified Euler's method was contrived that uses equation (5) below instead of equation (3).
x(t+Δt)=x(t)+(1/2·(v(t)+v(t+Δt)))·Δt (5)
In modified Euler's method the dynamic analysis is performed replacing v(t+Δt) with (1/2·(v(t)+v(t+Δt)) that is the average velocity between the reference time t and (t+Δt) which is the time after elapse of the analytical small time increment Δt from the reference time t. By modified Euler's method degree of divergence can be reduced but it may still result in a dynamic analysis with divergence and low accuracy. Therefore, when a dynamic analysis by Euler's method or modified Euler's method is applied it has been inevitable to make the analytical small time increment Δt extremely small and to consume a lot of time to obtain necessary accuracy. In contrast to this, by the dynamic analysis process where equations (1) to (3) are repeated, a result of high accuracy can be obtained even if the analytical small time increment Δt is very small.
FIGS. 3 to 5 are diagrams of the concepts of dynamic analysis processes according to the present invention, conventional Euler's method, and modified Euler's method respectively. FIG. 3 is a diagram of the concept of a dynamic analysis process according to the first embodiment of the present invention, FIG. 4 is a diagram of the concept of the dynamic analysis process according to conventional Euler's method, and FIG. 5 is a diagram of the concept of the dynamic analysis process according to modified Euler's method. Referring FIGS. 3 to 5, in the dynamic analysis process according to the present invention, acceleration is obtained from the displacement and the velocity at time t, velocity at time (t+Δt) that is time after elapse of the analytical small time increment Δt is obtained from the velocity and the acceleration at time t as shown in FIG. 3. Then displacement at time (t+Δt) is obtained from the velocity at time (t+Δt) and the displacement at time t, and it is specifically characterized in that displacement at time (t+Δt) is obtained by using the velocity at time (t+Δt). On the other hand, in Euler's method displacement at time (t+Δt) is obtained by using the velocity at time t and the displacement at time t as shown in FIG. 4. Further in conventional modified Euler's method displacement at time (t+Δt) is obtained by using velocity (1/2·(v(t)+v(t+Δt))), which is the average velocity at time t and time (t+Δt) as shown in FIG. 5.
A result of processing an analysis equation for a concrete simple harmonic motion is now explained. FIG. 6 is a schematic diagram of a conceptual structure of a spring-mass system to which an analysis equation is applied. A simple harmonic motion appears in the spring-mass system shown in FIG. 6. An equation of motion for the spring-mass system shown in FIG. 6 is given as equation (6) below.
m·α(t)=−k·x(t) (6)
In this spring-mass system a spring constant k=4π²(N/m), mass m=1.0 (kg), initial displacement x(0)=2.0 (m), and initial velocity v(0)=0.00 (m) are assumed.
An analysis equation for this equation of motion is given as equation (7) below.
α(t)=−k/m·x(t) (7)
After the parameters are assigned the analysis equation becomes equation (8) below.
α(t)=−4π² ·x(t) (8)
The algebraic solution of this spring-mass system is expressed in equation (9) below,
x(t)=cos(2π·t) (9)
where the analytical small time increment Δt=0.01(s) is applied.
FIG. 7 is a diagram of a dynamic analysis result of the spring-mass system shown in FIG. 6. This dynamic analysis is performed with the analytical small time increment Δt=0.01(s). As shown in FIG. 7, divergence occurs in Euler's method as well as in modified Euler's method with time lapse, while a result near the algebraic solution without divergence is obtained from the dynamic analytical process according to the present invention. Therefore, according to the first embodiment of the present invention can perform a dynamic analysis of high accuracy in a short period of time comparing with Euler's method or modified Euler's method. In other words, the dynamic analysis process according to the first embodiment of the present invention can perform a dynamic analysis of much higher accuracy than Euler's method or modified Euler's method in the same analysis time period.
The dynamic analysis process of the first embodiment is not a complicate program and therefore it can easily be implemented by software using a spreadsheet program. FIG. 8 is a view of an example of display screen of a spreadsheet program for a dynamic analysis of the spring-mass system shown in FIG. 7. In FIG. 8 initial time t=0.00 is firstly set in cell C11 of the time column and then the analytical small time increment Δt=0.01 is added to C21 to Cn1. The initial displacement x(0)=2.00 is set in cell C12 of the displacement column of the beginning of the interval that represents reference time point t, and the initial velocity v(0)=0.00 is set in cell C13 of the velocity column of the beginning of the interval. Acceleration is calculated using the displacement and the velocity given in cell C12 and C13 respectively and the result is shown in cell C14 of the acceleration column. The calculation result corresponding to equation (2) is shown in cell C15 of the velocity column of the end of the interval that represents after lapse of the analytical small time increment Δt from the reference time t, and the calculation result corresponding to equation (2) is shown in cell C16 of the velocity column of the end of the interval. The results indicated in cells C15 and C16 are copied into cells C22 and C23 of the next row for the time after lapse of the analytical small time increment Δt as the displacement and velocity at the beginning of the next interval respectively. Further, acceleration is calculated using the displacement and velocity indicated in C22 and C23, the result of the calculation is shown in cell C24, and then the velocity and displacement after next lapse of the analytical small time increment Δt are calculated using the displacement, velocity, and acceleration shown in cell C22, C23 and C24 respectively, and the results are shown in cells C25 and C26. In the same way hereafter, the process, in which the velocity and displacement at the end of the interval are utilized as the velocity and displacement at the beginning of the next interval after lapse of the analytical small time increment Δt and acceleration at this time point is calculated and displayed, and also the velocity and displacement at the end of the next interval after lapse of the analytical small time increment Δt are calculated and displayed, is repeated for repeat count that is set in the time row, namely number of the rows. These results can also be displayed in a graph by graphing process. The graphing process is carried out by the output processing unit 13.
FIG. 9 is a chart of an example that a dynamic analysis for the spring-mass system shown in FIG. 6 is performed by a spreadsheet program. Here an initial displacement x(0)=1.000 (m) is assumed. FIG. 10 is a graph of the result shown in FIG. 9, where the horizontal axis indicates time and the vertical axis indicates displacement. As shown in FIG. 10, the algebraic solution and the result of a numerical analysis (dynamic analysis) are nearly in coincidence and the error is ±0.042 m. FIG. 10 shows the result of a dynamic analysis for up to 10(s), but even for 100(s) the error of the dynamic analysis is only ±0.134 m.
FIG. 11 is a chart of an example that a dynamic analysis for the same spring-mass system is performed by Euler's method, FIG. 12 is a graph of the result shown in FIG. 11, and FIG. 13 is a graph of an example that a dynamic analysis for the same spring-mass system is performed by modified Euler's method. As shown in FIGS. 11 and 12, when a dynamic analysis is performed by Euler's method using the same analytical small time increment, the error between the dynamic analysis and the algebraic solution increases with time and diverges. In this case, the error increases up to ±6 m in 10(s). When modified Euler's method is employed, as shown in FIG. 13, the error between the dynamic analysis and the algebraic solution also increases with time and diverges. In this case, the error is smaller than that by Euler's method, but still reaches about ±2 m in 10(s).

Second Embodiment

Now the second embodiment is explained. In the first embodiment, using the present first-order time derivative like velocity and the present second derivative like acceleration, the first-order time derivative at the next time point is calculated, and then using the first-order time derivative at the next time point and the present physical value like displacement the next physical value at the next time point is obtained, in this second embodiment, however, the physical value at the next time point is obtained at once using the present first-order time derivative, the present second-order time derivative, and the present physical value. In this regard, the first-order time derivative at the next time point that is used to calculate the second-order time derivative at the next time point shall already be calculated at present time point.
FIG. 14 is a block diagram of an outline of configuration of a dynamic analysis device according to the second embodiment of the present invention. As shown in FIG. 14, the dynamic analysis device includes a physical-value calculation module 24 a instead of the physical-value calculation module 24. Otherwise, the configuration is the same as in the first embodiment, and therefore the same reference numerals are applied to the same components in the configuration.
The second-order-derivative calculation module 22 calculates acceleration α(t) using velocity v(t) that was calculated by the first-order-derivative calculation module 23 and displacement x(t) at the reference time t that was calculated by the physical-value calculation module 24, and then the physical-value calculation module 24 a calculates directly displacement x(t+Δt) after lapse of the analytical small time increment Δt using the acceleration α(t), the velocity v(t), and the displacement x(t) at the reference time t. It means that the displacement x(t+Δt) can be directly obtained from equation (10) below that is derived by combining equation (2) with equation (3).
x(t+Δt)=x(t)+v(t)·Δt+α(t)·Δt ² (10)
FIG. 15 is a flow chart of a dynamic analysis procedure by the dynamic analysis device shown in FIG. 14. In the dynamic analysis procedure shown in FIG. 15 process step S204 is applied instead of process step S104 in the first embodiment. Otherwise steps S201 to S203 are same as S101 to S103, and S205 to S207 are same as S105 to S107, respectively.
Namely at step S204 the acceleration α(t) calculated at step S204 is multiplied by square value of the analytical small time increment Δt, the velocity v(t) at the reference time t is multiplied by the analytical small time increment Δt, each multiplied value is added to the displacement x(t) at the reference time t, and this added value is equated to the displacement x(t+Δt) after lapse of the analytical small time increment Δt. The reason why v(t) at the reference time is calculated at step S203 is that, different from the first embodiment, the v(t) necessary for calculation of the displacement at the next time point has not yet been calculated.
Thus in the second embodiment, different from the first embodiment, the velocity v(t+Δt) after lapse of the analytical small time increment Δt can be calculated at once using the present velocity, acceleration, and displacement.
The concept of dynamic analysis process in the second embodiment according to the present invention is shown in FIG. 16. FIG. 16 is a diagram of the concept of the dynamic analysis process according to the second embodiment of the present invention. Referring to FIG. 16, in the dynamic analysis process in the second embodiment according to the present invention, acceleration is calculated from the displacement and the velocity at time t, and using the velocity and the acceleration at time t displacement at time (t+Δt) after lapse of the analytical small time increment Δt is calculated.
In the first and the second embodiment analysis equation α(t)=f(x(t), v(t)) is applied for explanation, but it can also be α(t)=f(x(t)) without velocity term v(t). Also, analysis equation a (t) may include time term. Namely α(t) can be α(t)=f(x(t), v(t), t) or α(t)=f(x(t), t). Further, depending on the object of analysis, a plurality of analysis equations can be simultaneous equations. In other words, the analysis equations can be equations with multiple unknowns. By generalization these equations can be expressed as simultaneous differential equations below, αi(t)=f(x1(t), . . . , xn(t), v1(t), . . . , vn(t), t) where i=1 to n.
In the first and the second embodiments the dynamic analysis device is used as an example of the numerical analysis device for explanation, however, it can be applied not only to the dynamic analysis where time is a variable but also can be applied to a general numerical analysis. Namely the numerical analysis can be performed where the variable is not time t but temperature T or distance L. For example, a numerical analysis with temperature change can be performed when temperature T and analytical small temperature increment ΔT are applied instead of time t and analytical small time increment Δt.
Following are explained as examples of the first and the second embodiment: a dynamic analysis process of a spring-mass dashpot system, a dynamic analysis process of a transporting machine, a dynamic analysis process of a plurality of bodies of revolution connected by power transmission and the like, a dynamic analysis process of a speed governor, a dynamic analysis process of a cutting tool of a turning machine and the like, a dynamic analysis process of a reciprocating machine having a crank, a dynamic analysis process of a brake utilizing friction, a dynamic analysis process of a manipulator-like device connected with a rigid body link, a dynamic analysis process of particles or celestial bodies influencing force like gravity each other, a dynamic analysis process of a three-body problem of the Sun, the Earth and the Moon, a dynamic analysis process of surging characteristics of a system having a centrifugal pump, and a dynamic analysis process of current or electric charge generated in an electric or electronic circuit having inductance, capacitance, and resistance.
(First Application)
FIG. 17 is a schematic diagram of an outline of configuration of a spring-mass dashpot system to which dynamic analysis is applied. The equation of motion specific for the spring-mass dashpot system is given as an equation (11) below.
mα+cv+kx=f·sin((2π/T)·t) (11)
where,
m: mass of the mass point (1500.00),
c: viscous damping coefficient of the dashpot (4242.64),
k: spring coefficient (12000),
f: amplitude of the external force (1000.0),
T: period of the external force (1.00),
t: reference time,
x: displacement of the mass point,
v: velocity of the mass point,
α: acceleration of the mass point,
t₀: beginning time of analysis,
x₀: initial displacement of the mass point (0.10),
v₀: initial velocity of the mass point (1.00),
α₀: initial acceleration of the mass point, and these are shown in the International System of Units.
In this case, the analysis equation is derived from the equation (11) and given as an equation (12) below.
α(t)=(l/m)·(f·sin((2π/T)·t)−cv(t)−kx(t)) (12)
α(t₀) is calculated by assigning parameters, t=t₀, x(t₀)=x₀, and v(t₀)=v₀to the α(t) given by the equation (12), and v(t+Δt) and x(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equations (2) and (3). Acceleration α(t+Δt) after lapse of the analytical small time increment Δt using v(t+Δt) and x(t+Δt), and the process is repeated specified number of times.
FIG. 18 a graph of a dynamic analysis result of the mass dashpot system shown in FIG. 17, and the horizontal axis indicates time and the vertical axis indicates displacement. From this dynamic analysis result, as shown in FIG. 18, almost the same dynamic characteristic has been obtained as the output result of the dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including a feedback processing.
(Second Application)
FIG. 19 is a schematic diagram of an outline of configuration of a transporting machine, to which the dynamic analysis is applied, accompanied by a simple harmonic motion and rotational motion. The equations of motion specific for this transporting machine are given as equations (13) and (14) below.
mα+k ₁(x+L ₁θ)+k ₂(x−L ₂θ)=0 (13)
J(d ² θ/dt ²)+L ₁ k ₁(x+L ₁θ)−L ₂ k ₂(x−L ₂θ) (14)
where,
m: mass of the transporting machine (1500.00),
J: moment of inertia of the transporting machine (3125),
k₁: first spring constant (3000),
k₂: second spring constant (3000),
L₁: first distance from the center of gravity (2.00),
L₂: second distance from the center of gravity (2.00),
x: displacement of the transporting machine (initial value x₀=0.10),
v: velocity of the transporting machine (initial value v₀=0.00),
α: acceleration of the transporting machine,
θ: angular displacement of the transporting machine (initial value θ₀=0.26),
ω: angular velocity of the transporting machine (initial value ω₀=0.00),
d²θ/dt²: angular acceleration of the transporting machine, and these are shown in the International System of Units.
In this case, analysis equations are derived from the equations (13) and (14), and are given as equations (15) and (16) below.
α(t)=−(k ₁ /m)(x(t)+L ₁θ(t))−(k ₂ /m)(x(t)−L ₂θ(t)) (15)
(d ² θ/dt ²)=−(L ₁ k ₁ /J)(x(t)+L ₁θ(t))+(L ₂ k ₂ /J)(x(t)−L ₂θ(t)) (16)
Then α(t₀) and d²θ/dt²(t₀) are calculated by assigning parameters, initial displacement v(t₀)=0.00 and initial angular displacement θ₀(t₀)=0.25 to the equations (15) and (16) respectively, further x(t+Δt) and θ(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by equation (10). Acceleration α(t+Δt) and (d²θ(t+Δt)/dt²) after the next lapse of the analytical small time increment Δt are calculated using x(t+Δt) and θ(t+Δt), and the process is repeated specified number of times. When the equations (2) and (3) are utilized as in the first embodiment, velocity of the transporting machine v(t+Δt) and angular velocity of the transporting machine ω(t+Δt) are calculated by the equation (3). In this case initial velocity of the transporting machine v(t₀)=0.00 and angular velocity of the transporting machine is ω(t₀)=0.00 are assumed.
FIG. 20 is a graph of a dynamic analysis result of displacement of the transporting machine shown in FIG. 19, and FIG. 21 is a graph of a dynamic analysis result of angular displacement of the transporting machine shown in FIG. 19. From each dynamic analysis result, as shown in FIGS. 20 and 21, almost the same dynamic characteristics have been obtained as the output results of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Third Application)
FIG. 22 is a schematic diagram of an outline of configuration of multiple bodies of revolution to which dynamic analysis is applied, connected by a power transmission such as a torsion axis, joint, and gear. The equations of motion for the object of analysis are given as equations (17) and (18) below.
I ₁(d ²θ₁ /dt ²)+k ₁₂(θ₁−θ₂)=0 (17)
I ₂(d ²θ₂ /dt ²)+k ₁₂(θ₂−θ₁)=0 (18)
where
I₁: moment of inertia of the first body of revolution (9500),
I₂: moment of inertia of the second body of revolution (53900),
k₁₂: torsion spring constant (98980000),
θ₁: angular displacement of the first body of revolution (initial value=0.000),
θ₂: angular displacement of the second body of revolution (initial value=0.026),
ω₁: angular velocity of the first body of revolution (initial value=0.000),
ω₂: angular velocity of the second body of revolution (initial value=0.000),
and these are shown in the International System of Units.
In this case, analysis equations are derived from the equations (17) and (18), and are given as equations (19) and (20) below.
(d ²θ₁ /dt ²)=−(k ₁₂ /I ₁)·(θ₁−θ₂) (19)
(d ²θ₂ /dt ²)=−(k ₁₂ /I ₂)·(θ₂−θ₁) (20)
(d²θ₁/dt²) and (d²θ₂/dt²) are respectively calculated assigning parameters, initial angular displacement θ₁₀=0.000, initial angular displacement θ₂₀=0.026 to the equations (19) and (20), then θ₁(t+Δt) and θ₂(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equation (10). (d²θ₁(t+Δt)/dt²) and (d²θ₂(t+Δt)/dt²) after lapse of the next analytical small time increment Δt are calculated using θ₁(t+Δt) and θ₂(t+Δt), and the process is repeated specified number of times. When the equations (2) and (3) are utilized as in the first embodiment, angular velocity ω₁(t+Δt) and ω₂(t+Δt) are respectively calculated by the equation (3). Here initial angular velocity ω₁(t+Δt)=0.000 and ω₂(t+Δt)=0.000 are assumed.
FIG. 23 and FIG. 24 are graphs of dynamic analysis results of the angular displacement of each body of revolution connected by a power transmission shown in FIG. 22. From each dynamic analysis result, as shown in FIGS. 23 and 24, almost the same dynamic characteristics have been obtained as the output results of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Fourth Application)
FIG. 25 is a schematic diagram of an outline of configuration of a slider for a speed governor and a turbine, to which dynamic analysis is applied. The equations of motion specific for this speed governor are given as equations (21) and (22) below.
m _gvα_gv +c _gv v _gv +k _gv x _gv =A _gvω_gv (21)
J _gv(d ²φ_gv /dt ²)=−B _gv x _gv (22)
where
m_gv: equivalent mass of the slider for the speed governor (10.00),
c_gv: viscous damping coefficient of the speed governor (50.66),
k_gv: spring constant of the speed governor (1579.14),
A_gv: upward force working on the slider for the speed governor by every 1 rad/s increment of rotation of the turbine (10.00),
J_gv: equivalent moment of inertia of the turbine (2.50),
B_gv: increase of torque by every unit length descent of the slider for the speed governor (1000.00),
x_gv: displacement of the slider for the speed governor (initial value=0.00),
v_gv: velocity of the slider for the speed governor (initial value=0.00),
α_gv: acceleration of the slider for the speed governor,
φ_gv: angular displacement of the turbine (initial value=0.00),
ω_gv: angular velocity of the turbine (initial value=0.20), and these are shown in the International System of Units.
In this case, analysis equations are derived from the equations (21) and (22), and are given as equations (23) and (24) below.
α_gv=(−c _gv v _gv −k _gv x _gv +A _gvω_gv)/m _gv (23)
(d ²φ_gv /dt ²)=−(B _gv /J _gv)·x _gv (24)
Acceleration α_gvand angular acceleration (d²φ_gv/dt²) are calculated assigning parameters, initial displacement x_gv(t0)=0.00, initial angular displacement φ_gv(t0)=0.00, initial velocity v_gv(t0)=0.00, initial angular velocity ω_gv(t0)=^0.20to the equations (23) and (24) respectively. Then velocity v_gv(t+Δt) and ω_gv(t+Δt) after lapse of the analytical small time increment Δt(=0.0005(s)) are calculated by the equation (2). Further displacement x_gv(t+Δt) and angular displacement φ_gv(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equation (3) respectively. Acceleration α_gv(t+Δt) and angular acceleration (d²φ_gv(t+Δt)/dt²) after lapse of the analytical small time increment Δt(=0.0005 (s)) are calculated using the obtained v_gv(t+Δt), ω_gv(t+Δt), x_gv(t+Δt), and φ_gv(t+Δt), and the process above is repeated specified number of times.
FIG. 26 is a graph of a dynamic analysis result of the displacement of the slider for the speed governor shown in FIG. 25, and FIG. 27 is a graph of a dynamic analysis result of the angular displacement of the turbine shown in FIG. 25. From each dynamic analysis result, as shown in FIGS. 26 and 27, almost the same dynamic characteristics have been obtained as the output results of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Fifth Application)
FIG. 28 is a schematic diagram of an outline of configuration of a cutting tool of a turning machine such as a lathe, to which the dynamic analysis is applied. The equation of motion specific for the cutting tool is given as an equation (25) below.
m _bα_b +c _b v _b +k _b x _b=μ(V−v _b)P _b (25)
where,
m_b: mass of the cutting tool (10.00),
c_b: viscous damping coefficient (20.00),
k_b: spring constant (1000.00),
V: cutting speed (2.00),
P_b: radial force (100.00),
x_b: displacement of the cutting tool (initial value=0.00),
v_b: velocity of the cutting tool (initial value=1.00),
α_b: acceleration of the cutting tool,
μ(V−v _b): 0.0075(V−v _b)²−0.055(V−v _b)+0.400,
and these are shown in the International System of Units.
In this case, analysis equation is derived from the equation (25) and is given as an equation (26) below.
α_b=(−c _b v _b −k _b x _b+μ(V−v _b)P _b)/m _b (26)
Acceleration α_bis calculated assigning parameters, initial displacement x_b(t₀)=0.00, initial velocity v_b(t₀)=1.00 to the equation (26). Then velocity v_b(t+Δt) after lapse of the analytical small time increment Δt(=0.02(s)) is calculated by the equation (2). Further displacement x_b(t+Δt) after lapse of the analytical small time increment Δt(=0.02(s)) is calculated by the equation (3). Acceleration α_b(t+Δt) after lapse of the analytical small time increment Δt(=0.02(s)) is calculated using these obtained v_b(t+Δt) and x_b(t+Δt), and the process is repeated specified number of times.
FIG. 29 is a graph of a dynamic analysis result of the displacement of the cutting tool of the turning machine such as a lathe shown in FIG. 28. From this dynamic analysis result, as shown in FIG. 29, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Sixth Application)
FIG. 30 is a schematic diagram of an outline of configuration of a reciprocating machine, to which the dynamic analysis is applied, with a crank. The equation of motion specific for the reciprocating machine with a crank is given as the equation (27) below.
(I _cr +m _pi R _cr ²(sin θ_cr+(ε/2)sin 2θ_cr)²)(d ²θ_cr /dt ²)=−m _pi R _cr ²(sin θ_cr+(ε/2)sin 2θ_cr)(cos θ_cr+εcos 2θ_cr)(dθ _cr /dt)²−(Tr _cr +P _pi R _cr(sin θ_cr+(ε/2)sin 2θ_cr)) (27)
where,
m_pi: reciprocating mass (15.00),
I_cr: moment of inertia of total rotating mass around the crankshaft (7.73),
R_cr: arm length of the crankshaft (0.13),

ε: R_cr/L_con,

L_con: length of the connection (0.25),
P_pi: external force (0.00),
Tr_cr: resistive torque (1.00)
θ_cr: angular displacement of the crank (initial value=0.26),
dθ_cr/dt: angular velocity of the crank (initial value=12.57),
and these are shown in the International System of Units.
In this case the analysis equation is derived from the equation (27) and is given as equation (28) below.
(d ²θ_cr /dt ²)=−{m _pi R _cr ²(sin θ_cr+(ε/2)sin 2θ_cr)(cos θ_cr+εcos 2θ_cr)(dθ _cr /dt)²}/(I _cr +m _pi R _cr ²(sin θ_cr+(ε/2)sin 2θ_cr)²)−{(Tr _cr +P _pi R _cr(sin θ_cr+(ε/2)sin 2θ_cr))}(I _cr +m _pi R _cr ²(sin θ_cr+(ε/2)sin 2θ_cr)²) (28)
Angular acceleration (d²θ_cr/dt²) is calculated assigning parameters, initial angular displacement θ_cr(t₀)=0.26, and initial angular velocity (dθ_cr/dt(t₀))=12.57 to the equation (28). Then angular velocity (dθ_cr(t+Δt)/dt) after lapse of the analytical small time increment Δt(=0.001(s)) is calculated by the equation (2). Further angular displacement θ_cr(t+Δt) after lapse of the analytical small time increment Δt(=0.001(s)) is calculated by the equation (3). Angular acceleration (d²θ_cr(t+Δt)/dt²) after lapse of the analytical small time increment Δt(=0.001(s)) is calculated using these obtained (dθ_cr(t+Δt)/dt) and θ_cr(t+Δt), and the process is repeated specified number of times.
FIG. 31 is a graph of a dynamic analysis result of the reciprocating machine with a crank shown in FIG. 30. From this dynamic analysis result, as shown in FIG. 31, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
The displacement of the piston can structurally be obtained from the crank angle and other parameters as shown in an equation (29) below.
x _pi=(R _cr +L _con)−(R _crcos θ_cr +L _con(1−ε²sin²θ_cr)^1/2) (29)
where,
x_pi: displacement of the piston,
R_cr: arm length of the crankshaft (0.13),

ε: R_cr/L_con,

L_con: length of the connection (0.25),
θ_cr: angular displacement of the crank,
and these are shown in the International System of Units. FIG. 32 is a graph of a dynamic analysis result of the piston of the reciprocating machine with the crank obtained from the equation (29). From this dynamic analysis result, as shown in FIG. 32, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Seventh Application)
FIG. 33 is a schematic diagram of an outline of configuration of a brake disc to which the dynamic analysis is applied. The equation of motion specific for the brake disc is given as an equation (30) below.
J _B(d ²θ_B /dt ²)=−u _B f _B (30)
where,
J_B: moment of inertia of the brake disc (3.00),
f_B: press force (500.00),
u_B: friction coefficient (0.30),
θ_B: angular displacement of the brake disc (initial value=0.00),
dθ_B(t)/dt: angular velocity of the brake disc (initial value=139.63),
and these are shown in the International System of Units.
In this case, the analysis equation is derived from the equation (30) and is given as an equation (31) below.
(d ²θ_B /dt ²)=−u _B f _B /J _B (31)
Angular acceleration (d²θ_B/dt²(t₀)) is calculated assigning parameters, initial angular displacement θ_B(t₀)=0.00, and initial angular velocity (dθ_B/dt(t₀))=139.63 to the equation (30). Then angular velocity (dθ_B(t+Δt)/dt) after lapse of the analytical small time increment Δt(=0.01(s)) is calculated by the equation (2). Further angular displacement. θ_B(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) is calculated by the equation (3). Angular acceleration (d²θ_B(t+Δt)/dt²) after lapse of the analytical small time increment Δt(=0.01(s)) is calculated using these obtained (dθ_B(t+Δt)/dt) and θ_B(t+Δt), and the process is repeated specified number of times.
FIG. 34 is a graph of a dynamic analysis result of the brake disc shown in FIG. 33. From this dynamic analysis result, as shown in FIG. 34, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Eighth Application)
FIG. 35 is a schematic diagram of an outline of configuration of a manipulator-like device, to which the dynamic analysis is applied, connected by a rigid body link. The equations of motion specific for the manipulator-like device connected by a rigid body link are given as equations (32) and (33) below.
(4/3)(d ²θ_link1 /dt ²)+(1/2)(d ²θ_link2 /dt ²)cos(θ_link1−θ_link2)−(1/2)(dθ _link2 /dt)²sin(θ_link2−θ_link1)+(3/2)(g/L _link)sin θ_link1=0 (32)
(1/3)(d ²θ_link2 /dt ²)+(1/2)(d ²θ_link1 /dt ²)cos(θ_link2−θ_link1)+(1/2)(dθ _link1 /dt)²sin(θ_link2−θ_link1)+(1/2)(g/L _link)sin θ_link2=0 (33)
where,
L_link: length of link 1 and link 2 (5.000),
g: acceleration of gravity (9.807),
θ_link1: angular displacement of link 1 (initial value=0.087),
dθ_link1/dt: angular velocity of link 1 (initial value=0.000),
θ_link2: angular displacement of link 2 (initial value=−0.087),
dθ_link2/dt:angular velocity of link 2 (initial value=0.000), and these are shown in the International System of Units.
In this case, analysis equations are derived from the equations (32) and (33) and are given as equations (34) and (35) below.
(d ²θ_link1 /dt ²)=(1/(4/3−(3/4)cos²(θ_link2−θ_link1)))·{(3/4)(dθ _link1 /dt)²cos(θ_link2−θ_link1)sin(θ_link2−θ_link1)+(1/2)(dθ _link2 /dt)²sin(θ_link2−θ_link1)+(3g/4L _link)(sin θ_link2cos(θ_link2−θ_link1)−2 sin θ_link1)} (34)
(d ²θ_link2 /dt ²)=(1/(1/3−(3/16)cos²(θ_link2−θ_link1)))·{−(1/2)(dθ _link1 /dt)²sin(θ_link2−θ_link1)−(3/16)(dθ _link2 /dt)²cos(θ_link2−Θ_link1)sin(θ_link2−θ_link1)+(9g/16L _link)(sin θ_link1cos(θ_link2−θ_link1)−(8/9)sin θ_link2)} (35)
Angular acceleration (d²θ_link1/dt²(t₀)) and (d²θ_link2/dt²(t₀)) are calculated assigning parameters, initial angular displacement θ_link1(t₀)=0.087, initial angular displacement θ_link2(t₀)=−0.087, initial angular velocity (dθ_link1/dt(t₀))=0.000, and initial angular velocity (dθ_link2/dt(t₀))=0.000 to the equations (34) and (35). Then angular velocity (dθ_link1(t+Δt)/dt) and (dθ_link2(t+Δt)/dt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equation (2). Further angular displacement θ_link1(t+Δt) and θ_link2(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equation (3). Angular acceleration (d²θ_link1(t+Δt)/dt²) and (d²θ_link2(t+Δt)/dt²) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated using these obtained (dθ_link1(t+Δt)/dt), (dθ_link2(t+Δt)/dt), θ_link1(t+Δt), and θ_link2(t+Δt), and the process is repeated specified number of times.
FIGS. 36 and 37 are graphs of dynamic analysis results of the manipulator-like device connected by rigid body link shown in FIG. 35. From each dynamic analysis results, as shown in FIGS. 36 and 37, almost the same dynamic characteristics have been obtained as the output results of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
In the equations (32) and (33), when angular displacements θ_link1and θ_link2are small, the equations (34) and (35) can be approximated equations (36) and (37) as below.
(d ²θ_link1(t)/dt ²)≈(12g/7L _link)((3/4)θ_link2−(3/2)θ_link1) (36)
(d ²θ_link2(t)/dt ²)≈(48g/7L _link)((9/16)θ_link1−(1/2)θ_link2) (37)
These equations can give almost same results as those shown in FIGS. 36 and 37. In this case, the similar process to that in the second embodiment can be applied.
(Ninth Application)
FIG. 38 is a schematic diagram of an outline of configuration of particles or celestial bodies, to which the dynamic analysis is applied, influencing force like gravity each other. The equations of motion specific for these particles or celestial bodies influencing force like gravity each other are given as equations (38) to (41) below.
m _1gα_x1 =G(m _1g m _2g /r ²)cos θ (38)
m _1gα_y1 =G(m _1g m _2g /r ²)sin θ (39)
m _2gα_x2 =−G(m _1g m _2g /r ²)cos θ (40)
m _2gα_y2 =−G(m _1g m _2g /r ²)sin θ (41)
where,
m_1g: mass of the Earth (5.97×10²⁴),
m_2g: mass of the Moon (7.35×10²²),
G: gravitational constant (6.67×10⁻¹¹),
α_x1: acceleration of the Earth in x direction,
α_y1: acceleration of the Earth in y direction,
α_x2: acceleration of the Moon in x direction,
α_y2: acceleration of the Moon in y direction,
r: distance between the Earth and the Moon,
and these are shown in the International System of Units.
The r, cos θ, and sin θ in equations (38) to (41) are defined in the equations (42) to (44) below.
r=√((x ₁ −x ₂)²+(y ₁ −y ₂)²) (42)
cos θ=|x₁ −x ₂ |/r (43)
sin θ=|y ₁ −y ₂ |/r (44)
where,
x₁: displacement of the Earth in x direction (initial value=0 0.00×10⁰⁰),
x₂: displacement of the Earth in y direction (initial value=0 0.00×10⁰⁰),
y₁: displacement of the Moon in x direction (initial value=3.84×10⁰⁸),
y₂: displacement of the Moon in y direction (initial value=0.00×10⁰⁰),
v_x1: velocity of the Earth in x direction (initial value=0.00×10⁰⁰),
v_y1: velocity of the Earth in y direction (initial value=0.00×10⁰⁰),
v_y2: velocity of the Moon in x direction (initial value=3.84×10⁰⁰), and
v_y2: velocity of the Moon in y direction (initial value=1.04×10⁰³)
In this case, analysis equations are derived from the equations (38) to (41) and are given as equations (45) to (48) below.
α_x1 =G(m _2g /r ²)cos θ (45)
α_y1 =G(m _2g /r ²)sin θ (46)
α_x2 =−G(m _1g /r ²)cos θ (47)
α_y2 =−G(m _1g /r ²)sin θ (48)
Acceleration α_x1, α_y1, α_x2, and α_y2are calculated assigning parameters, initial displacement of the Earth in x direction x₁(t₀)=0.00×10⁰⁰, initial displacement of the Earth in y direction y₁(t₀)=0.00×10⁰⁰, initial displacement of the Moon in x direction x₂(t₀)=3.84×10⁰⁸, and initial displacement of the Moon in y direction y₂(t₀)=0.00×10⁰⁰to the equations (45) to (48). Then velocity v_x1(t+Δt), v_y1(t+Δt), v_x2(t+Δt), and v_y2(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equation (2). Further displacement x₁(t+Δt), x₂(t+Δt), y₁(t+Δt), and y₂(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated by the equation (3). Acceleration α_x1(t+Δt), a_y1(t+Δt) a_x2(t+Δt), and α_y2(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) are calculated using these obtained v_x1(t+Δt), v_y1(t+Δt), v_x2(t+Δt), v_y2(t+Δt), x₁(t+Δt), x₂(t+Δt), y₁(t+Δt), and y₂(t+Δt), and the process is repeated specified number of times.
FIG. 39 is a graph of a dynamic analysis result of displacement of the Earth in x direction shown in FIG. 38. From this dynamic analysis result, as shown in FIG. 39, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Tenth Application)
The present invention can process a numerical analysis of high precision about three-body problem in a short period of time. Three-body problem means a problem to solve behavior of each of three existing bodies using universal gravitation, and here dynamic analysis is applied to the Sun, the Earth, and the Moon respectively.
FIG. 40 is a diagram of configuration relationship to analyze three-body problem of the Sun, the Earth, and the Moon. Analysis is performed assuming the Sun 101, the Earth 102, and the Moon 103 move in the identical x-y plane (orbital plane of the Earth).
The equation of motion for the Sun 101 is given as below.
F1=M ₁ *a ₁ (49)
where M₁is the mass of the Sun 101 and a₁is an acceleration of the Sun 101. Here a₁is a vector a₁(a_x1, a_y1). Universal gravitation f₁₂by the Earth 102 is given as below.
f ₁₂ ={G·M ₂ ·M ₁/(|P ₂ −P ₁|²)}*{(P ₂ −P ₁)/|P ₂ −P ₁|} (50)
where G is universal gravitation constant, M₂is the mass of the Earth 102, P₁and P₂are position vectors P₁(x₁,y₁) and P₂(x₂,y₂) of the Sun 101 and the Earth 102, respectively. Universal gravitation f₁₂is a vector. Here {G·M₂·M₁/(|P₂−P₁|²)} means absolute value of universal gravitation and {(P₂−P₁)/|P₂−P₁|} means unit vector from the Sun 101 to the Earth 102.
An equation (50) can be transformed into the equation below.
f ₁₂ ={G·M ₂ ·M ₁/(|P ₂ −P ₁|³)}*(P ₂ −P ₁) (51)
Similarly, the universal gravitation constant by the Moon 103 is given as below.
f ₁₃ ={G·M ₃ ·M ₁/(|P ₃ −P ₁|³)}*(P ₃ −P ₁) (52)
where M₃is the mass of the Moon and P₃is a position vector of the Moon P₃(x₃,y₃). The universal gravitation f₁₃is also a vector.
Using the relationship
F ₁ =f ₁₂ +f ₁₃ (53)
, and erasing F₁, f₁₂, and f₁₃from equations (49), and from (51) to (53), the following equation is derived.
a ₁ ={G·M ₂/(|P ₂ −P ₁|³)}*(P ₂ −P ₁)+{G·M ₃/(|P ₃ −P ₁|³)}*(P ₃ −P ₁) (54)
where {G·M₂/(|P₂−P₁|³)}*(P₂−P₁) is an acceleration part by gravitation of the Earth 102 and {G·M₃/(|P₃−P₁|³)}*(P₃−P₁) is an acceleration part by gravitation of the Moon 103.
Similarly the acceleration of the Earth a₂is given as below.
a ₂ ={G·M ₃/(|P ₃ −P ₂|³)}*(P ₃ −P ₂)+{G·M ₁/(|P ₁ −P ₂|³)}*(P ₁ −P ₂) (55)
where {G·M₃/(|P₃−P₂|³)}*(P₃−P₂) is an acceleration part by gravitation of the Moon 103 and {G·M₁/(|P₁−P₂|³)}*(P₁−P₂) is an acceleration part by gravitation of the Sun 101. Here a₂is a vector a₂(a_x2,a_y2). Further the acceleration of the Moon a₃is similarly given as below.
a ₃ ={G·M ₁/(|P ₁ −P ₃|³)}*(P ₁ −P ₃)+{G·M ₂/(|P ₂ −P ₃|³)}*(P ₂ −P ₃) (56)
where {G·M₁/(|P₁−P₃|³)}*(P₁−P₃) is an acceleration part by gravitation of the Sun 101 and {G·M₂/(|P₂−P₃|³)}*(P₂−P₃) is an acceleration part by gravitation of the Earth 102. Here a₃is a vector a₃(a_x3,a_y3).
The simultaneous differential equations (54), (55), and (56) can be solved using a dynamic analysis as three-body problem by incrementing the analytical small time increment Δt by one hour and finally simultaneous differential equations with six unknowns are to be solved all together because every two unknowns are included in each simultaneous equation. Now initial values are explained. Firstly time t=0 is taken at the moment when the Sun 101, the Earth 102, and the Moon 103 are on the straight-line, the Sun 101 is deemed as the origin, the directions of the Earth 102 and the Moon 103 are taken in +x direction, and the velocity vectors of the Earth 102 and the Moon 103 are taken in +y direction.
Conditions and parameters are given below.


Initial distance between the Sun 101 and	1.496 × 10¹¹	(m)
the Earth 102:
Initial distance between the Earth 102	3.844 × 10⁸	(m)
and the Moon 103:
Initial velocity of the Earth 102:	2.978 × 10⁴	(m/s)
Initial period of the Moon 103:	27.3217	(day)
Initial relative velocity of the Moon 103	1023	(m/s)
to the Earth 102:
Mass of the Sun 101 M¹:	1.989 × 10³⁰	(kg)
Mass of the Earth 102 M²:	5.974 × 10²⁴	(kg)
Mass of the Moon 103 M³:	7.347 × 10²²	(kg)
Universal gravitation constant G:	6.67259 × 10⁻¹¹	(N · m²/kg²).

To solve this three-body problem, accelerations a_x1, a_y1, a_x2, a_y2, a_x3, a_y3are calculated assigning each initial value to the equations (54), (55), and (56). Then velocities v_x1(t+Δt), v_y1(t+Δt), v_x2(t+Δt), v_y2(t+Δt) v_x3(t+Δt), and v_y3(t+Δt) after lapse of the analytical small time increment Δt(=3600(s)) are calculated by the equation (2). Further displacements x₁(t+Δt), y₁(t+Δt), x₂(t+Δt), y₂(t+Δt), x₃(t+Δt), and y₃(t+Δt) after lapse of the analytical small time increment Δt(=3600(s)) are calculated by the equation (3). Accelerations (x₁(t+Δt), a_y1(t+Δt), a_x2(t+Δt), a_y2(t+Δt), α_x3(t+Δt), and α_y3(t+Δt) after lapse of the analytical small time increment Δt(=3600 (s)) are calculated using these obtained v_x1(t+Δt), v_y1(t+Δt), v_x2(t+Δt), v_y2(t+Δt), v_x3(t+Δt), v_y3(t+Δt), x₁(t+Δt), y₁(t+Δt), x₂(t+Δt), y₂(t+Δt), x₃(t+Δt), and y₃(t+Δt), and the process is repeated specified number of times.
FIG. 41 is a graph of a dynamic analysis result of the orbit of the Moon in three-body problem above, placing the Earth in the center. From this dynamic analysis result, as shown in FIG. 41, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
FIG. 42 is a graph of a dynamic analysis result of the orbit of the Moon using conventional Euler's method, placing the Earth in the center. The dynamic analysis whose result is shown in FIG. 42 was performed with the same analytical small time increment Δt=3600 second as in tenth application and it shows that the orbit of the Moon diverges and no numerical analysis result of high precision could be obtained.
(Eleventh Application)
FIG. 43 is a schematic diagram of an outline of configuration of a surge tank, to which the dynamic analysis is applied, with a centrifugal pump. The equation of motion specific for water level change of the surge tank is given as an equation (57) below.
((l _surge A _surge)/(g·A _tank))α_surge=(A _tank h′(Q _surge)−2V _surge Q _surge A _surge /A _tank ²)v _surge −x _surge (57)
where,
g: acceleration of gravity (9.807),
l_surge: length of the conduit line (10.000),
A_surge: cross-section area of the conduit line (0.008),
A_tank: cross-section area of the surge tank (0.100),
Q_surge: flow rate (0.012),
h′(Q_surge)=2.67×10⁹Q_surge ³−1.2×10⁸Q_surge ²+1.2×10⁵Q_surge+472.487: (head-discharge curve differentiated by flow rate),
V_surge: pressure loss coefficient (0.144),
x_surge: water level change of the surge tank (initial value=0.02),
v_surge: velocity of water level change of the surge tank (initial value=0.00),
α_surge: acceleration of water level change of the surge tank, and these are shown in the International System of Units.
In this case, an analysis equation is derived from the equation (57) and is given as an equation (58) below.
α_surge=((g·A _tank)/(l _surge A _surge)){(A _tank h′(Q _surge)−2V _surge Q _surge A _surge /A _tank ²)v _surge −x _surge} (58)
Acceleration of water level change of the surge tank α_surgeis calculated assigning parameters, initial water level change of the surge tank x_surge(t₀)=0.02, and initial velocity of water level of the surge tank v_surge(t₀)=0.00 to the equation (58). Then velocity of water level change v_surge(t+Δt) after lapse of analytical small time increment Δt(=0.01(s)) is calculated by the equation (2). Further water level change x_surge(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) is calculated by the equation (3). Acceleration of water level change α_surge(t+Δt) after lapse of the analytical small time increment Δt(=0.01(s)) is calculated using these obtained v_surge(t+Δt) and x_surge(t+Δt), and the process is repeated specified number of times.
FIG. 44 is a graph of a dynamic analysis result of the change of water level of the surge tank shown in FIG. 43. From this dynamic analysis result, as shown in FIG. 44, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
(Twelfth Application)
FIG. 45 is a schematic diagram of an outline of configuration of an LCR circuit to which the dynamic analysis is applied. The equation of motion specific for electric charge of the LCR circuit is given as equation (59) below.
L _lcr(d ² q/dt ²)+R _lcr i+q/c _lcr =E _lcr (59)
where,
L_lcr: inductance (0.20),
R_lcr: resistance (5.00),
C_lcr: capacitance (0.002),
E_lcr: voltage (initial value=5.00),
q: electric charge (initial value=0.00),
i: current (initial value=0.00),
and these are shown in the International System of Units.
In this case, the analysis equation is derived from the equation (59) and is given as an equation (60) below.
(d ² q/dt ²)=(1/L _lcr)(−R _lcr i−q/c _lcr +E _lcr) (60)
Second-order time derivative of the electric charge (d²q/dt²) is calculated assigning parameters, initial electric charge q(t₀)=0.00, and initial current i(t₀)=0.00 to the equation (60). Then current i(t+Δt) after lapse of analytical small time increment Δt(=0.001(s)) is calculated by the equation (2). Further electric charge q(t+Δt) after lapse of the analytical small time increment Δt(=0.001(s)) is calculated by the equation (3). Second-order time derivative of the electric charge (d²q(t+Δt)/dt²) after lapse of the analytical small time increment Δt(=0.001(s)) is calculated using these obtained i(t+Δt) and q(t+Δt), and the process is repeated specified number of times.
FIG. 46 is a graph of a dynamic analysis result of change of electric charge in the LCR circuit shown in FIG. 45. From this dynamic analysis result, as shown in FIG. 46, almost the same dynamic characteristic has been obtained as the output result of a dynamic analysis device that performs time consuming analysis of high precision combining many processes like correction process including the feedback processing.
The dynamic analysis device that is an example of the numerical analysis device can be implemented by running prepared programs on a computer system like a personal computer or workstation. Here a computer system is explained that runs dynamic analysis programs that have similar function as the said dynamic analysis devices.
FIG. 47 is a schematic diagram of a computer system configuration that implements the dynamic analysis device and FIG. 48 is a block diagram of configuration of the main body of the computer system shown in FIG. 47.
Referring to FIG. 47, a computer system 200 includes a main unit 201, a display 202 to display information such as images on a display screen 202 a according to indication by the main unit 201, a keyboard 203 to input various information to the computer system 200, and a mouse 204 to point arbitrary location on the display screen 202 a of the display 202.
The main unit 201 of the computer system 200 includes, as shown in FIG. 48, a CPU 221, a RAM 222, a ROM 223, a hard disk drive (HDD) 224, a CD-ROM drive 225 that accepts a CD-ROM 209, an FD drive 226 that accepts a flexible disk (FD) 208, an I/O interface 227 to connect to the display 202, to the keyboard 203, and to the mouse 204, and a LAN interface 228 to connect to a local area network/wide area network (LAN/WAN) 206.
A modem 205 for connection to a public telephone line 207 for internet and the like is connected to the computer system 200 and furthermore a computer system (PC) 211, a server 212, and a printer 213 and the like are connected via the LAN interface 228 and the LAN/WAN 206 to the computer system 200.
The computer system 200 implements a dynamic analysis device in the first and the second embodiment by reading out and executing dynamic analysis programs recorded on specified recording media.
The specified recording media include, besides “removable physical media” such as the flexible disk (FD) 208, the CD-ROM 209, a magneto-optical disk, a DVD, an optical-magnetic disk, and an IC card, all recording media to store dynamic analysis programs that can be read out by the computer system 200, such as the hard disk drive (HDD) 224 that is installed inside or outside of the computer system 200, “irremovable physical media” like the RAM 222 and the ROM 223, “communication media” that can retain programs for a short period of time for transmitting the programs, such as the public telephone line 207 connected via the modem 205, or the LAN/WAN 206 to which another the computer system 211 and the server 212 are connected.
Namely, the dynamic analysis programs are recorded on the “removable physical media”, “irremovable physical media” and “communication media” to be legible to a computer, and the computer system 200 reads out the dynamic analysis programs from the recording media and runs them to implement dynamic analysis devices in the first and the second embodiment.
The dynamic analysis programs are not limited to be run on the computer system 200 but also can be run on another computer system 211 or the server 212, and even when they work together the present invention can be applied as well.
Each component of the devices shown in the figures only represents the concept of the function and shall not necessarily be implemented as shown in the figures. In other words, a concrete aspect of partition or integration of each device shall not be limited as illustrated but it is allowed to configure all or a part of the device being functionally or physically partitioned or integrated in arbitrary unit depending on various conditions of loads and use.
Each process function that is performed on each device can all or partly be implemented by CPU and analysis programs that run on the CPU, and can also be implemented as wired logic hardware.

INDUSTRIAL APPLICABILITY

Thus the present invention is useful for numerical analysis devices as well as numerical analysis programs that are applied in control field to analyze various aspects of motion in natural phenomena, mechanical motion and the like.

Claims

1-36. (canceled)

37. A numerical analysis device that performs a numerical analysis of second-order differential equation for a target object to be analyzed, the numerical analysis device comprising:

a setting unit that performs a setting of parameters including an initial condition of the second-order differential equation and a setting of an analytical small variable value for the numerical analysis;

a second-derivative calculation unit that calculates a second derivative of the second-order differential equation at a reference variable value that is a variable at the start of the numerical analysis;

a first-derivative calculation unit that performs a multiplication of the second derivative calculated by the second-derivative calculation unit by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputs a value obtained by the addition as a first derivative after an increment of the analytical small variable value;

a physical-value calculation unit that performs a multiplication of the first derivative calculated by the first-derivative calculation unit by the analytical small variable value and an addition of a physical value of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputs a value obtained by the addition as a physical value after an increment of the analytical small variable value; and

an operation control unit that repeats a control of setting a variable value obtained by adding the analytical small variable value to the reference variable value as a new reference variable value and causing the second-derivative calculation unit, the first-derivative calculation unit, and the physical-value calculation unit to calculate the second derivative, the first derivative, and the physical value, using the new reference variable value, respectively.

38. The numerical analysis device according to claim 37, wherein

the analytical small variable value and the reference variable value are an analytical small time and reference time, respectively, and

the first derivative and the second derivative are a first-order time derivative and a second-order time derivative, respectively.

39. The numerical analysis device according to claim 37, wherein

the setting unit sets a termination condition of the numerical analysis, and

when the termination condition of the numerical analysis is satisfied, the operation control unit performs a control of terminating the numerical analysis.

40. The numerical analysis device according to claim 37, further comprising an output processing unit that outputs the physical value of each reference variable value calculated by the physical-value calculation unit, in association with the increment of the analytical small variable value.

41. The numerical analysis device according to claim 40, wherein the output processing unit generates a graph of the physical value for each reference variable value calculated by the physical-value calculation unit corresponding to the increment of the analytical small variable value.

42. The numerical analysis device according to claim 37, wherein the second-order differential equation includes a variable value term.

43. The numerical analysis device according to claim 38, wherein the second-order differential equation includes a time term.

44. The numerical analysis device according to claim 37, wherein the second-order differential equation is a simultaneous equation with multiple unknowns.

45. The numerical analysis device according to claim 37, wherein the numerical analysis is performed using a spreadsheet program.

46. A numerical analysis device that performs a numerical analysis of second-order differential equation for a target object to be analyzed, the numerical analysis device comprising:

a physical-value calculation unit that performs a multiplication of the second derivative calculated by the second-derivative calculation unit by a square value of the analytical small variable value, a multiplication of the first derivative calculated by the first-derivative calculation unit by the analytical small variable value, and an addition of values obtained by the multiplications and a physical value of the second-order differential equation at the reference variable value, and outputs a value obtained by the addition as a physical value after an increment of the analytical small variable value; and

47. The numerical analysis device according to claim 46, wherein

48. The numerical analysis device according to claim 46, wherein

the setting unit sets a termination condition of the numerical analysis, and

49. The numerical analysis device according to claim 46, further comprising an output processing unit that outputs the physical value of each reference variable value calculated by the physical-value calculation unit, in association with the increment of the analytical small variable value.

50. The numerical analysis device according to claim 49, wherein the output processing unit generates a graph of the physical value for each reference variable value calculated by the physical-value calculation unit corresponding to the increment of the analytical small variable value.

51. The numerical analysis device according to claim 46, wherein the second-order differential equation includes a variable value term.

52. The numerical analysis device according to claim 47, wherein the second-order differential equation includes a time term.

53. The numerical analysis device according to claim 46, wherein the second-order differential equation is a simultaneous equation with multiple unknowns.

54. The numerical analysis device according to claim 46, wherein the numerical analysis is performed using a spreadsheet program.

55. A computer-readable recording medium that stores therein a numerical analysis program for a numerical analysis of second-order differential equation for a target object to be analyzed, the numerical analysis program causing a computer to execute:

parameter setting including performing a setting of parameters including an initial condition of the second-order differential equation and a setting of an analytical small variable value for the numerical analysis;

second-derivative calculating including calculating a second derivative of the second-order differential equation at a reference variable value that is a variable at the start of the numerical analysis;

first-derivative calculating including

performing a multiplication of the second derivative calculated at the second-derivative calculating by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and

outputting a value obtained by the addition as a first derivative after an increment of the analytical small variable value;

physical-value calculating including

performing a multiplication of the first derivative calculated at the first-derivative calculating by the analytical small variable value and an addition of a physical value of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and

outputting a value obtained by the addition as a physical value after an increment of the analytical small variable value; and

operation controlling including repeating a control of

setting a variable value obtained by adding the analytical small variable value to the reference variable value as a new reference variable value, and

causing the second-derivative calculating, the first-derivative calculating, and the physical-value calculating to calculate the second derivative, the first derivative, and the physical value, using the new reference variable value, respectively.

56. The computer-readable recording medium according to claim 55, wherein

57. The computer-readable recording medium according to claim 55, wherein

the parameter setting includes setting a termination condition of the numerical analysis, and

when the termination condition of the numerical analysis is satisfied, the operation controlling includes performing a control of terminating the numerical analysis.

58. The computer-readable recording medium according to claim 55, wherein the numerical analysis program further causes the computer to execute output processing including outputting the physical value of each reference variable value calculated at the physical-value calculating, in association with the increment of the analytical small variable value.

59. The computer-readable recording medium according to claim 58, wherein the output processing includes generating a graph of the physical value for each reference variable value calculated at the physical-value calculating corresponding to the increment of the analytical small variable value.

60. The computer-readable recording medium according to claim 55, wherein the second-order differential equation includes a variable value term.

61. The computer-readable recording medium according to claim 56, wherein the second-order differential equation includes a time term.

62. The computer-readable recording medium according to claim 55, wherein the second-order differential equation is a simultaneous equation with multiple unknowns.

63. The computer-readable recording medium according to claim 55, wherein the numerical analysis is performed using a spreadsheet program.

64. A computer-readable recording medium that stores therein a numerical analysis program for a numerical analysis of second-order differential equation for a target object to be analyzed, the numerical analysis program causing a computer to execute:

first-derivative calculating including performing a multiplication of the second derivative calculated at the second-derivative calculating by the analytical small variable value and an addition of the first derivative of the second-order differential equation at the reference variable value to a value obtained by the multiplication, and outputting a value obtained by the addition as a first derivative after an increment of the analytical small variable value;

physical-value calculating including performing a multiplication of the second derivative calculated at the second-derivative calculating by a square value of the analytical small variable value, a multiplication of the first derivative calculated at the first-derivative calculating by the analytical small variable value, and an addition of values obtained by the multiplications and a physical value of the second-order differential equation at the reference variable value, and outputting a value obtained by the addition as a physical value after an increment of the analytical small variable value; and

operation controlling including repeating a control of

65. The computer-readable recording medium according to claim 64, wherein

66. The computer-readable recording medium according to claim 64, wherein

67. The computer-readable recording medium according to claim 64, wherein the numerical analysis program further causes the computer to execute output processing including outputting the physical value of each reference variable value calculated at the physical-value calculating, in association with the increment of the analytical small variable value.

68. The computer-readable recording medium according to claim 67, wherein the output processing includes generating a graph of the physical value for each reference variable value calculated at the physical-value calculating corresponding to the increment of the analytical small variable value.

69. The computer-readable recording medium according to claim 64, wherein the second-order differential equation includes a variable value term.

70. The computer-readable recording medium according to claim 65, wherein the second-order differential equation includes a time term.

71. The computer-readable recording medium according to claim 64, wherein the second-order differential equation is a simultaneous equation with multiple unknowns.

72. The computer-readable recording medium according to claim 64, wherein the numerical analysis is performed using a spreadsheet program.