JPWO2006103996A1 - 数値解析装置および数値解析プログラム - Google Patents

Abstract

解析対象の２階微分方程式を動解析する動解析装置１であって、前記２階微分方程式の初期条件を含むパラメータ設定および前記動解析の解析微小時間Δｔを設定する算入モジュール２１と、解析開始時刻である基準時刻における前記２階微分方程式の２階時間微分量α（ｔ）を算出する二階時間微分量算出モジュール２２と、前記解析微小時間経過後の１階時間微分量ｖ（ｔ＋Δｔ）＝ｖ（ｔ）＋α（ｔ）・Δｔを算出する一階時間微分量算出モジュール２３と、前記解析微小時間経過後の物理量ｘ（ｔ＋Δｔ）＝ｘ（ｔ）＋ｖ（ｔ＋Δｔ）・Δｔを算出する物理量算出モジュール２４と、時刻（ｔ＋Δｔ）を新たな基準時刻として設定し、各モジュールに対してそれぞれ設定された基準時刻を用いて前記２階時間微分量、前記１階時間微分量、および前記物理量を算出させる制御を繰り返し行う演算制御モジュール２０とを備える。

Description

この発明は、自然現象や機械運動などにおける種々の運動態様を解析することができる数値解析装置および数値解析プログラムに関するものである。

従来から、常微分方程式を用いた数値解析法には、区間の初期変位、初期速度、初期加速度などを用いて、つぎの時刻の変位や速度を求めるオイラー法が広く用いられている（特許文献１，２参照）。

特開２００４−０６２６７６号公報 特開２００２−２５８９３３号公報

しかしながら、上述したオイラー法は、微小時間Δｔ毎に処理ステップを重ねる毎にその解析結果に比較的大きな誤差が生じるという問題点があった。さらに、この誤差を抑えるために微小時間Δｔを小さくすると、解析結果を得るために膨大な処理時間を要し、実用に供しないという問題点があった。

この発明は、上記に鑑みてなされたものであって、解析結果に大きな誤差を生じさせず、簡易かつ短時間に数値解析処理を行うことができる数値解析装置および数値解析プログラムを提供することを目的とする。

上述した課題を解決し、目的を達成するために、この発明にかかる数値解析装置は、解析対象の２階微分方程式を数値解析する数値解析装置であって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、を備えたことを特徴とする。

また、この発明にかかる数値解析装置は、解析対象の２階微分方程式を数値解析する数値解析装置であって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、を備えたことを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記設定手段は、前記数値解析の終了条件を設定し、前記演算制御手段は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記物理量算出手段によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手段を備えたことを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記出力処理手段は、前記解析微小変数量の増分に対応する、物理量算出手段によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記２階微分方程式は、変数量項を含むことを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記２階微分方程式は、時間項を含むことを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、前記２階微分方程式は、多元連立方程式であることを特徴とする。

また、この発明にかかる数値解析装置は、上記の発明において、表計算プログラムを用いて数値解析を行うことを特徴とする。

また、この発明にかかる数値解析プログラムは、解析対象の２階微分方程式を数値解析する数値解析プログラムであって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、をコンピュータに実行させることを特徴とする。

また、この発明にかかる数値解析プログラムは、解析対象の２階微分方程式を数値解析する数値解析プログラムであって、前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、をコンピュータに実行させることを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記設定手順は、前記数値解析の終了条件を設定し、前記演算制御手順は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記物理量算出手順によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手順を備えたことを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記出力処理手順は、前記解析微小変数量の増分に対応する、物理量算出手順によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記２階微分方程式は、変数量項を含むことを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記２階微分方程式は、時間項を含むことを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、前記２階微分方程式は、多元連立方程式であることを特徴とする。

また、この発明にかかる数値解析プログラムは、上記の発明において、表計算プログラムを用いて数値解析を行うことを特徴とする。

この発明にかかる数値解析装置および数値解析プログラムによれば、設定手段が、２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定し、２階微分量算出手段が、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出し、１階微分量算出手段が、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出し、物理量算出手段が、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出し、演算制御手段が、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行うようにしているので、解析結果に大きな誤差を生じさせず、簡易かつ短時間に数値解析処理を行うことができるという効果を奏する。

また、この発明にかかる数値解析装置および数値解析プログラムによれば、設定手段が、２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定し、２階微分量算出手段が、解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出し、１階微分量算出手段が、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出し、物理量算出手段が、前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出し、演算制御手段が、前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行うようにしているので、解析結果に大きな誤差を生じさせず、簡易かつ短時間に数値解析処理を行うことができるという効果を奏する。

図１は、この発明の実施の形態１にかかる数値解析装置である動解析装置の概要構成を示すブロック図である。 図２は、図１に示した動解析装置による動解析処理手順を示すフローチャートである。 図３は、この発明の実施の形態１による動解析処理の概念を示す図である。 図４は、従来のオイラー法による動解析処理の概念を示す図である。 図５は、従来の修正オイラー法による動解析処理の概念を示す図である。 図６は、解析方程式の対象である、ばね−質点系の概念構成を示す模式図である。 図７は、図６に示したばね−質点系の動解析結果を示す図である。 図８は、図７に示したばね−質点系を動解析する表計算プログラムの表示画面の一例を示す図である。 図９は、図６に示したばね−質点系に対する動解析処理を表計算プログラムによって行った一例を示す図である。 図１０は、図９に示した結果をグラフ化した図である。 図１１は、同じばね−質点系に対する動解析処理をオイラー法によって行った一例を示す図である。 図１２は、図１１に示した動解析結果をグラフ化した図である。 図１３は、同じばね−質点系に対する動解析を修正オイラー法によって行った結果をグラフ化した図である。 図１４は、この発明の実施の形態２である動解析装置の概要構成を示すブロック図である。 図１５は、図１４に示した動解析装置による動解析処理手順を示すフローチャートである。 図１６は、この発明の実施の形態２による動解析処理の概念を示す図である。 図１７は、動解析対象であるばね質点ダッシュポット系の概要構成を示す模式図である。 図１８は、図１７に示した質点ダッシュポット系の動解析結果を示す図である。 図１９は、動解析対象である、単振動と回転運動とが伴った運搬機の概要構成を示す模式図である。 図２０は、図１９に示した運搬機の変位の動解析結果を示す図である。 図２１は、図１９に示した運搬機の角変位の動解析結果を示す図である。 図２２は、動解析対象である、ねじり軸、継手、ギア等の動力伝達器で連結された複数の回転体の概要構成を示す模式図である。 図２３は、図２２に示した動力伝達器で結合された各回転体の角変位の動解析結果を示す図である。 図２４は、図２２に示した動力伝達器で結合された各回転体の角変位の動解析結果を示す図である。 図２５は、動解析対象である、調速機のすべり子およびタービンの概要構成を示す模式図である。 図２６は、図２５に示した調速機のすべり子の変位の動解析結果を示す図である。 図２７は、図２５に示したタービンの角変位の動解析結果を示す図である。 図２８は、動解析対象である、旋削機などのバイトの概要構成を示す模式図である。 図２９は、図２８に示した旋削機などのバイト変位の動解析結果を示す図である。 図３０は、動解析対象である、クランクを有する往復機械の概要構成を示す模式図である。 図３１は、図３０に示したクランクを有する往復機械の動解析結果を示す図である。 図３２は、図３０に示した往復機械のピストンの動解析結果を示す図である。 図３３は、動解析対象である、ブレーキディスクの概要構成を示す模式図である。 図３４は、図３３に示したブレーキディスクの動解析結果を示す図である。 図３５は、動解析対象である、剛体リンクで接続されたマニピュレータ状機器の概要構成を示す模式図である。 図３６は、図３５に示した剛体リンクで接続されたマニピュレータ状機器の動解析結果を示す図である。 図３７は、図３５に示した剛体リンクで接続されたマニピュレータ状機器の動解析結果を示す図である。 図３８は、動解析対象である、重力等を相互に及ぼし合う粒子あるいは天体の概要構成を示す模式図である。 図３９は、図３８に示した地球のｘ方向変位の動解析結果を示す図である。 図４０は、太陽、地球、月の三体問題を解析する座標系および配置関係を示す図である。 図４１は、３体問題における地球を中心とした月の軌道の動解析結果を示す図である。 図４２は、従来のオイラー法を用いて地球を中心とした月の軌道の動解析結果を示す図である。 図４３は、動解析対象である、遠心ポンプを有するシステムのサージタンクの概要構成を示す模式図である。 図４４は、図４３に示したサージタンクの水位変動の動解析結果を示す図である。 図４５は、動解析対象である、ＬＣＲ回路の概要構成を示す模式図である。 図４６は、図４５に示したＬＣＲ回路の電荷変動の動解析結果を示す図である。 図４７は、上述した動解析装置を実現するコンピュータシステムの構成を示すシステム構成図である。 図４８は、図４７に示したコンピュータシステムにおける本体部の構成を示すブロック図である。

