(57)【要約】解析対象の２階微分方程式を動解析する動解析装置１であって、前記２階微分方程式の初期条件を含むパラメータ設定および前記動解析の解析微小時間Δｔを設定する算入モジュール２１と、解析開始時刻である基準時刻における前記２階微分方程式の２階時間微分量α（ｔ）を算出する二階時間微分量算出モジュール２２と、前記解析微小時間経過後の１階時間微分量ｖ（ｔ＋Δｔ）＝ｖ（ｔ）＋α（ｔ）・Δｔを算出する一階時間微分量算出モジュール２３と、前記解析微小時間経過後の物理量ｘ（ｔ＋Δｔ）＝ｘ（ｔ）＋ｖ（ｔ＋Δｔ）・Δｔを算出する物理量算出モジュール２４と、時刻（ｔ＋Δｔ）を新たな基準時刻として設定し、各モジュールに対してそれぞれ設定された基準時刻を用いて前記２階時間微分量、前記１階時間微分量、および前記物理量を算出させる制御を繰り返し行う演算制御モジュール２０とを備える。【特許請求の範囲】
【請求項１】
解析対象の２階微分方程式を数値解析する数値解析装置であって、
前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、
解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、
前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、
前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、
前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、
を備えたことを特徴とする数値解析装置。
【請求項２】
前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項１に記載の数値解析装置。
【請求項３】
前記設定手段は、前記数値解析の終了条件を設定し、
前記演算制御手段は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項１に記載の数値解析装置。
【請求項４】
前記物理量算出手段によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手段を備えたことを特徴とする請求項１に記載の数値解析装置。
【請求項５】
前記出力処理手段は、前記解析微小変数量の増分に対応する、物理量算出手段によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項４に記載の数値解析装置。
【請求項６】
前記２階微分方程式は、変数量項を含むことを特徴とする請求項１に記載の数値解析装置。
【請求項７】
前記２階微分方程式は、時間項を含むことを特徴とする請求項２に記載の数値解析装置。
【請求項８】
前記２階微分方程式は、多元連立方程式であることを特徴とする請求項１に記載の数値解析装置。
【請求項９】
表計算プログラムを用いて数値解析を行うことを特徴とする請求項１に記載の数値解析装置。
【請求項１０】
解析対象の２階微分方程式を数値解析する数値解析装置であって、
前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、
解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、
前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、
前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、
前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、
を備えたことを特徴とする数値解析装置。
【請求項１１】
前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項１０に記載の数値解析装置。
【請求項１２】
前記設定手段は、前記数値解析の終了条件を設定し、
前記演算制御手段は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項１０に記載の数値解析装置。
【請求項１３】
前記物理量算出手段によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手段を備えたことを特徴とする請求項１０に記載の数値解析装置。
【請求項１４】
前記出力処理手段は、前記解析微小変数量の増分に対応する、物理量算出手段によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項１３に記載の数値解析装置。
【請求項１５】
前記２階微分方程式は、変数量項を含むことを特徴とする請求項１０に記載の数値解析装置。
【請求項１６】
前記２階微分方程式は、時間項を含むことを特徴とする請求項１１に記載の数値解析装置。
【請求項１７】
前記２階微分方程式は、多元連立方程式であることを特徴とする請求項１０に記載の数値解析装置。
【請求項１８】
表計算プログラムを用いて数値解析を行うことを特徴とする請求項１０に記載の数値解析装置。
【請求項１９】
解析対象の２階微分方程式を数値解析する数値解析プログラムであって、
前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、
解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、
前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、
前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、
前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、
をコンピュータに実行させることを特徴とする数値解析プログラム。
【請求項２０】
前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項１９に記載の数値解析プログラム。
【請求項２１】
前記設定手順は、前記数値解析の終了条件を設定し、
前記演算制御手順は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項１９に記載の数値解析プログラム。
【請求項２２】
前記物理量算出手順によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手順を備えたことを特徴とする請求項１９に記載の数値解析プログラム。
【請求項２３】
前記出力処理手順は、前記解析微小変数量の増分に対応する、物理量算出手順によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項２２に記載の数値解析プログラム。
【請求項２４】
前記２階微分方程式は、変数量項を含むことを特徴とする請求項１９に記載の数値解析プログラム。
【請求項２５】
前記２階微分方程式は、時間項を含むことを特徴とする請求項２０に記載の数値解析プログラム。
【請求項２６】
前記２階微分方程式は、多元連立方程式であることを特徴とする請求項１９に記載の数値解析プログラム。
【請求項２７】
表計算プログラムを用いて数値解析を行うことを特徴とする請求項１９に記載の数値解析プログラム。
【請求項２８】
解析対象の２階微分方程式を数値解析する数値解析プログラムであって、
前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、
解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、
前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、
前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、
前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、
をコンピュータに実行させることを特徴とする数値解析プログラム。
【請求項２９】
前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項２８に記載の数値解析プログラム。
【請求項３０】
前記設定手順は、前記数値解析の終了条件を設定し、
前記演算制御手順は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項２８に記載の数値解析プログラム。
【請求項３１】
前記物理量算出手順によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手順を備えたことを特徴とする請求項２８に記載の数値解析プログラム。
【請求項３２】
前記出力処理手順は、前記解析微小変数量の増分に対応する、物理量算出手順によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項３１に記載の数値解析プログラム。
【請求項３３】
前記２階微分方程式は、変数量項を含むことを特徴とする請求項２８に記載の数値解析プログラム。
【請求項３４】
前記２階微分方程式は、時間項を含むことを特徴とする請求項２９に記載の数値解析プログラム。
【請求項３５】
前記２階微分方程式は、多元連立方程式であることを特徴とする請求項２８に記載の数値解析プログラム。
【請求項３６】
表計算プログラムを用いて数値解析を行うことを特徴とする請求項２８に記載の数値解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】  この発明は、自然現象や機械運動などにおける種々の運動態様を解析することができる数値解析装置および数値解析プログラムに関するものである。
【背景技術】
【０００２】  従来から、常微分方程式を用いた数値解析法には、区間の初期変位、初期速度、初期加速度などを用いて、つぎの時刻の変位や速度を求めるオイラー法が広く用いられている（特許文献１，２参照）。【０００３】
【特許文献１】特開２００４−０６２６７６号公報
【特許文献２】特開２００２−２５８９３３号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】  しかしながら、上述したオイラー法は、微小時間Δｔ毎に処理ステップを重ねる毎にその解析結果に比較的大きな誤差が生じるという問題点があった。さらに、この誤差を抑えるために微小時間Δｔを小さくすると、解析結果を得るために膨大な処理時間を要し、実用に供しないという問題点があった。【０００５】  この発明は、上記に鑑みてなされたものであって、解析結果に大きな誤差を生じさせず、簡易かつ短時間に数値解析処理を行うことができる数値解析装置および数値解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】  上述した課題を解決し、目的を達成するために、この発明にかかる数値解析装置は、解析対象の２階微分方程式を数値解析する数値解析装置であって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、を備えたことを特徴とする。【０００７】  また、この発明にかかる数値解析装置は、解析対象の２階微分方程式を数値解析する数値解析装置であって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、を備えたことを特徴とする。【０００８】  また、この発明にかかる数値解析装置は、上記の発明において、前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする。【０００９】  また、この発明にかかる数値解析装置は、上記の発明において、前記設定手段は、前記数値解析の終了条件を設定し、前記演算制御手段は、前記数値解析の終了条���を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする。【００１０】  また、この発明にかかる数値解析装置は、上記の発明において、前記物理量算出手段によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手段を備えたことを特徴とする。【００１１】  また、この発明にかかる数値解析装置は、上記の発明において、前記出力処理手段は、前記解析微小変数量の増分に対応する、物理量算出手段によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする。【００１２】  また、この発明にかかる数値解析装置は、上記の発明において、前記２階微分方程式は、変数量項を含むことを特徴とする。【００１３】  また、この発明にかかる数値解析装置は、上記の発明において、前記２階微分方程式は、時間項を含むことを特徴とする。【００１４】  また、この発明にかかる数値解析装置は、上記の発明において、前記２階微分方程式は、多元連立方程式であることを特徴とする。【００１５】  また、この発明にかかる数値解析装置は、上記の発明において、表計算プログラムを用いて数値解析を行うことを特徴とする。【００１６】  また、この発明にかかる数値解析プログラムは、解析対象の２階微分方程式を数値解析する数値解析プログラムであって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、をコンピュータに実行させることを特徴とする。【００１７】  また、この発明にかかる数値解析プログラムは、解析対象の２階微分方程式を数値解析する数値解析プログラムであって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、をコンピュータに実行させることを特徴とする。【００１８】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする。【００１９】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記設定手順は、前記数値解析の終了条件を設定し、前記演算制御手順は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする。【００２０】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記物理量算出手順によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手順を備えたことを特徴とする。【００２１】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記出力処理手順は、前記解析微小変数量の増分に対応する、物理量算出手順によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする。【００２２】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記２階微分方程式は、変数量項を含むことを特徴とする。【００２３】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記２階微分方程式は、時間項を含むことを特徴とする。【００２４】  また、この発明にかかる数値解析プログラムは、上記の発明において、前記２階微分方程式は、多元連立方程式であることを特徴とする。【００２５】  また、この発明にかかる数値解析プログラムは、上記の発明において、表計算プログラムを用いて数値解析を行うことを特徴とする。
【発明の効果】
【００２６】  この発明にかかる数値解析装置および数値解析プログラムによれば、設定手段が、２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定し、２階微分量算出手段が、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出し、１階微分量算出手段が、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出し、物理量算出手段が、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出し、演算制御手段が、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行うようにしているので、解析結果に大きな誤差を生じさせず、簡易かつ短時間に数値解析処理を行うことができるという効果を奏する。【００２７】  また、この発明にかかる数値解析装置および数値解析プログラムによれば、設定手段が、２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定し、２階微分量算出手段が、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出し、１階微分量算出手段が、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出し、物理量算出手段が、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出し、演算制御手段が、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行うようにしているので、解析結果に大きな誤差を生じさせず、簡易かつ短時間に数値解析処理を行うことができるという効果を奏する。
【図面の簡単な説明】
【００２８】
【図１】図１は、この発明の実施の形態１にかかる数値解析装置である動解析装置の概要構成を示すブロック図である。
【図２】図２は、図１に示した動解析装置による動解析処理手順を示すフローチャートである。
【図３】図３は、この発明の実施の形態１による動解析処理の概念を示す図である。
【図４】図４は、従来のオイラー法による動解析処理の概念を示す図である。
【図５】図５は、従来の修正オイラー法による動解析処理の概念を示す図である。
【図６】図６は、解析方程式の対象である、ばね−質点系の概念構成を示す模式図である。
【図７】図７は、図６に示したばね−質点系の動解析結果を示す図である。
【図８】図８は、図７に示したばね−質点系を動解析する表計算プログラムの表示画面の一例を示す図である。
【図９】図９は、図６に示したばね−質点系に対する動解析処理を表計算プログラムによって行った一例を示す図である。
【図１０】図１０は、図９に示した結果をグラフ化した図である。
【図１１】図１１は、同じばね−質点系に対する動解析処理をオイラー法によって行った一例を示す図である。
【図１２】図１２は、図１１に示した動解析結果をグラフ化した図である。
【図１３】図１３は、同じばね−質点系に対する動解析を修正オイラー法によって行った結果をグラフ化した図である。
【図１４】図１４は、この発明の実施の形態２である動解析装置の概要構成を示すブロック図である。
【図１５】図１５は、図１４に示した動解析装置による動解析処理手順を示すフローチャートである。
【図１６】図１６は、この発明の実施の形態２による動解析処理の概念を示す図である。
【図１７】図１７は、動解析対象であるばね質点ダッシュポット系の概要構成を示す模式図である。
【図１８】図１８は、図１７に示した質点ダッシュポット系の動解析結果を示す図である。
【図１９】図１９は、動解析対象である、単振動と回転運動とが伴った運搬機の概要構成を示す模式図である。
【図２０】図２０は、図１９に示した運搬機の変位の動解析結果を示す図である。
【図２１】図２１は、図１９に示した運搬機の角変位の動解析結果を示す図である。
【図２２】図２２は、動解析対象である、ねじり軸、継手、ギア等の動力伝達器で連結された複数の回転体の概要構成を示す模式図である。
【図２３】図２３は、図２２に示した動力伝達器で結合された各回転体の角変位の動解析結果を示す図である。
【図２４】図２４は、図２２に示した動力伝達器で結合された各回転体の角変位の動解析結果を示す図である。
【図２５】図２５は、動解析対象である、調速機のすべり子およびタービンの概要構成を示す模式図である。
【図２６】図２６は、図２５に示した調速機のすべり子の変位の動解析結果を示す図である。
【図２７】図２７は、図２５に示したタービンの角変位の動解析結果を示す図である。
【図２８】図２８は、動解析対象である、旋削機などのバイトの概要構成を示す模式図である。
【図２９】図２９は、図２８に示した旋削機などのバイト変位の動解析結果を示す図である。
【図３０】図３０は、動解析対象である、クランクを有する往復機械の概要構成を示す模式図である。
【図３１】図３１は、図３０に示したクランクを有する往復機械の動解析結果を示す図である。
【図３２】図３２は、図３０に示した往復機械のピストンの動解析結果を示す図である。
【図３３】図３３は、動解析対象である、ブレーキディスクの概要構成を示す模式図である。
【図３４】図３４は、図３３に示したブレーキディスクの動解析結果を示す図である。
【図３５】図３５は、動解析対象である、剛体リンクで接続されたマニピュレータ状機器の概要構成を示す模式図である。
【図３６】図３６は、図３５に示した剛体リンクで接続されたマニピュレータ状機器の動解析結果を示す図である。
【図３７】図３７は、図３５に示した剛体リンクで接続されたマニピュレータ状機器の動解析結果を示す図である。
【図３８】図３８は、動解析対象である、重力等を相互に及ぼし合う粒子あるいは天体の概要構成を示す模式図である。
【図３９】図３９は、図３８に示した地球のｘ方向変位の動解析結果を示す図である。
【図４０】図４０は、太陽、地球、月の三体問題を解析する座標系および配置関係を示す図である。
【図４１】図４１は、３体問題における地球を中心とした月の軌道の動解析結果を示す図である。
【図４２】図４２は、従来のオイラー法を用いて地球を中心とした月の軌道の動解析結果を示す図である。
【図４３】図４３は、動解析対象である、遠心ポンプを有するシステムのサージタンクの概要構成を示す模式図である。
【図４４】図４４は、図４３に示したサージタンクの水位変動の動解析結果を示す図である。
【図４５】図４５は、動解析対象である、ＬＣＲ回路の概要構成を示す模式図である。
【図４６】図４６は、図４５に示したＬＣＲ回路の電荷変動の動解析結果を示す図である。
【図４７】図４７は、上述した動解析装置を実現するコンピュータシステムの構成を示すシステム構成図である。
【図４８】図４８は、図４７に示したコンピュータシステムにおける本体部の構成を示すブロック図である。
【符号の説明】
【００２９】    １  動解析装置
１０  入力部
１１  出力部
１２  演算部
１３  出力処理部
１４  記憶部
２０  演算制御モジュール
２１  算入モジュール
２２  二階時間微分量算出モジュール
２３  一階時間微分量算出モジュール
２４，２４ａ  物理量算出モジュール
【発明を実施するための最良の形態】
【００３０】  以下、この発明を実施するための最良の形態である数値解析装置の一例としての動解析装置および数値解析プログラムの一例としての動解析プログラムについて説明する。【００３１】（実施の形態１）
図１は、この発明の実施の形態１にかかる数値解析装置の一例である動解析装置の概要構成を示すブロック図である。図１に示すように、この動解析装置１は、キーボードやポインティングデバイスによって実現され、各種情報を入力する入力部１０と、液晶ディスプレイやプリンタなどによって実現され、処理結果を出力する出力部１１と、各種プログラムやデータを格納するとともに、演算結果などを一時的に格納する記憶部１４と、動解析処理の実行制御を行う制御部Ｃとを有する。【００３２】  制御部Ｃは、演算部１２と出力処理部１３とを有し、演算部１２が動解析処理を行い、出力処理部１３は、動解析処理された結果を出力部１１に出力させるための処理を行う。演算部１２は、演算制御モジュール２０、算入モジュール２１、二階時間微分量算出モジュール２２、一階時間微分量算出モジュール２３、および物理量算出モジュール２４を有する。【００３３】  ここで、図２に示したフローチャートを参照して制御部Ｃの制御処理手順について説明する。図１および図２において、まず入力部１０および出力部１１を用い、算入モジュール２２による解析情報の設定入力を行う（ステップＳ１０１）。解析情報とは、解析対象となる解析方程式α（ｔ）＝ｆ（ｘ（ｔ），ｖ（ｔ））、解析微小時間Δｔ、さらには、パラメータ、初期条件、終了条件などが含まれる。ここでいう解析方程式は、二階常微分方程式であり、α（ｔ）は、加速度などの二階時間微分量を示し、ｘ（ｔ）は、位置などの物理量を示し、ｖ（ｔ）は、速度などの一階時間微分量を示している。【００３４】  二階時間微分量算出モジュール２２は、初期条件をもとに基準時刻ｔにおける加速度α（ｔ）を算出する（ステップＳ１０２）。その後、一階時間微分量算出モジュール２３は、算出された加速度α（ｔ）に解析微小時間Δｔを乗算し、この乗算値に基準時刻ｔの速度ｖ（ｔ）を加算し、この加算値を解析微小時間Δｔ経過後の速度ｖ（ｔ＋Δｔ）として算出する（ステップＳ１０３）。すなわち、ｖ（ｔ＋Δｔ）＝ｖ（ｔ）＋α（ｔ）・Δｔを算出する。【００３５】  その後、物理量算出モジュール２４は、算出された速度ｖ（ｔ＋Δｔ）に解析微小時間Δｔを乗算し、この乗算値に基準時刻ｔの変位ｘ（ｔ）を加算し、この加算値を解析微小時間Δｔ経過後の変位ｘ（ｔ＋Δｔ）として算出する。すなわち、ｘ（ｔ＋Δｔ）＝ｘ（ｔ）＋ｖ（ｔ＋Δｔ）・Δｔを算出する。【００３６】  その後、演算部１２は、ステップＳ１０２〜Ｓ１０４で算出した値などを記憶部１４に記憶させるとともに、出力処理部１３を介して出力部１１から出力できる出力処理を行う（ステップＳ１０５）。さらに、演算部１２は、ｔが所定値を超えたか否かを判断する（ステップＳ１０６）。ここで、所定値とは基準時刻ｔからの経過時間である。この経過時間は、この後繰り返す処理回数をｎとすると、Δｔ×ｎとして示すことができる。ｔが所定値を超えていない場合（ステップＳ１０６，Ｎｏ）には、基準時刻ｔに解析微小時間Δｔを加えた値（ｔ＋Δｔ）を新たな基準時刻として設定してステップＳ１０２に移行し（ステップＳ１０７）、上述した処理を繰り返す。一方、ｔが所定値を超えた場合（ステップＳ１０６，Ｙｅｓ）には、本処理を終了する。【００３７】  すなわち、この動解析装置１による動解析処理は、解析微小時間Δｔ後における変位ｘ（ｔ＋Δｔ）を次のようにして求める。まず、基準時刻ｔにおける解析方程式である加速度α（ｔ）を、基準時刻ｔにおける物理量である変位ｘ（ｔ）および基準時刻ｔにおける一階時間微分量である速度ｖ（ｔ）をもとに、式（１）によって算出する。その後、この解析微小時間Δｔ経過後における速度ｖ（ｔ＋Δｔ）を、式（２）をもとに求め、さらに解析微小時間Δｔ経過後における変位ｘ（ｔ＋Δｔ）を、式（３）をもとに求める。
2006103996.tif 000002
【００３８】  ところで、従来の動解析装置では、オイラー法を用いるものがあり、このオイラー法では、上述した式（３）におけるｖ（ｔ＋Δｔ）をｖ（ｔ）とする次式（４）によって求めていた。すなわち、
2006103996.tif 000003
によって変位（ｔ＋Δｔ）を求めていた。しかし、このオイラー法では、単振動などの動解析を行うと発散してしまう場合があった。そこで、式（３）を次式（５）によって求める修正オイラー法が考え出された。すなわち、
2006103996.tif 000004
である。この修正オイラー法では、ｖ（ｔ＋Δｔ）を、基準時刻ｔと、基準時刻ｔから解析微小時間Δｔ経過後の時刻（ｔ＋Δｔ）との間の平均速度（１／２・（ｖ（ｔ）＋ｖ（ｔ＋Δｔ）））に置き換えて動解析を行うようにしている。この修正オイラー法によれば、オイラー法に比較して発散の程度を小さくすることができるが、依然として発散した動解析となってしまい、精度が低い動解析結果となってしまっていた。したがって、オイラー法や修正オイラー法を用いた動解析では、解析微小時間Δｔを極端に小さくし、多大な時間をかけて精度を確保せざるを得なかった。これに対し、上述した式（１）〜（３）を繰り返し行う動解析処理では、解析微小時間Δｔが短くても精度の高い動解析結果を得ることができる。【００３９】  この発明による動解析処理、従来のオイラー法、修正オイラー法の概念を図化すると図３〜図５に示すようになる。図３は、この発明の実施の形態１による動解析処理の概念を示す図であり、図４は、従来のオイラー法による動解析処理の概念を示す図であり、図５は、従来の修正オイラー法による動解析処理の概念を示す図である。図３〜図５において、この発明による動解析処理では、図３に示すように、時刻ｔにおける変位および速度から加速度を求め、この時刻ｔにおける速度と加速度とを用いて解析微小時間Δｔ経過時刻（ｔ＋Δｔ）の速度を求める。その後、この時刻（ｔ＋Δｔ）の速度と時刻ｔの変位を用いて時刻（ｔ＋Δｔ）の変位を求めるようにしており、特に、時刻（ｔ＋Δｔ）時の速度を用いて時刻（ｔ＋Δｔ）の変位を求めることが特徴である。一方、オイラー法では、図４に示すように、時刻ｔにおける速度と時刻ｔにおける変位を用いて、時刻（ｔ＋Δｔ）の変位を求めるようにしている。さらに、従来の修正オイラー法では、図５に示すように、時刻ｔと時刻（ｔ＋Δｔ）との間の平均速度（１／２・（ｖ（ｔ）＋ｖ（ｔ＋Δｔ）））を用いて、時刻（ｔ＋Δｔ）における変位を求めるようにしている。【００４０】  ここで、具体的な単振動の解析方程式に対する処理結果について説明する。図６は、解析方程式の対象である、ばね−質点系の概念構成を示す模式図である。図６に示したばね−質点系では、単振動の運動状態を呈する。図６に示すばね−質点系の運動方程式は、次式（６）に示すように、
2006103996.tif 000005
となる。なお、このばね−質点系では、ばね定数ｋ＝４π２（Ｎ／ｍ）、質量ｍ＝１．０（ｋｇ）、初期変位ｘ（０）＝２．０（ｍ）、初期速度ｖ（０）＝０．０（ｍ）である。この運動方程式の解析方程式は、次式（７）に示すように、
2006103996.tif 000006
となり、パラメータを代入すると、次式（８）に示す解析方程式となる。すなわち、
2006103996.tif 000007
となる。なお、このばね−質点系の代数解は、次式（９）に示すように、
2006103996.tif 000008
として表される。解析微小時間Δｔ＝０．０１（ｓ）としている。【００４１】  図７は、図６に示したばね−質点系の解析結果を示す図である。この動解析は、解析微小時間Δｔ＝０．０１（ｓ）で行っている。図７に示すように、オイラー法および修正オイラー法では、時間の経過とともに発散しているのに対し、この発明による動解析処理では、発散もなくほぼ代数解に一致した結果を得ることができている。この結果、この発明の実施の形態１によれば、オイラー法や修正オイラー法に比して短時間に精度の高い動解析を行うことができる。換言すれば、この実施の形態１による動解析処理は、同じ解析時間では、オイラー法や修正オイラー法に比して格段に精度の高い動解析を行うことができる。【００４２】  この実施の形態１による動解析処理は、複雑なプログラムではないため、たとえば表計算プログラムを用いたソフトウェアによって容易に実現することができる。図８は、図７に示したばね−質点系を動解析する表計算プログラムの表示画面の一例を示す図である。図８において、まず時間列のセルＣ１１には初期時刻ｔ＝０．００を設定し、その後のセルＣ２１〜Ｃｎ１には、微小解析時間Δｔ＝０．０１が加算されるようにしている。基準時刻ｔ時点を示す区間初期の変位列のセルＣ１２には、初期変位ｘ（０）＝２．００が設定され、区間初期の速度列のセルＣ１３には、初期速度ｖ（０）＝０．００が設定される。加速度列のセルＣ１４には、セルＣ１２，Ｃ１３に示された変位および速度をもとに加速度が演算され、その結果が表示される。基準時刻ｔから微小解析時間Δｔ経過後を示す区間終期の速度列のセルＣ１５には、式（２）に対応した演算結果が表示され、区間終期の速度列のセルＣ１６には、式（３）に対応した演算結果が表示される。そして、このセルＣ１５，Ｃ１６に示された結果は、次の微小解析時間Δｔ経過後の時間行のセルＣ２２，Ｃ２３には、次の区間初期の変位および速度としてコピーされる。さらに、セルＣ２２，Ｃ２３に示された変位および速度をもとに、加速度が求められ、結果がセルＣ２４に表示され、これらセルＣ２２，Ｃ２３，Ｃ２４に示された変位、速度、加速度をもとに次の微小解析時間Δｔ経過後の速度および変位が求められ、セルＣ２５，Ｃ２６に表示される。以後、同様にして、区間終期における速度および変位が次の微小解析時間Δｔ経過後の区間初期の速度および変位として用いられ、この時点における加速度が求められて表示されるとともに、次の微小解析時間Δｔ経過後の区間終期における速度および変位が求められて表示される。この処理は、時間列で設定した繰り返し数すなわち行数分、繰り返される。また、これらの結果は、表計算ソフトウェアに通常付加されたグラフ化処理によってグラフ化することができる。このグラフ化処理は、出力処理部１３が行うことになる。【００４３】  図９は、図６に示したばね−質点系に対する動解析処理を表計算プログラムによって行った一例を示す図である。なお、ここでは、初期変位ｘ（０）＝１．０００（ｍ）としている。また、図１０は、図９に示した結果をグラフ化したものであり、横軸に時間をとり、縦軸に変位をとっている。図１０に示すように、代数解と数値解析（動解析）結果は��ぼ一致しており、その誤差は、±０．０４２ｍであった。また、図１０では、１０（ｓ）までの動解析であったが、さらに１００（ｓ）までの動解析を行っても、誤差は、±０．１３４ｍであった。【００４４】  なお、図１１は、同じばね−質点系に対する動解析処理をオイラー法によって行った一例を示す図であり、図１２は、図１１に示した動解析結果をグラフ化したものである。さらに、図１３は、同じばね−質点系に対する動解析を修正オイラー法によって行った結果をグラフ化したものである。図１１および図１２に示すように、同じ微小解析時間によるオイラー法によって動解析した場合、時間の経過とともに動解析結果と代数解との誤差が大きくなり、発散している。この場合における誤差は、１０（ｓ）で±６ｍ程度まで広がっている。さらに、修正オイラー法によっても、図１３に示すように、時間の経過とともに動解析結果と代数解との誤差が大きくなり、発散している。この場合における誤差は、オイラー法よりも小さいが、１０（ｓ）で±２ｍ程度となっている。【００４５】（実施の形態２）
つぎに、実施の形態２について説明する。上述した実施の形態１では、現時点の速度などの一階時間微分量および現時点の加速度などの二階微分量をもとに次の時点の一階時間微分量を求め、さらにこの次の時点の一階時間微分量および現時点の変位などの物理量を用いて次の時点の物理量を求めるようにしていたが、この実施の形態２では、現時点の一階時間微分量、現時点の二階時間微分量および現時点の物理量をもとに、一度で次の時点の物理量を求めるようにしている。ただし、次の時点の二階時間微分量を算出する際に用いる次の時点の一階時間微分量は、現時点において求めておく。【００４６】  図１４は、この発明の実施の形態２である数値解析装置の一例である動解析装置の概要構成を示すブロック図である。図１４に示すように、この動解析装置は、物理量算出モジュール２４に替えて、物理量算出モジュール２４ａを設けている。その他の構成は、実施の形態１と同じであり、同一構成部分には同一符号を付している。【００４７】  二階時間微分量算出モジュール２２は、一階時間微分量算出モジュール２３が求めておいた速度ｖ（ｔ）と、物理量算出モジュール２４が求めておいた基準時刻ｔにおける変位ｘ（ｔ）とを用いて、加速度α（ｔ）を求め、物理量算出モジュール２４ａは、この加速度α（ｔ）と、速度ｖ（ｔ）と、基準時刻ｔにおける変位ｘ（ｔ）とを用いて、直接、解析微小時間Δｔ経過後の変位ｘ（ｔ＋Δｔ）を求めている。これは、式（３）に式（２）を代入した次式（１０）によって直接、変位ｘ（ｔ＋Δｔ）を求めることができることを意味する。すなわち、
2006103996.tif 000009
によって求めることができる。【００４８】  図１５は、図１４に示した動解析装置による動解析処理手順を示すフローチャートである。図１５に示した動解析処理手順では、実施の形態１におけるステップ１０４の処理をステップＳ２０４の処理に置き換えた処理を行っている。その他のステップＳ２０１〜Ｓ２０３，Ｓ２０５〜Ｓ２０７は、実施の形態１に示したステップＳ１０１〜Ｓ１０３，Ｓ１０５〜Ｓ１０７と同じ処理である。【００４９】  すなわち、ステップＳ２０４では、ステップＳ２０２によって算出された加速度α（ｔ）に解析微小時間Δｔの２乗値を乗算し、さらに基準時刻ｔの速度ｖ（ｔ）に解析微小時間Δｔを乗算し、これら乗算した各値と基準時刻ｔの変位ｘ（ｔ）とを加算し、この加算値を解析微小時間Δｔ経過後の変位ｘ（ｔ＋Δｔ）として算出するようにしている。なお、ステップＳ２０３で基準時刻ｖ（ｔ）の算出処理を行うのは、実施の形態１と異なり、次の時点の変位を算出する場合に速度ｖ（ｔ）が求められていないためである。【００５０】  これによって、この実施の形態２では、実施の形態１とは異なり、現時点の速度、加速度および変位をもとに、一度に解析微小時間Δｔ経過後の速度ｖ（ｔ＋Δｔ）を算出することができる。【００５１】  なお、この発明の実施の形態２による動解析処理の概念を図化すると図１６に示すようになる。図１６は、この発明の実施の形態２による動解析処理の概念を示す図である。図１６において、この発明の実施の形態２による動解析処理では、時刻ｔにおける変位および速度から加速度を求め、この時刻ｔにおける速度と加速度とを用いて解析微小時間Δｔ経過時刻（ｔ＋Δｔ）の変位を求めるようにしている。【００５２】  ここで、上述した実施の形態１，２では、解析方程式α（ｔ）＝ｆ（ｘ（ｔ），ｖ（ｔ））として説明したが、速度項ｖ（ｔ）を含まない解析方程式α（ｔ）＝ｆ（ｘ（ｔ））であってもよい。また、解析方程式α（ｔ）には、時間項を含んでも良い。すなわち、α（ｔ）＝ｆ（ｘ（ｔ），ｖ（ｔ），ｔ）あるいはα（ｔ）＝ｆ（ｘ（ｔ），ｔ）であってもよい。さらに、解析対象によっては、複数の解析方程式が連立方程式であってもよい。すなわち、解析方程式は、多元方程式であってもよい。この多元方程式は、一般化すると、αｉ（ｔ）＝ｆ（ｘ１（ｔ），…，ｘｎ（ｔ），ｖ１（ｔ），…，ｖｎ（ｔ），ｔ）とするｎ元の連立微分方程式として表すことができる。ただし、ｉ＝１〜ｎである。【００５３】  なお、上述した実施の形態１，２では、数値解析装置の一例として動解析装置をあげて説明したが、時間を変数量とする動解析に限らず、一般の数値解析にも適用することができる。すなわち、時間ｔを変数量とするのではなく、温度Ｔや距離Ｌを変数量とした数値解析を行うことができる。たとえば、時刻ｔおよび解析微小時間Δｔに替えて、温度Ｔおよび解析微小温度ΔＴとすれば、温度変化に伴う数値解析を行うことができる。【００５４】  つぎに、上述した実施の形態１，２の応用例として、ばね質点ダッシュポット系の動解析処理、運搬機器の動解析処理、動力伝達機等で連結された複数の回転体の動解析処理、調速機の動解析処理、旋削機等のバイトの動解析処理、クランクを有する往復機械の動解析処理、摩擦を利用したブレーキの動解析処理、剛体リンクで接続されたマニピュレータ状機器等の動解析処理、重力等を相互に及ぼし合う粒子あるいは天体の動解析処理、太陽・地球・月の三体問題の動解析処理、遠心ポンプを有するシステムのサージング特性解析処理、およびインダクタンス素子、コンデンサ素子、抵抗素子を含む電気・電子回路に生じる電流あるいは電荷の動解析処理について説明する。【００５５】（応用例１）
図１７は、動解析対象であるばね質点ダッシュポット系の概要構成を示す模式図である。このばね質点ダッシュポット系に固有の運動方程式は、次式（１１）で示される。
2006103996.tif 000010
ただし、ｍ：質点の質量（１５００．００）
ｃ：ダッシュポットの粘性減衰係数（４２４２．６４）
ｋ：バネ係数（１２０００）
ｆ：外力の振幅（１０００．０）
Ｔ：外力の周期（１．００）
ｔ：基準時刻
ｘ：質点の変位
ｖ：質点の速度
α：質点の加速度
ｔ０：解析開始時刻
ｘ０：質点の初期変位（０．１０）
ｖ０：質点の初期速度（１．００）
α０：質点の初期加速度
であり、これらはＳＩ単位系で示している。【００５６】  この場合、解析方程式は、式（１１）から、次式（１２）で示される。すなわち、
2006103996.tif 000011
によって示される。したがって、式（１２）で示したα（ｔ）に、各種パラメータおよび、ｔ＝ｔ０、ｘ（ｔ０）＝ｘ０、ｖ（ｔ０）＝ｖ０を代入してα（ｔ０）を求め、式（２）および式（３）で示す演算を行って解析微小時間Δｔ（＝０．０１（ｓ））経過後のｖ（ｔ＋Δｔ）およびｘ（ｔ＋Δｔ）を求める。このｖ（ｔ＋Δｔ）およびｘ（ｔ＋Δｔ）を用いて次の解析微小時間Δｔ経過後の加速度α（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。【００５７】  図１８は、図１７に示した質点ダッシュポット系の動解析結果を示す図であり、横軸に時間をとり、縦軸に変位をとっている。この動解析結果は、図１８に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００５８】（応用例２）
図１９は、動解析対象である、単振動と回転運動とが伴った運搬機の概要構成を示す模式図である。この運搬機に固有の運動方程式は、次式（１３），（１４）で示される。
2006103996.tif 000012
ただし、ｍ：運搬機の質量（１５００．００）
Ｊ：運搬機の慣性モーメント（３１２５）
ｋ１：バネ定数１（３０００）
ｋ２：バネ係数２（３０００）
Ｌ１：重心からの距離１（２．００）
Ｌ２：重心からの距離２（２．００）
ｘ：運搬機の変位（初期値ｘ０＝０．１０）
ｖ：運搬機の速度（初期値ｖ０＝０．００）
α：運搬機の加速度
θ：運搬機の角変位（初期値θ０＝０．２６）
ω：運搬機の角速度（初期値ω０＝０．００）
ｄ２θ／ｄｔ２：運搬機の角加速度
であり、これらはＳＩ単位系で示している。【００５９】  この場合、解析方程式は、式（１３），（１４）から、次式（１５），（１６）で示される。すなわち、
2006103996.tif 000013
によって示される。したがって、各式（１５），（１６）に、各種パラメータおよび、初期変位ｖ（ｔ０）＝０．００、初期角変位θ０（ｔ０）＝０．２５を代入して、それぞれα（ｔ０）および（ｄ２θ／ｄｔ０２）を求め、式（１０）で示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後のｘ（ｔ＋Δｔ）およびθ（ｔ＋Δｔ）を求める。このｘ（ｔ＋Δｔ）およびθ（ｔ＋Δｔ）を用いて次の解析微小時間Δｔ経過後の加速度α（ｔ＋Δｔ）および（ｄ２θ（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。なお、実施の形態１のように、式（２）および式（３）を用いる場合、それぞれ運搬機の速度ｖ（ｔ＋Δｔ）および運搬機の角速度ω（ｔ＋Δｔ）を式（３）によって求める。この場合、運搬機の初期速度ｖ（ｔ０）＝０．００および運搬機の角速度ω（ｔ０）＝０．００とする。【００６０】  図２０は、図１９に示した運搬機の変位の動解析結果を示す図であり、図２１は、図１９に示した運搬機の角変位の動解析結果を示す図である。それぞれの動解析結果は、図２０および図２１に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００６１】（応用例３）
図２２は、動解析対象である、ねじり軸、継手、ギア等の動力伝達器で連結された複数の回転体の概要構成を示す模式図である。この解析対象の運動方程式は、次式（１７），（１８）で示される。
2006103996.tif 000014
ただし、Ｉ１：回転体１の慣性モーメント（９５００）
Ｉ２：回転体２の慣性モーメント（５３９００）
ｋ１２：捩りバネ定数（９８９８００００）
θ１：回転体１の角変位（初期値＝０．０００）
θ２：回転体２の角変位（初期値＝０．０２６）
ω１：回転体１の角速度（初期値＝０．０００）
ω２：回転体２の角速度（初期値＝０．０００）
であり、これらはＳＩ単位系で示している。【００６２】  この場合、解析方程式は、式（１７），（１８）から、次式（１９），（２０）で示される。すなわち、
2006103996.tif 000015
によって示される。したがって、各式（１９），（２０）に、各種パラメータおよび、初期角変位θ１０＝０．０００、初期角変位θ２０＝０．０２６を代入して、それぞれ（ｄ２θ１／ｄｔ２）および（ｄ２θ２／ｄｔ２）を求め、式（１０）で示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後のθ１（ｔ＋Δｔ）およびθ２（ｔ＋Δｔ）を求める。このθ１（ｔ＋Δｔ）およびθ２（ｔ＋Δｔ）を用いて次の解析微小時間Δｔ経過後の（ｄ２θ１（ｔ＋Δｔ）／ｄｔ２）および（ｄ２θ２（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。なお、実施の形態１のように、式（２）および式（３）を用いる場合、それぞれの角速度ω１（ｔ＋Δｔ），ω２（ｔ＋Δｔ）を式（３）によって求める。この場合、初期角速度ω１（ｔ＋Δｔ），ω２（ｔ＋Δｔ）＝０．０００とする。【００６３】  図２３および図２４は、図２２に示した動力伝達器で結合された各回転体の角変位の動解析結果を示す図である。それぞれの動解析結果は、図２３および図２４に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００６４】（応用例４）
図２５は、動解析対象である、調速機のすべり子およびタービンの概要構成を示す模式図である。この調速機に固有の運動方程式は、次式（２１），（２２）で示される。
2006103996.tif 000016
ただし、ｍｇｖ：調速機のすべり子の等価質量（１０．００）
ｃｇｖ：調速機の粘性減衰係数（５０．６６）
ｋｇｖ：調速機のバネ定数（１５７９．１４）
Ａｇｖ：タービンの回転速度が１ｒａｄ／ｓ増すことによって調速機のすべり子に働く上向きの力（１０．００）
Ｊｇｖ：タービンの等価慣性モーメント（２．５０）
Ｂｇｖ：調速機のすべり子が単位長下がるために生じるタービンのトルク増加（１０００．００）
ｘｇｖ：調速機のすべり子の変位（初期値＝０．００）
ｖｇｖ：調速機のすべり子の速度（初期値＝０．００）
αｇｖ：調速機のすべり子の加速度
φｇｖ：タービンの角変位（初期値＝０．００）
ωｇｖ：タービンの角速度（初期値＝０．２０）
であり、これらはＳＩ単位系で示している。【００６５】  この場合、解析方程式は、式（２１），（２２）から、次式（２３），（２４）で示される。すなわち、
2006103996.tif 000017
によって示される。したがって、各式（２３），（２４）に、各種パラメータおよび、初期変位ｘｇｖ（ｔ０）＝０．００、初期角変位φｇｖ（ｔ０）＝０．００、初期速度ｖｇｖ（ｔ０）＝０．００、初期角速度ωｇｖ（ｔ０）＝０．２０を代入して、加速度αｇｖおよび角加速度（ｄ２φｇｖ／ｄｔ２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０００５（ｓ））経過後の速度ｖｇｖ（ｔ＋Δｔ）およびωｇｖ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の変位ｘｇｖ（ｔ＋Δｔ）および角変位φｇｖ（ｔ＋Δｔ）を求める。これら求められたｖｇｖ（ｔ＋Δｔ）、ωｇｖ（ｔ＋Δｔ）、ｘｇｖ（ｔ＋Δｔ）、およびφｇｖ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０００５（ｓ））経過後の加速度αｇｖ（ｔ＋Δｔ）および角加速度（ｄ２φｇｖ（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。【００６６】  図２６は、図２５に示した調速機のすべり子の変位の動解析結果を示す図であり、図２７は、図２５に示したタービンの角変位の動解析結果を示す図である。それぞれの動解析結果は、図２６および図２７に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００６７】（応用例５）
図２８は、動解析対象である、旋削機などのバイトの概要構成を示す模式図である。この旋削機に固有の運動方程式は、次式（２５）で示される。
2006103996.tif 000018
ただし、ｍｂ：バイトの質量（１０．００）
ｃｂ：粘性減衰係数（２０．００）
ｋｂ：バネ定数（１０００．００）
Ｖ：切削速度（２．００）
Ｐｂ：背分力（１００．００）
ｘｂ：バイトの変位（初期値＝０．００）
ｖｂ：バイトの速度（初期値＝１．００）
αｂ：バイトの加速度
μ（Ｖ−ｖｂ）Ｐｂ：０．００７５（Ｖ−ｖｂ）２−０．０５５（Ｖ−ｖｂ）＋０．４００
であり、これらはＳＩ単位系で示している。【００６８】  この場合、解析方程式は、式（２５）から、次式（２６）で示される。すなわち、
2006103996.tif 000019
によって示される。したがって、式（２６）に、各種パラメータおよび、初期変位ｘｂ（ｔ０）＝０．００、初期速度ｖｂ（ｔ０）＝１．００を代入して、加速度αｂを求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０２（ｓ））経過後の速度ｖｂ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０２（ｓ））経過後の変位ｘｂ（ｔ＋Δｔ）を求める。これら求められたｖｂ（ｔ＋Δｔ）、ｘｂ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０２（ｓ））経過後の加速度αｂ（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。【００６９】  図２９は、図２８に示した旋削機などのバイト変位の動解析結果を示す図である。この動解析結果は、図２９に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００７０】（応用例６）
図３０は、動解析対象である、クランクを有する往復機械の概要構成を示す模式図である。このクランクを有する往復機械に固有の運動方程式は、次式（２７）で示される。
2006103996.tif 000020
ただし、ｍｐｉ：往復質量（１５．００）
Ｉｃｒ：クランク軸周りの全回転質量の慣性モーメント（７．７３）
Ｒｃｒ：クランク軸の腕の長さ（０．１３）
ε：Ｒｃｒ／Ｌｃｏｎ
Ｌｃｏｎ：接続の長さ（０．２５）
Ｐｐｉ：外力（０．００）
Ｔｒｃｒ：抵抗トルク
θｃｒ：クランクの角変位（初期値＝０．２６）
ｄθｃｒ／ｄｔ：クランクの角速度（初期値＝１２．５７）
であり、これらはＳＩ単位系で示している。【００７１】  この場合、解析方程式は、式（２７）から、次式（２８）で示される。すなわち、
2006103996.tif 000021
2006103996.tif 000022
によって示される。したがって、式（２８）に、各種パラメータおよび、初期角変位θｃｒ（ｔ０）＝０．２６、初期角速度（ｄθｃｒ／ｄｔ０）＝１２．５７を代入して、角加速度（ｄ２θｃｒ／ｄｔ０２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の角速度（ｄθｃｒ（ｔ＋Δｔ）／ｄｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の角変位θｃｒ（ｔ＋Δｔ）を求める。これら求められた（ｄθｃｒ（ｔ＋Δｔ）／ｄｔ）、θｃｒ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．００１（ｓ））経過後の角加速度（ｄ２θｃｒ（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。【００７２】  図３１は、図３０に示したクランクを有する往復機械の動解析結果を示す図である。この動解析結果は、図３１に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００７３】  なお、ピストンの変位は、機構的に次式（２９）によって示すように、クランク角度などをもとに求めることができる。すなわち、
2006103996.tif 000023
ただし、ｘｐｉ：ピストンの変位
Ｒｃｒ：クランク軸の腕の長さ（０．１３）
ε：Ｒｃｒ／Ｌｃｏｎ
Ｌｃｏｎ：接続の長さ（０．２５）
θｃｒ：クランクの角変位
であり、これらはＳＩ単位系で示している。図３２は、この式（２９）をもとに求めた往復機械のピストンの動解析結果を示す図である。この動解析結果は、図３２に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００７４】（応用例７）
図３３は、動解析対象である、ブレーキディスクの概要構成を示す模式図である。このブレーキディスクに固有の運動方程式は、次式（３０）で示される。
2006103996.tif 000024
ただし、ＪＢ：ブレーキディスクの慣性モーメント（３．００）
ｆＢ：押付力（５００．００）
ｕＢ：摩擦係数（０．３０）
θＢ：ブレーキディスクの角変位（初期値＝０．００）
ｄθＢ（ｔ）／ｄｔ：ブレーキディスクの角速度（初期値＝１３９．６３）
であり、これらはＳＩ単位系で示している。【００７５】  この場合、解析方程式は、式（３０）から、次式（３１）で示される。すなわち、
2006103996.tif 000025
によって示される。したがって、式（３０）に、各種パラメータおよび、初期角変位θＢ（ｔ０）＝０．００、初期角速度（ｄθＢ（ｔ）／ｄｔ０）＝１３９．６３を代入して、角加速度（ｄ２θＢ／ｄｔ０２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角速度（ｄθＢ（ｔ＋Δｔ）／ｄｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角変位θＢ（ｔ＋Δｔ）を求める。これら求められた（ｄθＢ（ｔ＋Δｔ）／ｄｔ）、θＢ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の角加速度（ｄ２θＢ（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。【００７６】  図３４は、図３３に示したブレーキディスクの動解析結果を示す図である。この動解析結果は、図３４に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００７７】（応用例８）
図３５は、動解析対象である、剛体リンクで接続されたマニピュレータ状機器の概要構成を示す模式図である。この剛体リンクで接続されたマニピュレータ状機器に固有の運動方程式は、次式（３２），（３３）で示される。
2006103996.tif 000026
ただし、Ｌｌｉｎｋ：リンク１，２の長さ（５．０００）
ｇ：重力加速度（９．８０７）
θｌｉｎｋ１：リンク１の角変位（初期値＝０．０８７）
ｄθｌｉｎｋ１／ｄｔ：リンク１の角速度（初期値＝０．０００）
θｌｉｎｋ２：リンク２の角変位（初期値＝−０．０８７）
ｄθｌｉｎｋ２／ｄｔ：リンク２の角速度（初期値＝０．０００）
であり、これらはＳＩ単位系で示している。【００７８】  この場合、解析方程式は、式（３２），（３３）から、次式（３４），（３５）で示される。すなわち、
2006103996.tif 000027
によって示される。したがって、式（３４），（３５）に、各種パラメータおよび、初期角変位θｌｉｎｋ１（ｔ０）＝０．０８７、初期角変位θｌｉｎｋ２（ｔ０）＝−０．０８７、初期角速度（ｄθｌｉｎｋ１（ｔ）／ｄｔ０）＝０．０００、初期角速度（ｄθｌｉｎｋ２（ｔ）／ｄｔ０）＝０．０００を代入して、角加速度（ｄ２θｌｉｎｋ１／ｄｔ０２）、（ｄ２θｌｉｎｋ２／ｄｔ０２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角速度（ｄθｌｉｎｋ１（ｔ＋Δｔ）／ｄｔ）、（ｄθｌｉｎｋ２（ｔ＋Δｔ）／ｄｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角変位θｌｉｎｋ１（ｔ＋Δｔ）、θｌｉｎｋ２（ｔ＋Δｔ）を求める。これら求められた（ｄθｌｉｎｋ１（ｔ＋Δｔ）／ｄｔ）、（ｄθｌｉｎｋ２（ｔ＋Δｔ）／ｄｔ）、θｌｉｎｋ（ｔ＋Δｔ）、θｌｉｎｋ２（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の角加速度（ｄ２θｌｉｎｋ１（ｔ＋Δｔ）／ｄｔ２）、（ｄ２θｌｉｎｋ２（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。【００７９】  図３６および図３７は、図３５に示した剛体リンクで接続されたマニピュレータ状機器の動解析結果を示す図である。この動解析結果は、図３６および図３７に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００８０】  なお、式（３２），（３３）において、角変位θｌｉｎｋ１，θｌｉｎｋ２が微小な場合、式（３４），（３５）は、次式（３６），（３７）に示すように近似することができる。すなわち、
2006103996.tif 000028
として近似できる。この式をもとに動解析を行っても図３６，図３７とほぼ同じ結果を得ることができる。この場合、実施の形態２と同様な処理を適用することができる。【００８１】（応用例９）
図３８は、動解析対象である、重力等を相互に及ぼし合う粒子あるいは天体の概要構成を示す模式図である。この重力等を相互に及ぼし合う粒子あるいは天体に固有の運動方程式は、次式（３８）〜（４１）で示される。
2006103996.tif 000029
ただし、ｍ１ｇ：地球の質量（５．９７×１０２４）
ｍ２ｇ：月の質量（７．３５×１０２２）
Ｇ：重力定数（６．６７×１０−１１）
αｘ１：地球のｘ方向加速度
αｙ１：地球のｙ方向加速度
αｘ２：月のｘ方向加速度
αｙ２：月のｙ方向加速度
ｒ：地球と月の距離
であり、これらはＳＩ単位系で示している。【００８２】  なお、式（３８）〜（４１）に示したｒ，ｃｏｓθ，ｓｉｎθは、次式（４２）〜（４４）によって定義される。
2006103996.tif 000030
ただし、ｘ１：地球のｘ方向変位（初期値＝０．００×１０００）
ｘ２：地球のｙ方向変位（初期値＝０．００×１０００）
ｙ１：月のｘ方向変位（初期値＝３．８４×０８）
ｙ２：月のｙ方向変位（初期値＝０．００×１０００）
ｖｘ１：地球のｘ方向速度（初期値＝０．００×１０００）
ｖｙ１：地球のｙ方向速度（初期値＝０．００×１０００）
ｖｘ２：月のｘ方向速度（初期値＝０．００×１０００）
ｖｙ２：月のｙ方向速度（初期値＝１．０４×１００３）【００８３】  この場合、解析方程式は、式（３８）〜（４１）から、次式（４５）〜（４８）で示される。すなわち、
2006103996.tif 000031
によって示される。したがって、式（４５）〜（４８）に、各種パラメータおよび、地球のｘ方向初期変位ｘ１（ｔ０）＝０．００×１０００、地球のｙ方向初期変位ｙ１（ｔ０）＝０．００×１０００、月のｘ方向初期変位ｘ２（ｔ０）＝３．８４×１００８、月のｙ方向初期変位ｙ２（ｔ０）＝０．００×１０００を代入して、加速度αｘ１，αｙ１，αｘ２，αｙ２を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の速度ｖｘ１（ｔ＋Δｔ），ｖｙ１（ｔ＋Δｔ），ｖｘ２（ｔ＋Δｔ），ｖｙ２（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の変位ｘ１（ｔ＋Δｔ），ｘ２（ｔ＋Δｔ），ｙ１（ｔ＋Δｔ），ｙ２（ｔ＋Δｔ）を求める。これら求められたｖｘ１（ｔ＋Δｔ），ｖｙ１（ｔ＋Δｔ），ｖｘ２（ｔ＋Δｔ），ｖｙ２（ｔ＋Δｔ），ｘ１（ｔ＋Δｔ），ｘ２（ｔ＋Δｔ），ｙ１（ｔ＋Δｔ），ｙ２（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の加速度αｘ１（ｔ＋Δｔ），αｙ１（ｔ＋Δｔ），αｘ２（ｔ＋Δｔ），αｙ２（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。【００８４】  図３９は、図３８に示した地球のｘ方向変位の動解析結果を示す図である。この動解析結果は、図３９に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００８５】（応用例１０）
この発明では、三体問題についても精度の高い数値解析処理を短時間に行うことができる。三体問題とは、文字通り、３つの物体があるとき、万有引力を用いてそれぞれの物体の振る舞いを求める問題であり、ここでは、太陽、地球、および月のそれぞれに対する動解析を行う。【００８６】  図４０は、三体問題の動解析対象である太陽、地球、および月の位置関係を示す模式図である。なお、太陽１０１，地球１０２，月１０３はそれぞれｘ−ｙ平面（地球の公転面）の同一平面上で運動するものとして解析する。【００８７】  まず、太陽１０１の運動方程式は、
2006103996.tif 000032
として表すことができる。ただし、Ｍ１は太陽１０１の質量であり、ａ１は太陽１０１の加速度である。なお、ａ１はベクトルａ１（ａｘ１，ａｙ１）である。
地球１０２による万有引力ｆ１２は、
2006103996.tif 000033
で表される。Ｇは、万有引力定数を示し、Ｍ２は地球１０２の質量であり、Ｐ１およびＰ２はそれぞれ太陽１０１および地球１０２の位置ベクトルＰ１（ｘ１，ｙ１），Ｐ２（ｘ２，ｙ２）である。また、万有引力ｆ１２は、ベクトルである。なお、｛Ｇ・Ｍ２・Ｍ１／（｜Ｐ２−Ｐ１｜２）｝は、万有引力の絶対値を意味し、｛（Ｐ２−Ｐ１）／｜Ｐ２−Ｐ１｜｝は、太陽１０１から地球１０２への単位ベクトルを意味する。
式（５０）は、次式に整理される。
2006103996.tif 000034
同様にして、月１０３による万有引力ｆ１３は、
2006103996.tif 000035
で表される。ここで、Ｍ３は月の質量であり、Ｐ３は月の位置ベクトルＰ３（ｘ３，ｙ３）である。また、万有引力ｆ１３は、ベクトルである。
ここで、
2006103996.tif 000036
であるので、式（４９），（５１）〜（５３）からＦ１，ｆ１２，ｆ１３を消去すると、
2006103996.tif 000037
となる。ここで、｛Ｇ・Ｍ２／（｜Ｐ２−Ｐ１｜３）｝＊（Ｐ２−Ｐ１）は地球１０２の引力による加速度分であり、｛Ｇ・Ｍ３／（｜Ｐ３−Ｐ１｜３）｝＊（Ｐ３−Ｐ１）は、月１０３の引力による加速度分である。【００８８】  同様にして地球の加速度ａ２は、
2006103996.tif 000038
となる。ここで、｛Ｇ・Ｍ３／（｜Ｐ３−Ｐ２｜３）｝＊（Ｐ３−Ｐ２）は月１０３の引力による加速度分であり、｛Ｇ・Ｍ１／（｜Ｐ１−Ｐ２｜３）｝＊（Ｐ１−Ｐ２）は、太陽１０１の引力による加速度分である。なお、ａ２はベクトルａ２（ａｘ２，ａｙ２）である。
さらに同様にして月の加速度ａ３は、
2006103996.tif 000039
となる。ここで、｛Ｇ・Ｍ１／（｜Ｐ１−Ｐ３｜３）｝＊（Ｐ１−Ｐ３）は太陽１０１の引力による加速度分であり、｛Ｇ・Ｍ２／（｜Ｐ２−Ｐ３｜３）｝＊（Ｐ２−Ｐ３）は、地球１０２の引力による加速度分である。なお、ａ３はベクトルａ３（ａｘ３，ａｙ３）である。【００８９】  上述した式（５４），（５５），（５６）の３つの連立微分方程式を解析微小時間Δｔ＝１時間ごとに変化させて３体問題の数値解析を行うことができるが、各連立微分方程式は、それぞれ２次元であり、結果的に６元の連立微分方程式を同時に解析することになる。ここで、初期値について説明する。まず、初期状態としては太陽１０１，地球１０２，月１０３が一直線に並んだ瞬間を時刻ｔ＝０とし、太陽１０１を原点とし、地球１０２および月１０３の方向を＋ｘ軸方向にとり、地球１０２および月１０３の速度ベクトルを＋ｙ軸方向にとる。また、
太陽１０１と地球１０２との間の初期距離  ：１．４９６×１０１１（ｍ）
地球１０２と月１０３との間の初期距離  ：３．８４４×１０８（ｍ）
地球１０２の初期速度                ：２．９７８×１０４（ｍ／ｓ）
月１０３の初期周期                ：２７．３２１７（日）
月１０３の地球１０２に対する初期相対速度  ：１０２３（ｍ／ｓ）
太陽１０１の質量Ｍ１              ：１．９８９×１０３０（ｋｇ）
地球１０２の質量Ｍ２              ：５．９７４×１０２４（ｋｇ）
月１０３の質量Ｍ３                ：７．３４７×１０２２（ｋｇ）
万有引力定数Ｇ              ：６．６７２５９×１０−１１（Ｎ・ｍ２／ｋｇ２）
である。【００９０】  この三体問題の数値解析は、まず上述した初期値を式（５４），（５５），（５６）に代入して、加速度ａｘ１，ａｙ１，ａｘ２，ａｙ２，ａｘ３，ａｙ３を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝３６００（ｓ））経過後の速度ｖｘ１（ｔ＋Δｔ），ｖｙ１（ｔ＋Δｔ），ｖｘ２（ｔ＋Δｔ），ｖｙ２（ｔ＋Δｔ），ｖｘ３（ｔ＋Δｔ），ｖｙ３（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝３６００（ｓ））経過後の変位ｘ１（ｔ＋Δｔ），ｙ１（ｔ＋Δｔ），ｘ２（ｔ＋Δｔ），ｙ２（ｔ＋Δｔ），ｘ３（ｔ＋Δｔ），ｙ３（ｔ＋Δｔ）を求める。これら求められたｖｘ１（ｔ＋Δｔ），ｖｙ１（ｔ＋Δｔ），ｖｘ２（ｔ＋Δｔ），ｖｙ２（ｔ＋Δｔ），ｖｘ３（ｔ＋Δｔ），ｖｙ３（ｔ＋Δｔ），ｘ１（ｔ＋Δｔ），ｙ１（ｔ＋Δｔ），ｘ２（ｔ＋Δｔ），ｙ２（ｔ＋Δｔ），ｘ３（ｔ＋Δｔ），ｙ３（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝３６００（ｓ））経過後の加速度αｘ１（ｔ＋Δｔ），αｙ１（ｔ＋Δｔ），αｘ２（ｔ＋Δｔ），αｙ２（ｔ＋Δｔ），αｘ３（ｔ＋Δｔ），αｙ３（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。【００９１】  図４１は、上述した３体問題における地球を中心とした月の軌道の動解析結果を示す図である。図４１に示すように、この動解析結果は、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００９２】  なお、図４２は、従来のオイラー法を用いて地球を中心とした月の軌道の動解析結果を示す図である。この図４２に示した動解析結果は、この応用例１０と同じ解析微小時間Δｔ＝３６００ｓで行った結果であり、月の軌道が発散状態となり、高精度な数値解析結果を得ることができていない。【００９３】（応用例１１）
図４３は、動解析対象である、遠心ポンプを有するシステムのサージタンクの概要構成を示す模式図である。このサージタンクの水位変動に固有の運動方程式は、次式（５７）で示される。
2006103996.tif 000040
ただし、ｇ：重力加速度（９．８０７）
ｌｓｕｒｇｅ：管路の長さ（１０．０００）
Ａｓｕｒｇｅ：管路断面積（０．００８）
Ａｔａｎｋ：サージタンク断面積（０．１００）
Ｑｓｕｒｇｅ：流量（０．０１２）
ｈ′（Ｑｓｕｒｇｅ）＝２．６７×１０９Ｑｓｕｒｇｅ３−１．２×１０８Ｑｓｕｒｇｅ２＋１．２×１０５Ｑｓｕｒｇｅ＋４７２．４８７
：（揚程カーブを流量で微分したもの）
Ｖｓｕｒｇｅ：圧力損失係数（０．１４４）
ｘｓｕｒｇｅ：サージタンクの水位変動（初期値＝０．０２）
ｖｓｕｒｇｅ：サージタンクの水位変動速度（初期＝０．００）
αｓｕｒｇｅ：サージタンクの水位変動加速度
であり、これらはＳＩ単位系で示している。【００９４】  この場合、解析方程式は、式（５７）から、次式（５８）で示される。すなわち、
2006103996.tif 000041
によって示される。したがって、式（５８）に、各種パラメータおよび、サージタンクの初期水位変動ｘｓｕｒｇｅ（ｔ０）＝０．０２、サージタンクの初期水位変動速度ｖｓｕｒｇｅ（ｔ０）＝０．００を代入して、サージタンクの水位変動加速度αｓｕｒｇｅを求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の水位変動速度ｖｓｕｒｇｅ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の水位変動ｘｓｕｒｇｅ（ｔ＋Δｔ）を求める。これら求められたｖｓｕｒｇｅ（ｔ＋Δｔ），ｘｓｕｒｇｅ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の水位変動加速度αｓｕｒｇｅ（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。【００９５】  図４４は、図４３に示したサージタンクの水位変動の動解析結果を示す図である。この動解析結果は、図４４に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００９６】（応用例１２）
図４５は、動解析対象である、ＬＣＲ回路の概要構成を示す模式図である。このＬＣＲ回路の電荷に固有の運動方程式は、次式（５９）で示される。
2006103996.tif 000042
ただし、ｇ：重力加速度（９．８０７）
Ｌｌｃｒ：インダクタンス（０．２０）
Ｒｌｃｒ：抵抗（５．００）
ｃｌｃｒ：キャパシタンス（０．００２）
Ｅｌｃｒ：電圧（初期値＝５．００）
ｑ  ：電荷（初期値＝０．００）
ｉ  ：電流（初期値＝０．００）
であり、これらはＳＩ単位系で示している。【００９７】  この場合、解析方程式は、式（５９）から、次式（６０）で示される。すなわち、
2006103996.tif 000043
によって示される。したがって、式（６０）に、各種パラメータおよび、初期電荷ｑ（ｔ０）＝０．００、初期電流ｉ（ｔ０）＝０．００を代入して、電荷の二階時間微分量（ｄ２ｑ／ｄｔ２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の電流ｉ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の電荷ｑ（ｔ＋Δｔ）を求める。これら求められたｉ（ｔ＋Δｔ），ｑ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．００１（ｓ））経過後の電荷の二階時間微分量（ｄ２ｑ（ｔ＋Δｔ）／ｄｔ２）を求め、上述した処理を所定回数繰り返す。【００９８】  図４６は、図４５に示したＬＣＲ回路の電荷変動の動解析結果を示す図である。この動解析結果は、図４６に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。【００９９】ところで、上述した数値解析装置の一例である動解析装置は、あらかじめ用意されたプログラムをパーソナル・コンピュータやワークステーションなどのコンピュータシステムで実行することによって実現することができる。そこで、上述した動解析装置と同様の機能を有する動解析プログラムを実行するコンピュータシステムについて説明する。【０１００】  図４７は、上述した動解析装置を実現するコンピュータシステムの構成を示すシステム構成図であり、図４８は、図４７に示したコンピュータシステムにおける本体部の構成を示すブロック図である。【０１０１】  図４７に示すように、このコンピュータシステム２００は、本体部２０１と、本体部２０１からの指示によって表示画面２０２ａに画像などの情報を表示するためのディスプレイ２０２と、このコンピュータシステム２００に種々の情報を入力するためのキーボード２０３と、ディスプレイ２０２の表示画面２０２ａ上の任意の位置を指定するためのマウス２０４とを備える。【０１０２】  また、このコンピュータシステム２００における本体部２０１は、図４８に示すように、ＣＰＵ２２１、ＲＡＭ２２２、ＲＯＭ２２３、ハードディスクドライブ（ＨＤＤ）２２４、ＣＤ−ＲＯＭ２０９を受け入れるＣＤ−ＲＯＭドライブ２２５、フレキシブルディスク（ＦＤ）２０８を受け入れるＦＤドライブ２２６、ディスプレイ２０２、キーボード２０３ならびにマウス２０４を接続するＩ／Ｏインターフェース２２７、ローカルエリアネットワークまたはワイドエリアネットワーク（ＬＡＮ／ＷＡＮ）２０６に接続するＬＡＮインターフェース２２８を備える。【０１０３】  さらに、このコンピュータシステム２００には、インターネットなどの公衆回線２０７に接続するためのモデム２０５が接続されるとともに、ＬＡＮインターフェース２２８およびＬＡＮ／ＷＡＮ２０６を介して、他のコンピュータシステム（ＰＣ）２１１、サーバ２１２並びにプリンタ２１３などが接続される。【０１０４】  そして、このコンピュータシステム２００は、所定の記録媒体に記録された動解析プログラムを読み出して実行することで実施の形態１，２に示した動解析装置を実現する。【０１��５】  ここで、所定の記録媒体とは、フレキシブルディスク（ＦＤ）２０８、ＣＤ−ＲＯＭ２０９、ＭＯディスク、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」の他に、コンピュータシステム２００の内外に備えられるハードディスクドライブ（ＨＤＤ）２２４や、ＲＡＭ２２２、ＲＯＭ２２３などの「固定用の物理媒体」、さらに、モデム２０５を介して接続される公衆回線２０７や、他のコンピュータシステム２１１ならびにサーバ２１２が接続されるＬＡＮ／ＷＡＮ２０６などのように、プログラムの送信に際して短期にプログラムを保持する「通信媒体」など、コンピュータシステム２００によって読み取り可能な動解析プログラムを記録する、あらゆる記録媒体を含むものである。【０１０６】  すなわち、動解析プログラムは、上記した「可搬用の物理媒体」、「固定用の物理媒体」、「通信媒体」などの記録媒体に、コンピュータ読み取り可能に記録されるものであり、コンピュータシステム２００は、このような記録媒体から動解析プログラムを読み出して実行することで、実施の形態１，２に示した動解析装置を実現する。【０１０７】  なお、動解析プログラムは、コンピュータシステム２００によって実行されることに限定されるものではなく、他のコンピュータシステム２１１またはサーバ２１２が動解析プログラムを実行する場合や、これらが協働して動解析プログラムを実行するような場合にも、本発明を同様に適用することができる。【０１０８】  なお、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示した通りに構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。【０１０９】  また、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵおよび当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
【産業上の利用可能性】
【０１１０】  以上のように、この発明は、自然現象や機械運動などにおける種々の運動態様を解析することができる種々の制御分野に用いられる数値解析装置および数値解析プログラムに有用である。
【図１】
2006103996.tif 000044
【図２】
2006103996.tif 000045
【図３】
2006103996.tif 000046
【図４】
2006103996.tif 000047
【図５】
2006103996.tif 000048
【図６】
2006103996.tif 000049
【図７】
2006103996.tif 000050
【図８】
2006103996.tif 000051
【図９】
2006103996.tif 000052
【図１０】
2006103996.tif 000053
【図１１】
2006103996.tif 000054
【図１２】
2006103996.tif 000055
【図１３】
2006103996.tif 000056
【図１４】
2006103996.tif 000057
【図１５】
2006103996.tif 000058
【図１６】
2006103996.tif 000059
【図１７】
2006103996.tif 000060
【図１８】
2006103996.tif 000061
【図１９】
2006103996.tif 000062
【図２０】
2006103996.tif 000063
【図２１】
2006103996.tif 000064
【図２２】
2006103996.tif 000065
【図２３】
2006103996.tif 000066
【図２４】
2006103996.tif 000067
【図２５】
2006103996.tif 000068
【図２６】
2006103996.tif 000069
【図２７】
2006103996.tif 000070
【図２８】
2006103996.tif 000071
【図２９】
2006103996.tif 000072
【図３０】
2006103996.tif 000073
【図３１】
2006103996.tif 000074
【図３２】
2006103996.tif 000075
【図３３】
2006103996.tif 000076
【図３４】
2006103996.tif 000077
【図３５】
2006103996.tif 000078
【図３６】
2006103996.tif 000079
【図３７】
2006103996.tif 000080
【図３８】
2006103996.tif 000081
【図３９】
2006103996.tif 000082
【図４０】
2006103996.tif 000083
【図４１】
2006103996.tif 000084
【図４２】
2006103996.tif 000085
【図４３】
2006103996.tif 000086
【図４４】
2006103996.tif 000087
【図４５】
2006103996.tif 000088
【図４６】
2006103996.tif 000089
【図４７】
2006103996.tif 000090
【図４８】
2006103996.tif 000091
【国際調査報告】
2006103996.tif 000092
2006103996.tif 000093
2006103996.tif 000094
2006103996.tif 000095
2006103996.tif 000096
2006103996.tif 000097
2006103996.tif 000098