符号の説明

１ 動解析装置
１０ 入力部
１１ 出力部
１２ 演算部
１３ 出力処理部
１４ 記憶部
２０ 演算制御モジュール
２１ 算入モジュール
２２ 二階時間微分量算出モジュール
２３ 一階時間微分量算出モジュール
２４，２４ａ 物理量算出モジュール

以下、この発明を実施するための最良の形態である数値解析装置の一例としての動解析装置および数値解析プログラムの一例としての動解析プログラムについて説明する。

（実施の形態１）
図１は、この発明の実施の形態１にかかる数値解析装置の一例である動解析装置の概要構成を示すブロック図である。図１に示すように、この動解析装置１は、キーボードやポインティングデバイスによって実現され、各種情報を入力する入力部１０と、液晶ディスプレイやプリンタなどによって実現され、処理結果を出力する出力部１１と、各種プログラムやデータを格納するとともに、演算結果などを一時的に格納する記憶部１４と、動解析処理の実行制御を行う制御部Ｃとを有する。

制御部Ｃは、演算部１２と出力処理部１３とを有し、演算部１２が動解析処理を行い、出力処理部１３は、動解析処理された結果を出力部１１に出力させるための処理を行う。演算部１２は、演算制御モジュール２０、算入モジュール２１、二階時間微分量算出モジュール２２、一階時間微分量算出モジュール２３、および物理量算出モジュール２４を有する。

ここで、図２に示したフローチャートを参照して制御部Ｃの制御処理手順について説明する。図１および図２において、まず入力部１０および出力部１１を用い、算入モジュール２２による解析情報の設定入力を行う（ステップＳ１０１）。解析情報とは、解析対象となる解析方程式α（ｔ）＝ｆ（ｘ（ｔ），ｖ（ｔ））、解析微小時間Δｔ、さらには、パラメータ、初期条件、終了条件などが含まれる。ここでいう解析方程式は、二階常微分方程式であり、α（ｔ）は、加速度などの二階時間微分量を示し、ｘ（ｔ）は、位置などの物理量を示し、ｖ（ｔ）は、速度などの一階時間微分量を示している。

二階時間微分量算出モジュール２２は、初期条件をもとに基準時刻ｔにおける加速度α（ｔ）を算出する（ステップＳ１０２）。その後、一階時間微分量算出モジュール２３は、算出された加速度α（ｔ）に解析微小時間Δｔを乗算し、この乗算値に基準時刻ｔの速度ｖ（ｔ）を加算し、この加算値を解析微小時間Δｔ経過後の速度ｖ（ｔ＋Δｔ）として算出する（ステップＳ１０３）。すなわち、ｖ（ｔ＋Δｔ）＝ｖ（ｔ）＋α（ｔ）・Δｔを算出する。

その後、物理量算出モジュール２４は、算出された速度ｖ（ｔ＋Δｔ）に解析微小時間Δｔを乗算し、この乗算値に基準時刻ｔの変位ｘ（ｔ）を加算し、この加算値を解析微小時間Δｔ経過後の変位ｘ（ｔ＋Δｔ）として算出する。すなわち、ｘ（ｔ＋Δｔ）＝ｘ（ｔ）＋ｖ（ｔ＋Δｔ）・Δｔを算出する。

その後、演算部１２は、ステップＳ１０２〜Ｓ１０４で算出した値などを記憶部１４に記憶させるとともに、出力処理部１３を介して出力部１１から出力できる出力処理を行う（ステップＳ１０５）。さらに、演算部１２は、ｔが所定値を超えたか否かを判断する（ステップＳ１０６）。ここで、所定値とは基準時刻ｔからの経過時間である。この経過時間は、この後繰り返す処理回数をｎとすると、Δｔ×ｎとして示すことができる。ｔが所定値を超えていない場合（ステップＳ１０６，Ｎｏ）には、基準時刻ｔに解析微小時間Δｔを加えた値（ｔ＋Δｔ）を新たな基準時刻として設定してステップＳ１０２に移行し（ステップＳ１０７）、上述した処理を繰り返す。一方、ｔが所定値を超えた場合（ステップＳ１０６，Ｙｅｓ）には、本処理を終了する。

すなわち、この動解析装置１による動解析処理は、解析微小時間Δｔ後における変位ｘ（ｔ＋Δｔ）を次のようにして求める。まず、基準時刻ｔにおける解析方程式である加速度α（ｔ）を、基準時刻ｔにおける物理量である変位ｘ（ｔ）および基準時刻ｔにおける一階時間微分量である速度ｖ（ｔ）をもとに、式（１）によって算出する。その後、この解析微小時間Δｔ経過後における速度ｖ（ｔ＋Δｔ）を、式（２）をもとに求め、さらに解析微小時間Δｔ経過後における変位ｘ（ｔ＋Δｔ）を、式（３）をもとに求める。

ところで、従来の動解析装置では、オイラー法を用いるものがあり、このオイラー法では、上述した式（３）におけるｖ（ｔ＋Δｔ）をｖ（ｔ）とする次式（４）によって求めていた。すなわち、

によって変位（ｔ＋Δｔ）を求めていた。しかし、このオイラー法では、単振動などの動解析を行うと発散してしまう場合があった。そこで、式（３）を次式（５）によって求める修正オイラー法が考え出された。すなわち、

である。この修正オイラー法では、ｖ（ｔ＋Δｔ）を、基準時刻ｔと、基準時刻ｔから解析微小時間Δｔ経過後の時刻（ｔ＋Δｔ）との間の平均速度（１／２・（ｖ（ｔ）＋ｖ（ｔ＋Δｔ）））に置き換えて動解析を行うようにしている。この修正オイラー法によれば、オイラー法に比較して発散の程度を小さくすることができるが、依然として発散した動解析となってしまい、精度が低い動解析結果となってしまっていた。したがって、オイラー法や修正オイラー法を用いた動解析では、解析微小時間Δｔを極端に小さくし、多大な時間をかけて精度を確保せざるを得なかった。これに対し、上述した式（１）〜（３）を繰り返し行う動解析処理では、解析微小時間Δｔが短くても精度の高い動解析結果を得ることができる。

この発明による動解析処理、従来のオイラー法、修正オイラー法の概念を図化すると図３〜図５に示すようになる。図３は、この発明の実施の形態１による動解析処理の概念を示す図であり、図４は、従来のオイラー法による動解析処理の概念を示す図であり、図５は、従来の修正オイラー法による動解析処理の概念を示す図である。図３〜図５において、この発明による動解析処理では、図３に示すように、時刻ｔにおける変位および速度から加速度を求め、この時刻ｔにおける速度と加速度とを用いて解析微小時間Δｔ経過時刻（ｔ＋Δｔ）の速度を求める。その後、この時刻（ｔ＋Δｔ）の速度と時刻ｔの変位を用いて時刻（ｔ＋Δｔ）の変位を求めるようにしており、特に、時刻（ｔ＋Δｔ）時の速度を用いて時刻（ｔ＋Δｔ）の変位を求めることが特徴である。一方、オイラー法では、図４に示すように、時刻ｔにおける速度と時刻ｔにおける変位を用いて、時刻（ｔ＋Δｔ）の変位を求めるようにしている。さらに、従来の修正オイラー法では、図５に示すように、時刻ｔと時刻（ｔ＋Δｔ）との間の平均速度（１／２・（ｖ（ｔ）＋ｖ（ｔ＋Δｔ）））を用いて、時刻（ｔ＋Δｔ）における変位を求めるようにしている。

ここで、具体的な単振動の解析方程式に対する処理結果について説明する。図６は、解析方程式の対象である、ばね−質点系の概念構成を示す模式図である。図６に示したばね−質点系では、単振動の運動状態を呈する。図６に示すばね−質点系の運動方程式は、次式（６）に示すように、

となる。なお、このばね−質点系では、ばね定数ｋ＝４π^２（Ｎ／ｍ）、質量ｍ＝１．０（ｋｇ）、初期変位ｘ（０）＝２．０（ｍ）、初期速度ｖ（０）＝０．０（ｍ）である。この運動方程式の解析方程式は、次式（７）に示すように、

となり、パラメータを代入すると、次式（８）に示す解析方程式となる。すなわち、

となる。なお、このばね−質点系の代数解は、次式（９）に示すように、

として表される。解析微小時間Δｔ＝０．０１（ｓ）としている。

図７は、図６に示したばね−質点系の解析結果を示す図である。この動解析は、解析微小時間Δｔ＝０．０１（ｓ）で行っている。図７に示すように、オイラー法および修正オイラー法では、時間の経過とともに発散しているのに対し、この発明による動解析処理では、発散もなくほぼ代数解に一致した結果を得ることができている。この結果、この発明の実施の形態１によれば、オイラー法や修正オイラー法に比して短時間に精度の高い動解析を行うことができる。換言すれば、この実施の形態１による動解析処理は、同じ解析時間では、オイラー法や修正オイラー法に比して格段に精度の高い動解析を行うことができる。

この実施の形態１による動解析処理は、複雑なプログラムではないため、たとえば表計算プログラムを用いたソフトウェアによって容易に実現することができる。図８は、図７に示したばね−質点系を動解析する表計算プログラムの表示画面の一例を示す図である。図８において、まず時間列のセルＣ１１には初期時刻ｔ＝０．００を設定し、その後のセルＣ２１〜Ｃｎ１には、微小解析時間Δｔ＝０．０１が加算されるようにしている。基準時刻ｔ時点を示す区間初期の変位列のセルＣ１２には、初期変位ｘ（０）＝２．００が設定され、区間初期の速度列のセルＣ１３には、初期速度ｖ（０）＝０．００が設定される。加速度列のセルＣ１４には、セルＣ１２，Ｃ１３に示された変位および速度をもとに加速度が演算され、その結果が表示される。基準時刻ｔから微小解析時間Δｔ経過後を示す区間終期の速度列のセルＣ１５には、式（２）に対応した演算結果が表示され、区間終期の速度列のセルＣ１６には、式（３）に対応した演算結果が表示される。そして、このセルＣ１５，Ｃ１６に示された結果は、次の微小解析時間Δｔ経過後の時間行のセルＣ２２，Ｃ２３には、次の区間初期の変位および速度としてコピーされる。さらに、セルＣ２２，Ｃ２３に示された変位および速度をもとに、加速度が求められ、結果がセルＣ２４に表示され、これらセルＣ２２，Ｃ２３，Ｃ２４に示された変位、速度、加速度をもとに次の微小解析時間Δｔ経過後の速度および変位が求められ、セルＣ２５，Ｃ２６に表示される。以後、同様にして、区間終期における速度および変位が次の微小解析時間Δｔ経過後の区間初期の速度および変位として用いられ、この時点における加速度が求められて表示されるとともに、次の微小解析時間Δｔ経過後の区間終期における速度および変位が求められて表示される。この処理は、時間列で設定した繰り返し数すなわち行数分、繰り返される。また、これらの結果は、表計算ソフトウェアに通常付加されたグラフ化処理によってグラフ化することができる。このグラフ化処理は、出力処理部１３が行うことになる。

図９は、図６に示したばね−質点系に対する動解析処理を表計算プログラムによって行った一例を示す図である。なお、ここでは、初期変位ｘ（０）＝１．０００（ｍ）としている。また、図１０は、図９に示した結果をグラフ化したものであり、横軸に時間をとり、縦軸に変位をとっている。図１０に示すように、代数解と数値解析（動解析）結果はほぼ一致しており、その誤差は、±０．０４２ｍであった。また、図１０では、１０（ｓ）までの動解析であったが、さらに１００（ｓ）までの動解析を行っても、誤差は、±０．１３４ｍであった。

なお、図１１は、同じばね−質点系に対する動解析処理をオイラー法によって行った一例を示す図であり、図１２は、図１１に示した動解析結果をグラフ化したものである。さらに、図１３は、同じばね−質点系に対する動解析を修正オイラー法によって行った結果をグラフ化したものである。図１１および図１２に示すように、同じ微小解析時間によるオイラー法によって動解析した場合、時間の経過とともに動解析結果と代数解との誤差が大きくなり、発散している。この場合における誤差は、１０（ｓ）で±６ｍ程度まで広がっている。さらに、修正オイラー法によっても、図１３に示すように、時間の経過とともに動解析結果と代数解との誤差が大きくなり、発散している。この場合における誤差は、オイラー法よりも小さいが、１０（ｓ）で±２ｍ程度となっている。

（実施の形態２）
つぎに、実施の形態２について説明する。上述した実施の形態１では、現時点の速度などの一階時間微分量および現時点の加速度などの二階微分量をもとに次の時点の一階時間微分量を求め、さらにこの次の時点の一階時間微分量および現時点の変位などの物理量を用いて次の時点の物理量を求めるようにしていたが、この実施の形態２では、現時点の一階時間微分量、現時点の二階時間微分量および現時点の物理量をもとに、一度で次の時点の物理量を求めるようにしている。ただし、次の時点の二階時間微分量を算出する際に用いる次の時点の一階時間微分量は、現時点において求めておく。

図１４は、この発明の実施の形態２である数値解析装置の一例である動解析装置の概要構成を示すブロック図である。図１４に示すように、この動解析装置は、物理量算出モジュール２４に替えて、物理量算出モジュール２４ａを設けている。その他の構成は、実施の形態１と同じであり、同一構成部分には同一符号を付している。

二階時間微分量算出モジュール２２は、一階時間微分量算出モジュール２３が求めておいた速度ｖ（ｔ）と、物理量算出モジュール２４が求めておいた基準時刻ｔにおける変位ｘ（ｔ）とを用いて、加速度α（ｔ）を求め、物理量算出モジュール２４ａは、この加速度α（ｔ）と、速度ｖ（ｔ）と、基準時刻ｔにおける変位ｘ（ｔ）とを用いて、直接、解析微小時間Δｔ経過後の変位ｘ（ｔ＋Δｔ）を求めている。これは、式（３）に式（２）を代入した次式（１０）によって直接、変位ｘ（ｔ＋Δｔ）を求めることができることを意味する。すなわち、

によって求めることができる。

図１５は、図１４に示した動解析装置による動解析処理手順を示すフローチャートである。図１５に示した動解析処理手順では、実施の形態１におけるステップ１０４の処理をステップＳ２０４の処理に置き換えた処理を行っている。その他のステップＳ２０１〜Ｓ２０３，Ｓ２０５〜Ｓ２０７は、実施の形態１に示したステップＳ１０１〜Ｓ１０３，Ｓ１０５〜Ｓ１０７と同じ処理である。

すなわち、ステップＳ２０４では、ステップＳ２０２によって算出された加速度α（ｔ）に解析微小時間Δｔの２乗値を乗算し、さらに基準時刻ｔの速度ｖ（ｔ）に解析微小時間Δｔを乗算し、これら乗算した各値と基準時刻ｔの変位ｘ（ｔ）とを加算し、この加算値を解析微小時間Δｔ経過後の変位ｘ（ｔ＋Δｔ）として算出するようにしている。なお、ステップＳ２０３で基準時刻ｖ（ｔ）の算出処理を行うのは、実施の形態１と異なり、次の時点の変位を算出する場合に速度ｖ（ｔ）が求められていないためである。

これによって、この実施の形態２では、実施の形態１とは異なり、現時点の速度、加速度および変位をもとに、一度に解析微小時間Δｔ経過後の速度ｖ（ｔ＋Δｔ）を算出することができる。

なお、この発明の実施の形態２による動解析処理の概念を図化すると図１６に示すようになる。図１６は、この発明の実施の形態２による動解析処理の概念を示す図である。図１６において、この発明の実施の形態２による動解析処理では、時刻ｔにおける変位および速度から加速度を求め、この時刻ｔにおける速度と加速度とを用いて解析微小時間Δｔ経過時刻（ｔ＋Δｔ）の変位を求めるようにしている。

ここで、上述した実施の形態１，２では、解析方程式α（ｔ）＝ｆ（ｘ（ｔ），ｖ（ｔ））として説明したが、速度項ｖ（ｔ）を含まない解析方程式α（ｔ）＝ｆ（ｘ（ｔ））であってもよい。また、解析方程式α（ｔ）には、時間項を含んでも良い。すなわち、α（ｔ）＝ｆ（ｘ（ｔ），ｖ（ｔ），ｔ）あるいはα（ｔ）＝ｆ（ｘ（ｔ），ｔ）であってもよい。さらに、解析対象によっては、複数の解析方程式が連立方程式であってもよい。すなわち、解析方程式は、多元方程式であってもよい。この多元方程式は、一般化すると、αｉ（ｔ）＝ｆ（ｘ１（ｔ），…，ｘｎ（ｔ），ｖ１（ｔ），…，ｖｎ（ｔ），ｔ）とするｎ元の連立微分方程式として表すことができる。ただし、ｉ＝１〜ｎである。

なお、上述した実施の形態１，２では、数値解析装置の一例として動解析装置をあげて説明したが、時間を変数量とする動解析に限らず、一般の数値解析にも適用することができる。すなわち、時間ｔを変数量とするのではなく、温度Ｔや距離Ｌを変数量とした数値解析を行うことができる。たとえば、時刻ｔおよび解析微小時間Δｔに替えて、温度Ｔおよび解析微小温度ΔＴとすれば、温度変化に伴う数値解析を行うことができる。

つぎに、上述した実施の形態１，２の応用例として、ばね質点ダッシュポット系の動解析処理、運搬機器の動解析処理、動力伝達機等で連結された複数の回転体の動解析処理、調速機の動解析処理、旋削機等のバイトの動解析処理、クランクを有する往復機械の動解析処理、摩擦を利用したブレーキの動解析処理、剛体リンクで接続されたマニピュレータ状機器等の動解析処理、重力等を相互に及ぼし合う粒子あるいは天体の動解析処理、太陽・地球・月の三体問題の動解析処理、遠心ポンプを有するシステムのサージング特性解析処理、およびインダクタンス素子、コンデンサ素子、抵抗素子を含む電気・電子回路に生じる電流あるいは電荷の動解析処理について説明する。

（応用例１）
図１７は、動解析対象であるばね質点ダッシュポット系の概要構成を示す模式図である。このばね質点ダッシュポット系に固有の運動方程式は、次式（１１）で示される。

ただし、ｍ：質点の質量（１５００．００）
ｃ：ダッシュポットの粘性減衰係数（４２４２．６４）
ｋ：バネ係数（１２０００）
ｆ：外力の振幅（１０００．０）
Ｔ：外力の周期（１．００）
ｔ：基準時刻
ｘ：質点の変位
ｖ：質点の速度
α：質点の加速度
ｔ_０：解析開始時刻
ｘ_０：質点の初期変位（０．１０）
ｖ_０：質点の初期速度（１．００）
α_０：質点の初期加速度
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（１１）から、次式（１２）で示される。すなわち、

によって示される。したがって、式（１２）で示したα（ｔ）に、各種パラメータおよび、ｔ＝ｔ_０、ｘ（ｔ_０）＝ｘ_０、ｖ（ｔ_０）＝ｖ_０を代入してα（ｔ_０）を求め、式（２）および式（３）で示す演算を行って解析微小時間Δｔ（＝０．０１（ｓ））経過後のｖ（ｔ＋Δｔ）およびｘ（ｔ＋Δｔ）を求める。このｖ（ｔ＋Δｔ）およびｘ（ｔ＋Δｔ）を用いて次の解析微小時間Δｔ経過後の加速度α（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。

図１８は、図１７に示した質点ダッシュポット系の動解析結果を示す図であり、横軸に時間をとり、縦軸に変位をとっている。この動解析結果は、図１８に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例２）
図１９は、動解析対象である、単振動と回転運動とが伴った運搬機の概要構成を示す模式図である。この運搬機に固有の運動方程式は、次式（１３），（１４）で示される。

ただし、ｍ：運搬機の質量（１５００．００）
Ｊ：運搬機の慣性モーメント（３１２５）
ｋ_１：バネ定数１（３０００）
ｋ_２：バネ係数２（３０００）
Ｌ_１：重心からの距離１（２．００）
Ｌ_２：重心からの距離２（２．００）
ｘ：運搬機の変位（初期値ｘ_０＝０．１０）
ｖ：運搬機の速度（初期値ｖ_０＝０．００）
α：運搬機の加速度
θ：運搬機の角変位（初期値θ_０＝０．２６）
ω：運搬機の角速度（初期値ω_０＝０．００）
ｄ^２θ／ｄｔ^２：運搬機の角加速度
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（１３），（１４）から、次式（１５），（１６）で示される。すなわち、

によって示される。したがって、各式（１５），（１６）に、各種パラメータおよび、初期変位ｖ（ｔ_０）＝０．００、初期角変位θ_０（ｔ_０）＝０．２５を代入して、それぞれα（ｔ_０）および（ｄ^２θ／ｄｔ_０ ^２）を求め、式（１０）で示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後のｘ（ｔ＋Δｔ）およびθ（ｔ＋Δｔ）を求める。このｘ（ｔ＋Δｔ）およびθ（ｔ＋Δｔ）を用いて次の解析微小時間Δｔ経過後の加速度α（ｔ＋Δｔ）および（ｄ^２θ（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。なお、実施の形態１のように、式（２）および式（３）を用いる場合、それぞれ運搬機の速度ｖ（ｔ＋Δｔ）および運搬機の角速度ω（ｔ＋Δｔ）を式（３）によって求める。この場合、運搬機の初期速度ｖ（ｔ_０）＝０．００および運搬機の角速度ω（ｔ_０）＝０．００とする。

図２０は、図１９に示した運搬機の変位の動解析結果を示す図であり、図２１は、図１９に示した運搬機の角変位の動解析結果を示す図である。それぞれの動解析結果は、図２０および図２１に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例３）
図２２は、動解析対象である、ねじり軸、継手、ギア等の動力伝達器で連結された複数の回転体の概要構成を示す模式図である。この解析対象の運動方程式は、次式（１７），（１８）で示される。

ただし、Ｉ_１：回転体１の慣性モーメント（９５００）
Ｉ_２：回転体２の慣性モーメント（５３９００）
ｋ_１２：捩りバネ定数（９８９８００００）
θ_１：回転体１の角変位（初期値＝０．０００）
θ_２：回転体２の角変位（初期値＝０．０２６）
ω_１：回転体１の角速度（初期値＝０．０００）
ω_２：回転体２の角速度（初期値＝０．０００）
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（１７），（１８）から、次式（１９），（２０）で示される。すなわち、

によって示される。したがって、各式（１９），（２０）に、各種パラメータおよび、初期角変位θ_１０＝０．０００、初期角変位θ_２０＝０．０２６を代入して、それぞれ（ｄ^２θ_１／ｄｔ^２）および（ｄ^２θ_２／ｄｔ^２）を求め、式（１０）で示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後のθ_１（ｔ＋Δｔ）およびθ_２（ｔ＋Δｔ）を求める。このθ_１（ｔ＋Δｔ）およびθ_２（ｔ＋Δｔ）を用いて次の解析微小時間Δｔ経過後の（ｄ^２θ_１（ｔ＋Δｔ）／ｄｔ^２）および（ｄ^２θ_２（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。なお、実施の形態１のように、式（２）および式（３）を用いる場合、それぞれの角速度ω_１（ｔ＋Δｔ），ω_２（ｔ＋Δｔ）を式（３）によって求める。この場合、初期角速度ω_１（ｔ＋Δｔ），ω_２（ｔ＋Δｔ）＝０．０００とする。

図２３および図２４は、図２２に示した動力伝達器で結合された各回転体の角変位の動解析結果を示す図である。それぞれの動解析結果は、図２３および図２４に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例４）
図２５は、動解析対象である、調速機のすべり子およびタービンの概要構成を示す模式図である。この調速機に固有の運動方程式は、次式（２１），（２２）で示される。

ただし、ｍ_ｇｖ：調速機のすべり子の等価質量（１０．００）
ｃ_ｇｖ：調速機の粘性減衰係数（５０．６６）
ｋ_ｇｖ：調速機のバネ定数（１５７９．１４）
Ａ_ｇｖ：タービンの回転速度が１ｒａｄ／ｓ増すことによって調速機のすべり子に働く上向きの力（１０．００）
Ｊ_ｇｖ：タービンの等価慣性モーメント（２．５０）
Ｂ_ｇｖ：調速機のすべり子が単位長下がるために生じるタービンのトルク増加（１０００．００）
ｘ_ｇｖ：調速機のすべり子の変位（初期値＝０．００）
ｖ_ｇｖ：調速機のすべり子の速度（初期値＝０．００）
α_ｇｖ：調速機のすべり子の加速度
φ_ｇｖ：タービンの角変位（初期値＝０．００）
ω_ｇｖ：タービンの角速度（初期値＝０．２０）
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（２１），（２２）から、次式（２３），（２４）で示される。すなわち、

によって示される。したがって、各式（２３），（２４）に、各種パラメータおよび、初期変位ｘ_ｇｖ（ｔ_０）＝０．００、初期角変位φ_ｇｖ（ｔ_０）＝０．００、初期速度ｖ_ｇｖ（ｔ_０）＝０．００、初期角速度ω_ｇｖ（ｔ_０）＝０．２０を代入して、加速度α_ｇｖおよび角加速度（ｄ^２φ_ｇｖ／ｄｔ^２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０００５（ｓ））経過後の速度ｖ_ｇｖ（ｔ＋Δｔ）およびω_ｇｖ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の変位ｘ_ｇｖ（ｔ＋Δｔ）および角変位φ_ｇｖ（ｔ＋Δｔ）を求める。これら求められたｖ_ｇｖ（ｔ＋Δｔ）、ω_ｇｖ（ｔ＋Δｔ）、ｘ_ｇｖ（ｔ＋Δｔ）、およびφ_ｇｖ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０００５（ｓ））経過後の加速度α_ｇｖ（ｔ＋Δｔ）および角加速度（ｄ^２φ_ｇｖ（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。

図２６は、図２５に示した調速機のすべり子の変位の動解析結果を示す図であり、図２７は、図２５に示したタービンの角変位の動解析結果を示す図である。それぞれの動解析結果は、図２６および図２７に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例５）
図２８は、動解析対象である、旋削機などのバイトの概要構成を示す模式図である。この旋削機に固有の運動方程式は、次式（２５）で示される。

ただし、ｍ_ｂ：バイトの質量（１０．００）
ｃ_ｂ：粘性減衰係数（２０．００）
ｋ_ｂ：バネ定数（１０００．００）
Ｖ：切削速度（２．００）
Ｐ_ｂ：背分力（１００．００）
ｘ_ｂ：バイトの変位（初期値＝０．００）
ｖ_ｂ：バイトの速度（初期値＝１．００）
α_ｂ：バイトの加速度
μ（Ｖ−ｖ_ｂ）Ｐ_ｂ：０．００７５（Ｖ−ｖ_ｂ）^２−０．０５５（Ｖ−ｖ_ｂ）＋０．４００
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（２５）から、次式（２６）で示される。すなわち、

によって示される。したがって、式（２６）に、各種パラメータおよび、初期変位ｘ_ｂ（ｔ_０）＝０．００、初期速度ｖ_ｂ（ｔ_０）＝１．００を代入して、加速度α_ｂを求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０２（ｓ））経過後の速度ｖ_ｂ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０２（ｓ））経過後の変位ｘ_ｂ（ｔ＋Δｔ）を求める。これら求められたｖ_ｂ（ｔ＋Δｔ）、ｘ_ｂ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０２（ｓ））経過後の加速度α_ｂ（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。

図２９は、図２８に示した旋削機などのバイト変位の動解析結果を示す図である。この動解析結果は、図２９に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例６）
図３０は、動解析対象である、クランクを有する往復機械の概要構成を示す模式図である。このクランクを有する往復機械に固有の運動方程式は、次式（２７）で示される。

ただし、ｍ_ｐｉ：往復質量（１５．００）
Ｉ_ｃｒ：クランク軸周りの全回転質量の慣性モーメント（７．７３）
Ｒ_ｃｒ：クランク軸の腕の長さ（０．１３）
ε：Ｒ_ｃｒ／Ｌ_ｃｏｎ
Ｌ_ｃｏｎ：接続の長さ（０．２５）
Ｐ_ｐｉ：外力（０．００）
Ｔｒ_ｃｒ：抵抗トルク
θ_ｃｒ：クランクの角変位（初期値＝０．２６）
ｄθ_ｃｒ／ｄｔ：クランクの角速度（初期値＝１２．５７）
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（２７）から、次式（２８）で示される。すなわち、

によって示される。したがって、式（２８）に、各種パラメータおよび、初期角変位θ_ｃｒ（ｔ_０）＝０．２６、初期角速度（ｄθ_ｃｒ／ｄｔ_０）＝１２．５７を代入して、角加速度（ｄ^２θ_ｃｒ／ｄｔ_０ ^２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の角速度（ｄθ_ｃｒ（ｔ＋Δｔ）／ｄｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の角変位θ_ｃｒ（ｔ＋Δｔ）を求める。これら求められた（ｄθ_ｃｒ（ｔ＋Δｔ）／ｄｔ）、θ_ｃｒ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．００１（ｓ））経過後の角加速度（ｄ^２θ_ｃｒ（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。

図３１は、図３０に示したクランクを有する往復機械の動解析結果を示す図である。この動解析結果は、図３１に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

なお、ピストンの変位は、機構的に次式（２９）によって示すように、クランク角度などをもとに求めることができる。すなわち、

ただし、ｘ_ｐｉ：ピストンの変位
Ｒ_ｃｒ：クランク軸の腕の長さ（０．１３）
ε：Ｒ_ｃｒ／Ｌ_ｃｏｎ
Ｌ_ｃｏｎ：接続の長さ（０．２５）
θ_ｃｒ：クランクの角変位
であり、これらはＳＩ単位系で示している。図３２は、この式（２９）をもとに求めた往復機械のピストンの動解析結果を示す図である。この動解析結果は、図３２に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例７）
図３３は、動解析対象である、ブレーキディスクの概要構成を示す模式図である。このブレーキディスクに固有の運動方程式は、次式（３０）で示される。

ただし、Ｊ_Ｂ：ブレーキディスクの慣性モーメント（３．００）
ｆ_Ｂ：押付力（５００．００）
ｕ_Ｂ：摩擦係数（０．３０）
θ_Ｂ：ブレーキディスクの角変位（初期値＝０．００）
ｄθ_Ｂ（ｔ）／ｄｔ：ブレーキディスクの角速度（初期値＝１３９．６３）
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（３０）から、次式（３１）で示される。すなわち、

によって示される。したがって、式（３０）に、各種パラメータおよび、初期角変位θ_Ｂ（ｔ_０）＝０．００、初期角速度（ｄθ_Ｂ（ｔ）／ｄｔ_０）＝１３９．６３を代入して、角加速度（ｄ^２θ_Ｂ／ｄｔ_０ ^２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角速度（ｄθ_Ｂ（ｔ＋Δｔ）／ｄｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角変位θ_Ｂ（ｔ＋Δｔ）を求める。これら求められた（ｄθ_Ｂ（ｔ＋Δｔ）／ｄｔ）、θ_Ｂ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の角加速度（ｄ^２θ_Ｂ（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。

図３４は、図３３に示したブレーキディスクの動解析結果を示す図である。この動解析結果は、図３４に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例８）
図３５は、動解析対象である、剛体リンクで接続されたマニピュレータ状機器の概要構成を示す模式図である。この剛体リンクで接続されたマニピュレータ状機器に固有の運動方程式は、次式（３２），（３３）で示される。

ただし、Ｌ_ｌｉｎｋ：リンク１，２の長さ（５．０００）
ｇ：重力加速度（９．８０７）
θ_{ｌｉｎｋ１}：リンク１の角変位（初期値＝０．０８７）
ｄθ_{ｌｉｎｋ１}／ｄｔ：リンク１の角速度（初期値＝０．０００）
θ_{ｌｉｎｋ２}：リンク２の角変位（初期値＝−０．０８７）
ｄθ_{ｌｉｎｋ２}／ｄｔ：リンク２の角速度（初期値＝０．０００）
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（３２），（３３）から、次式（３４），（３５）で示される。すなわち、

によって示される。したがって、式（３４），（３５）に、各種パラメータおよび、初期角変位θ_{ｌｉｎｋ１}（ｔ_０）＝０．０８７、初期角変位θ_{ｌｉｎｋ２}（ｔ_０）＝−０．０８７、初期角速度（ｄθ_{ｌｉｎｋ１}（ｔ）／ｄｔ_０）＝０．０００、初期角速度（ｄθ_{ｌｉｎｋ２}（ｔ）／ｄｔ_０）＝０．０００を代入して、角加速度（ｄ^２θ_{ｌｉｎｋ１}／ｄｔ_０ ^２）、（ｄ^２θ_{ｌｉｎｋ２}／ｄｔ_０ ^２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角速度（ｄθ_{ｌｉｎｋ１}（ｔ＋Δｔ）／ｄｔ）、（ｄθ_{ｌｉｎｋ２}（ｔ＋Δｔ）／ｄｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の角変位θ_{ｌｉｎｋ１}（ｔ＋Δｔ）、θ_{ｌｉｎｋ２}（ｔ＋Δｔ）を求める。これら求められた（ｄθ_{ｌｉｎｋ１}（ｔ＋Δｔ）／ｄｔ）、（ｄθ_{ｌｉｎｋ２}（ｔ＋Δｔ）／ｄｔ）、θ_ｌｉｎｋ（ｔ＋Δｔ）、θ_{ｌｉｎｋ２}（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の角加速度（ｄ^２θ_{ｌｉｎｋ１}（ｔ＋Δｔ）／ｄｔ^２）、（ｄ^２θ_{ｌｉｎｋ２}（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。

図３６および図３７は、図３５に示した剛体リンクで接続されたマニピュレータ状機器の動解析結果を示す図である。この動解析結果は、図３６および図３７に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

なお、式（３２），（３３）において、角変位θ_{ｌｉｎｋ１}，θ_{ｌｉｎｋ２}が微小な場合、式（３４），（３５）は、次式（３６），（３７）に示すように近似することができる。すなわち、

として近似できる。この式をもとに動解析を行っても図３６，図３７とほぼ同じ結果を得ることができる。この場合、実施の形態２と同様な処理を適用することができる。

（応用例９）
図３８は、動解析対象である、重力等を相互に及ぼし合う粒子あるいは天体の概要構成を示す模式図である。この重力等を相互に及ぼし合う粒子あるいは天体に固有の運動方程式は、次式（３８）〜（４１）で示される。

ただし、ｍ_１ｇ：地球の質量（５．９７×１０^２４）
ｍ_２ｇ：月の質量（７．３５×１０^２２）
Ｇ：重力定数（６．６７×１０^−１１）
α_ｘ１：地球のｘ方向加速度
α_ｙ１：地球のｙ方向加速度
α_ｘ２：月のｘ方向加速度
α_ｙ２：月のｙ方向加速度
ｒ：地球と月の距離
であり、これらはＳＩ単位系で示している。

なお、式（３８）〜（４１）に示したｒ，ｃｏｓθ，ｓｉｎθは、次式（４２）〜（４４）によって定義される。

ただし、ｘ_１：地球のｘ方向変位（初期値＝０．００×１０^００）
ｘ_２：地球のｙ方向変位（初期値＝０．００×１０^００）
ｙ_１：月のｘ方向変位（初期値＝３．８４×^０８）
ｙ_２：月のｙ方向変位（初期値＝０．００×１０^００）
ｖ_ｘ１：地球のｘ方向速度（初期値＝０．００×１０^００）
ｖ_ｙ１：地球のｙ方向速度（初期値＝０．００×１０^００）
ｖ_ｘ２：月のｘ方向速度（初期値＝０．００×１０^００）
ｖ_ｙ２：月のｙ方向速度（初期値＝１．０４×１０^０３）

この場合、解析方程式は、式（３８）〜（４１）から、次式（４５）〜（４８）で示される。すなわち、

によって示される。したがって、式（４５）〜（４８）に、各種パラメータおよび、地球のｘ方向初期変位ｘ_１（ｔ_０）＝０．００×１０^００、地球のｙ方向初期変位ｙ_１（ｔ_０）＝０．００×１０^００、月のｘ方向初期変位ｘ_２（ｔ_０）＝３．８４×１０^０８、月のｙ方向初期変位ｙ_２（ｔ_０）＝０．００×１０^００を代入して、加速度α_ｘ１，α_ｙ１，α_ｘ２，α_ｙ２を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の速度ｖ_ｘ１（ｔ＋Δｔ），ｖ_ｙ１（ｔ＋Δｔ），ｖ_ｘ２（ｔ＋Δｔ），ｖ_ｙ２（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の変位ｘ_１（ｔ＋Δｔ），ｘ_２（ｔ＋Δｔ），ｙ_１（ｔ＋Δｔ），ｙ_２（ｔ＋Δｔ）を求める。これら求められたｖ_ｘ１（ｔ＋Δｔ），ｖ_ｙ１（ｔ＋Δｔ），ｖ_ｘ２（ｔ＋Δｔ），ｖ_ｙ２（ｔ＋Δｔ），ｘ_１（ｔ＋Δｔ），ｘ_２（ｔ＋Δｔ），ｙ_１（ｔ＋Δｔ），ｙ_２（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の加速度α_ｘ１（ｔ＋Δｔ），α_ｙ１（ｔ＋Δｔ），α_ｘ２（ｔ＋Δｔ），α_ｙ２（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。

図３９は、図３８に示した地球のｘ方向変位の動解析結果を示す図である。この動解析結果は、図３９に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例１０）
この発明では、三体問題についても精度の高い数値解析処理を短時間に行うことができる。三体問題とは、文字通り、３つの物体があるとき、万有引力を用いてそれぞれの物体の振る舞いを求める問題であり、ここでは、太陽、地球、および月のそれぞれに対する動解析を行う。

図４０は、三体問題の動解析対象である太陽、地球、および月の位置関係を示す模式図である。なお、太陽１０１，地球１０２，月１０３はそれぞれｘ−ｙ平面（地球の公転面）の同一平面上で運動するものとして解析する。

まず、太陽１０１の運動方程式は、

として表すことができる。ただし、Ｍ_１は太陽１０１の質量であり、ａ_１は太陽１０１の加速度である。なお、ａ_１はベクトルａ_１（ａ_ｘ１，ａ_ｙ１）である。
地球１０２による万有引力ｆ_１２は、

で表される。Ｇは、万有引力定数を示し、Ｍ_２は地球１０２の質量であり、Ｐ_１およびＰ_２はそれぞれ太陽１０１および地球１０２の位置ベクトルＰ_１（ｘ_１，ｙ_１），Ｐ_２（ｘ_２，ｙ_２）である。また、万有引力ｆ_１２は、ベクトルである。なお、｛Ｇ・Ｍ_２・Ｍ_１／（｜Ｐ_２−Ｐ_１｜^２）｝は、万有引力の絶対値を意味し、｛（Ｐ_２−Ｐ_１）／｜Ｐ_２−Ｐ_１｜｝は、太陽１０１から地球１０２への単位ベクトルを意味する。
式（５０）は、次式に整理される。

同様にして、月１０３による万有引力ｆ_１３は、

で表される。ここで、Ｍ_３は月の質量であり、Ｐ_３は月の位置ベクトルＰ_３（ｘ_３，ｙ_３）である。また、万有引力ｆ_１３は、ベクトルである。
ここで、

であるので、式（４９），（５１）〜（５３）からＦ_１，ｆ_１２，ｆ_１３を消去すると、

となる。ここで、｛Ｇ・Ｍ_２／（｜Ｐ_２−Ｐ_１｜^３）｝＊（Ｐ_２−Ｐ_１）は地球１０２の引力による加速度分であり、｛Ｇ・Ｍ_３／（｜Ｐ_３−Ｐ_１｜^３）｝＊（Ｐ_３−Ｐ_１）は、月１０３の引力による加速度分である。

同様にして地球の加速度ａ_２は、

となる。ここで、｛Ｇ・Ｍ_３／（｜Ｐ_３−Ｐ_２｜^３）｝＊（Ｐ_３−Ｐ_２）は月１０３の引力による加速度分であり、｛Ｇ・Ｍ_１／（｜Ｐ_１−Ｐ_２｜^３）｝＊（Ｐ_１−Ｐ_２）は、太陽１０１の引力による加速度分である。なお、ａ_２はベクトルａ_２（ａ_ｘ２，ａ_ｙ２）である。
さらに同様にして月の加速度ａ_３は、

となる。ここで、｛Ｇ・Ｍ_１／（｜Ｐ_１−Ｐ_３｜^３）｝＊（Ｐ_１−Ｐ_３）は太陽１０１の引力による加速度分であり、｛Ｇ・Ｍ_２／（｜Ｐ_２−Ｐ_３｜^３）｝＊（Ｐ_２−Ｐ_３）は、地球１０２の引力による加速度分である。なお、ａ_３はベクトルａ_３（ａ_ｘ３，ａ_ｙ３）である。

上述した式（５４），（５５），（５６）の３つの連立微分方程式を解析微小時間Δｔ＝１時間ごとに変化させて３体問題の数値解析を行うことができるが、各連立微分方程式は、それぞれ２次元であり、結果的に６元の連立微分方程式を同時に解析することになる。ここで、初期値について説明する。まず、初期状態としては太陽１０１，地球１０２，月１０３が一直線に並んだ瞬間を時刻ｔ＝０とし、太陽１０１を原点とし、地球１０２および月１０３の方向を＋ｘ軸方向にとり、地球１０２および月１０３の速度ベクトルを＋ｙ軸方向にとる。また、
太陽１０１と地球１０２との間の初期距離 ：１．４９６×１０^１１（ｍ）
地球１０２と月１０３との間の初期距離 ：３．８４４×１０^８（ｍ）
地球１０２の初期速度 ：２．９７８×１０^４（ｍ／ｓ）
月１０３の初期周期 ：２７．３２１７（日）
月１０３の地球１０２に対する初期相対速度 ：１０２３（ｍ／ｓ）
太陽１０１の質量Ｍ_１ ：１．９８９×１０^３０（ｋｇ）
地球１０２の質量Ｍ_２ ：５．９７４×１０^２４（ｋｇ）
月１０３の質量Ｍ_３ ：７．３４７×１０^２２（ｋｇ）
万有引力定数Ｇ ：６．６７２５９×１０^−１１（Ｎ・ｍ^２／ｋｇ^２）
である。

この三体問題の数値解析は、まず上述した初期値を式（５４），（５５），（５６）に代入して、加速度ａ_ｘ１，ａ_ｙ１，ａ_ｘ２，ａ_ｙ２，ａ_ｘ３，ａ_ｙ３を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝３６００（ｓ））経過後の速度ｖ_ｘ１（ｔ＋Δｔ），ｖ_ｙ１（ｔ＋Δｔ），ｖ_ｘ２（ｔ＋Δｔ），ｖ_ｙ２（ｔ＋Δｔ），ｖ_ｘ３（ｔ＋Δｔ），ｖ_ｙ３（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝３６００（ｓ））経過後の変位ｘ_１（ｔ＋Δｔ），ｙ_１（ｔ＋Δｔ），ｘ_２（ｔ＋Δｔ），ｙ_２（ｔ＋Δｔ），ｘ_３（ｔ＋Δｔ），ｙ_３（ｔ＋Δｔ）を求める。これら求められたｖ_ｘ１（ｔ＋Δｔ），ｖ_ｙ１（ｔ＋Δｔ），ｖ_ｘ２（ｔ＋Δｔ），ｖ_ｙ２（ｔ＋Δｔ），ｖ_ｘ３（ｔ＋Δｔ），ｖ_ｙ３（ｔ＋Δｔ），ｘ_１（ｔ＋Δｔ），ｙ_１（ｔ＋Δｔ），ｘ_２（ｔ＋Δｔ），ｙ_２（ｔ＋Δｔ），ｘ_３（ｔ＋Δｔ），ｙ_３（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝３６００（ｓ））経過後の加速度α_ｘ１（ｔ＋Δｔ），α_ｙ１（ｔ＋Δｔ），α_ｘ２（ｔ＋Δｔ），α_ｙ２（ｔ＋Δｔ），α_ｘ３（ｔ＋Δｔ），α_ｙ３（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。

図４１は、上述した３体問題における地球を中心とした月の軌道の動解析結果を示す図である。図４１に示すように、この動解析結果は、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

なお、図４２は、従来のオイラー法を用いて地球を中心とした月の軌道の動解析結果を示す図である。この図４２に示した動解析結果は、この応用例１０と同じ解析微小時間Δｔ＝３６００ｓで行った結果であり、月の軌道が発散状態となり、高精度な数値解析結果を得ることができていない。

（応用例１１）
図４３は、動解析対象である、遠心ポンプを有するシステムのサージタンクの概要構成を示す模式図である。このサージタンクの水位変動に固有の運動方程式は、次式（５７）で示される。

ただし、ｇ：重力加速度（９．８０７）
ｌ_{ｓｕｒｇｅ}：管路の長さ（１０．０００）
Ａ_{ｓｕｒｇｅ}：管路断面積（０．００８）
Ａ_ｔａｎｋ：サージタンク断面積（０．１００）
Ｑ_{ｓｕｒｇｅ}：流量（０．０１２）
ｈ′（Ｑ_{ｓｕｒｇｅ}）＝２．６７×１０^９Ｑ_{ｓｕｒｇｅ} ^３−１．２×１０^８Ｑ_{ｓｕｒｇｅ} ^２＋１．２×１０^５Ｑ_{ｓｕｒｇｅ}＋４７２．４８７
：（揚程カーブを流量で微分したもの）
Ｖ_{ｓｕｒｇｅ}：圧力損失係数（０．１４４）
ｘ_{ｓｕｒｇｅ}：サージタンクの水位変動（初期値＝０．０２）
ｖ_{ｓｕｒｇｅ}：サージタンクの水位変動速度（初期＝０．００）
α_{ｓｕｒｇｅ}：サージタンクの水位変動加速度
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（５７）から、次式（５８）で示される。すなわち、

によって示される。したがって、式（５８）に、各種パラメータおよび、サージタンクの初期水位変動ｘ_{ｓｕｒｇｅ}（ｔ_０）＝０．０２、サージタンクの初期水位変動速度ｖ_{ｓｕｒｇｅ}（ｔ_０）＝０．００を代入して、サージタンクの水位変動加速度α_{ｓｕｒｇｅ}を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の水位変動速度ｖ_{ｓｕｒｇｅ}（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．０１（ｓ））経過後の水位変動ｘ_{ｓｕｒｇｅ}（ｔ＋Δｔ）を求める。これら求められたｖ_{ｓｕｒｇｅ}（ｔ＋Δｔ），ｘ_{ｓｕｒｇｅ}（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．０１（ｓ））経過後の水位変動加速度α_{ｓｕｒｇｅ}（ｔ＋Δｔ）を求め、上述した処理を所定回数繰り返す。

図４４は、図４３に示したサージタンクの水位変動の動解析結果を示す図である。この動解析結果は、図４４に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

（応用例１２）
図４５は、動解析対象である、ＬＣＲ回路の概要構成を示す模式図である。このＬＣＲ回路の電荷に固有の運動方程式は、次式（５９）で示される。

ただし、ｇ：重力加速度（９．８０７）
Ｌ_ｌｃｒ：インダクタンス（０．２０）
Ｒ_ｌｃｒ：抵抗（５．００）
ｃ_ｌｃｒ：キャパシタンス（０．００２）
Ｅ_ｌｃｒ：電圧（初期値＝５．００）
ｑ ：電荷（初期値＝０．００）
ｉ ：電流（初期値＝０．００）
であり、これらはＳＩ単位系で示している。

この場合、解析方程式は、式（５９）から、次式（６０）で示される。すなわち、

によって示される。したがって、式（６０）に、各種パラメータおよび、初期電荷ｑ（ｔ_０）＝０．００、初期電流ｉ（ｔ_０）＝０．００を代入して、電荷の二階時間微分量（ｄ^２ｑ／ｄｔ^２）を求める。その後、式（２）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の電流ｉ（ｔ＋Δｔ）を求める。さらに、式（３）に示す演算をそれぞれ行って解析微小時間Δｔ（＝０．００１（ｓ））経過後の電荷ｑ（ｔ＋Δｔ）を求める。これら求められたｉ（ｔ＋Δｔ），ｑ（ｔ＋Δｔ）を用いて、解析微小時間Δｔ（＝０．００１（ｓ））経過後の電荷の二階時間微分量（ｄ^２ｑ（ｔ＋Δｔ）／ｄｔ^２）を求め、上述した処理を所定回数繰り返す。

図４６は、図４５に示したＬＣＲ回路の電荷変動の動解析結果を示す図である。この動解析結果は、図４６に示すように、フィードバック処理等からなる補正処理等が数多く組み合わされた精密かつ時間のかかる高精度な動解析装置の出力結果とほぼ同じ動特性が得られている。

ところで、上述した数値解析装置の一例である動解析装置は、あらかじめ用意されたプログラムをパーソナル・コンピュータやワークステーションなどのコンピュータシステムで実行することによって実現することができる。そこで、上述した動解析装置と同様の機能を有する動解析プログラムを実行するコンピュータシステムについて説明する。

図４７は、上述した動解析装置を実現するコンピュータシステムの構成を示すシステム構成図であり、図４８は、図４７に示したコンピュータシステムにおける本体部の構成を示すブロック図である。

図４７に示すように、このコンピュータシステム２００は、本体部２０１と、本体部２０１からの指示によって表示画面２０２ａに画像などの情報を表示するためのディスプレイ２０２と、このコンピュータシステム２００に種々の情報を入力するためのキーボード２０３と、ディスプレイ２０２の表示画面２０２ａ上の任意の位置を指定するためのマウス２０４とを備える。

また、このコンピュータシステム２００における本体部２０１は、図４８に示すように、ＣＰＵ２２１、ＲＡＭ２２２、ＲＯＭ２２３、ハードディスクドライブ（ＨＤＤ）２２４、ＣＤ−ＲＯＭ２０９を受け入れるＣＤ−ＲＯＭドライブ２２５、フレキシブルディスク（ＦＤ）２０８を受け入れるＦＤドライブ２２６、ディスプレイ２０２、キーボード２０３ならびにマウス２０４を接続するＩ／Ｏインターフェース２２７、ローカルエリアネットワークまたはワイドエリアネットワーク（ＬＡＮ／ＷＡＮ）２０６に接続するＬＡＮインターフェース２２８を備える。

さらに、このコンピュータシステム２００には、インターネットなどの公衆回線２０７に接続するためのモデム２０５が接続されるとともに、ＬＡＮインターフェース２２８およびＬＡＮ／ＷＡＮ２０６を介して、他のコンピュータシステム（ＰＣ）２１１、サーバ２１２並びにプリンタ２１３などが接続される。

そして、このコンピュータシステム２００は、所定の記録媒体に記録された動解析プログラムを読み出して実行することで実施の形態１，２に示した動解析装置を実現する。

ここで、所定の記録媒体とは、フレキシブルディスク（ＦＤ）２０８、ＣＤ−ＲＯＭ２０９、ＭＯディスク、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」の他に、コンピュータシステム２００の内外に備えられるハードディスクドライブ（ＨＤＤ）２２４や、ＲＡＭ２２２、ＲＯＭ２２３などの「固定用の物理媒体」、さらに、モデム２０５を介して接続される公衆回線２０７や、他のコンピュータシステム２１１ならびにサーバ２１２が接続されるＬＡＮ／ＷＡＮ２０６などのように、プログラムの送信に際して短期にプログラムを保持する「通信媒体」など、コンピュータシステム２００によって読み取り可能な動解析プログラムを記録する、あらゆる記録媒体を含むものである。

すなわち、動解析プログラムは、上記した「可搬用の物理媒体」、「固定用の物理媒体」、「通信媒体」などの記録媒体に、コンピュータ読み取り可能に記録されるものであり、コンピュータシステム２００は、このような記録媒体から動解析プログラムを読み出して実行することで、実施の形態１，２に示した動解析装置を実現する。

なお、動解析プログラムは、コンピュータシステム２００によって実行されることに限定されるものではなく、他のコンピュータシステム２１１またはサーバ２１２が動解析プログラムを実行する場合や、これらが協働して動解析プログラムを実行するような場合にも、本発明を同様に適用することができる。

なお、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示した通りに構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。

また、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵおよび当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。

以上のように、この発明は、自然現象や機械運動などにおける種々の運動態様を解析することができる種々の制御分野に用いられる数値解析装置および数値解析プログラムに有用である。

Claims (36)

1. 解析対象の２階微分方程式を数値解析する数値解析装置であって、
   前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、
   解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、
   前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、
   前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、
   前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、
   を備えたことを特徴とする数値解析装置。
2. 前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項１に記載の数値解析装置。
3. 前記設定手段は、前記数値解析の終了条件を設定し、
   前記演算制御手段は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項１に記載の数値解析装置。
4. 前記物理量算出手段によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手段を備えたことを特徴とする請求項１に記載の数値解析装置。
5. 前記出力処理手段は、前記解析微小変数量の増分に対応する、物理量算出手段によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項４に記載の数値解析装置。
6. 前記２階微分方程式は、変数量項を含むことを特徴とする請求項１に記載の数値解析装置。
7. 前記２階微分方程式は、時間項を含むことを特徴とする請求項２に記載の数値解析装置。
8. 前記２階微分方程式は、多元連立方程式であることを特徴とする請求項１に記載の数値解析装置。
9. 表計算プログラムを用いて数値解析を行うことを特徴とする請求項１に記載の数値解析装置。
10. 解析対象の２階微分方程式を数値解析する数値解析装置であって、
    前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手段と、
    解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手段と、
    前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手段と、
    前記２階微分量算出手段が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手段が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手段と、
    前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手段、前記１階微分量算出手段、および前記物理量算出手段に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手段と、
    を備えたことを特徴とする数値解析装置。
11. 前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項１０に記載の数値解析装置。
12. 前記設定手段は、前記数値解析の終了条件を設定し、
    前記演算制御手段は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項１０に記載の数値解析装置。
13. 前記物理量算出手段によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手段を備えたことを特徴とする請求項１０に記載の数値解析装置。
14. 前記出力処理手段は、前記解析微小変数量の増分に対応する、物理量算出手段によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項１３に記載の数値解析装置。
15. 前記２階微分方程式は、変数量項を含むことを特徴とする請求項１０に記載の数値解析装置。
16. 前記２階微分方程式は、時間項を含むことを特徴とする請求項１１に記載の数値解析装置。
17. 前記２階微分方程式は、多元連立方程式であることを特徴とする請求項１０に記載の数値解析装置。
18. 表計算プログラムを用いて数値解析を行うことを特徴とする請求項１０に記載の数値解析装置。
19. 解析対象の２階微分方程式を数値解析する数値解析プログラムであって、
    前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、
    解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、
    前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、
    前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の物理量を加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、
    前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、
    をコンピュータに実行させることを特徴とする数値解析プログラム。
20. 前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項１９に記載の数値解析プログラム。
21. 前記設定手順は、前記数値解析の終了条件を設定し、
    前記演算制御手順は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項１９に記載の数値解析プログラム。
22. 前記物理量算出手順によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手順を備えたことを特徴とする請求項１９に記載の数値解析プログラム。
23. 前記出力処理手順は、前記解析微小変数量の増分に対応する、物理量算出手順によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項２２に記載の数値解析プログラム。
24. 前記２階微分方程式は、変数量項を含むことを特徴とする請求項１９に記載の数値解析プログラム。
25. 前記２階微分方程式は、時間項を含むことを特徴とする請求項２０に記載の数値解析プログラム。
26. 前記２階微分方程式は、多元連立方程式であることを特徴とする請求項１９に記載の数値解析プログラム。
27. 表計算プログラムを用いて数値解析を行うことを特徴とする請求項１９に記載の数値解析プログラム。
28. 解析対象の２階微分方程式を数値解析する数値解析プログラムであって、
    前記２階微分方程式の初期条件を含むパラメータ設定および前記数値解析の解析微小変数量を設定する設定手順と、
    解析開始変数値である基準変数値における前記２階微分方程式の２階微分量を算出する２階微分量算出手順と、
    前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量を乗算し、この乗算した値に前記基準変数値における前記２階微分方程式の１階微分量を加算し、この加算した値を前記解析微小変数量増分後の１階微分量として算出する１階微分量算出手順と、
    前記２階微分量算出手順が算出した２階微分量に前記解析微小変数量の２乗値を乗算し、前記１階微分量算出手順が算出した１階微分量に前記解析微小変数量を乗算し、これらの乗算した値と前記基準変数値における前記２階微分方程式の物理量とを加算し、この加算した値を前記解析微小変数量増分後の物理量として算出する物理量算出手順と、
    前記基準変数値に前記解析微小変数量を加算した変数値を新たな基準変数値として設定し、２階微分量算出手順、前記１階微分量算出手順、および前記物理量算出手順に対してそれぞれ設定された基準変数値を用いて前記２階微分量、前記１階微分量、および前記物理量を算出させる制御を繰り返し行う演算制御手順と、
    をコンピュータに実行させることを特徴とする数値解析プログラム。
29. 前記解析微小変数量および前記基準変数値は、それぞれ解析微小時間および基準時刻であり、前記１階微分量および前記２階微分量は、それぞれ１階時間微分量および２階時間微分量であることを特徴とする請求項２８に記載の数値解析プログラム。
30. 前記設定手順は、前記数値解析の終了条件を設定し、
    前記演算制御手順は、前記数値解析の終了条件を満足する場合に前記数値解析を終了させる制御を行うことを特徴とする請求項２８に記載の数値解析プログラム。
31. 前記物理量算出手順によって算出された各基準変数値毎の物理量を、前記解析微小変数量の増分に対応させて出力する出力処理手順を備えたことを特徴とする請求項２８に記載の数値解析プログラム。
32. 前記出力処理手順は、前記解析微小変数量の増分に対応する、物理量算出手順によって算出された各基準変数値毎の物理量をグラフ化処理することを特徴とする請求項３１に記載の数値解析プログラム。
33. 前記２階微分方程式は、変数量項を含むことを特徴とする請求項２８に記載の数値解析プログラム。
34. 前記２階微分方程式は、時間項を含むことを特徴とする請求項２９に記載の数値解析プログラム。
35. 前記２階微分方程式は、多元連立方程式であることを特徴とする請求項２８に記載の数値解析プログラム。
36. 表計算プログラムを用いて数値解析を行うことを特徴とする請求項２８に記載の数値解析プログラム。

Images