JP5257962B1 - 高速演算装置、高速演算プログラム及び高速演算プログラムを記録した記録媒体、機器制御システム、並びにシミュレーションシステム - Google Patents

Abstract

分数階微積分演算を行う高速演算装置であって、関数ｆ（τ）の積分区間［ａ，ｔ］を、区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割し、区間［ａ，ｃ１］においては、計算してメモリに記憶しておいたｔに依存しない重み付き積分値を用いて微積分演算を行い、区間［ｃ１，ｔ］について計算した微積分演算の結果と合算することにより、微積分演算結果を得る。ｔが更新されるのにしたがって、新たな分割点ｃ２を設定し、区間［ｃ１，ｔ］における計算結果に基づいて区間［ａ，ｃ２］の重み付き積分値を算出し、メモリに記憶する。

Description

本発明は、分数階微積分値を計算する高速演算装置、高速演算プログラム及び高速演算プログラムを記録した記録媒体、この高速演算装置を備えた機器制御システム、並びにこの高速演算装置を備えたシミュレーションシステムに関する。

工学設計や現象解析、装置の制御などの作業において、コンピュータを用いた数値計算やシミュレーション、さらにはそれを応用した製品設計が多くの分野で盛んに実践されている。特に、近年、コンピュータの飛躍的な性能向上により、従来の単純化したモデル式から、より精度に優れたモデル式でのシミュレーションや制御方法が研究されている。その一手法として、従来の整数階の微積分要素でのモデルに代えて、非整数階である分数階微積分要素を含むモデルの適用が有効であり、例えば、振動現象などの挙動解析や、信号処理、画像処理など多くの分野で研究がされている。

分数階微分は、Ｎｅｗｔｏｎ力学で使用される整数階微分を一般化したものであり、現在では工学分野のみならず多くの自然科学分野でも使用されつつある。しかし、分数階微分は、整数階微分に比べて解析が困難なことから、一部では１以上の導関数を含む整数階微分モデルに変換して解析する試みも行われている。

分数階微分は、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義、及びＣａｐｕｔｏの定義で代表されるように、積分式で記述されるため、コンピュータによる分数階微分値の計算は、積分値から求める手法が用いられている。また、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義やＣａｐｕｔｏの定義のように分数階微分は、分数階積分の形で定義されることから明らかなように、分数階微分の計算手法は、分数階積分値の計算にも適用できる。

積分計算としては、従来から台形則やシンプソン則など公知の方法が使用されている。これら従来の方法は、積分計算に際して積分間隔を細かく分割して計算するため、高い精度の計算値を得るには変数軸の分割点を多くする必要がある。しかし、分割点の数が多くなるほど計算時間が増加し、計算精度を求めるほど処理時間が膨大となる。例えば、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義における閉区間［ａ，ｔ］についての積分において、変数の点（変数が時間の場合には時刻）τについての積分は、変数ｔに対する相対値（ｔ−τ）で表された変数ｔに依存の積分であるため、従来の分数階微分値や分数階積分値の計算では、全ての分割点のｆ（τ）をメモリに記憶しておき、変数ｔが進む毎にメモリのｆ（τ）を用いて更新前の積分計算を行う必要があることから、第ｎステップまでの積分計算ではｎ^２オーダーの計算が必要となる。これは、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義、及びＣａｐｕｔｏの定義など、積分形式を含む場合にも共通する課題である。

このような課題を解消する方法として、整数階微分の計算法であるＮｅｗｍａｒｋ−β法を分数階微分に拡張して応用する方法が提案されている（非特許文献１）。

前述の課題を解消するために他の方法として、変数ｔを時間として現時刻から過去の時間領域の分割幅が広くなるようにベキ時間を適用して計算ステップ数を減じることにより、計算時間の短縮化を試みる提案がされている（非特許文献２）。

那須野洋，清水信行：「非線形分数階微分方程式のための数値積分アルゴリズム」，日本機械学会論文集（Ｃ編）７２巻７２２号（２００６−１０）論文Ｎｏ．０６−１３７ 那須野洋，清水信行：「ベキ時間を用いた分数階微分方程式の数値積分法」，日本機械学会論文集（Ｃ編）７２巻７２４号（２００６−１２）論文Ｎｏ．０６−６０５ 那須野洋，清水信行，安野拓也：「分数階微分で記述される粘弾性体の幾何学的非線形静的・動的モデル」，日本機械学会論文集（Ｃ編）７２巻７１６号（２００６−４）論文Ｎｏ．０５−１３０ Ｒ．Ｄｅｎｇ， Ｐ．Ｄａｖｉｅｓ， Ａ．Ｋ．Ｂａｊａｊ：「Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｆｒａｃｔｉｏｎａｌ Ｄｅｒｉｖａｔｉｖｅｓ ｔｏ Ｍｏｄｅｌｉｎｇ ｔｈｅ Ｑｕａｓｉ−Ｓｔａｔｉｃ Ｒｅｓｐｏｎｃｅ ｏｆ Ｐｏｌｙｕｒｅｔｈａｎｅ Ｆｏａｍ」， Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ＡＳＭＥ ２００３ Ｄｅｓｉｇｎ Ｅｎｇｉｎｅｅｒｉｎｇ Ｔｅｃｈｎｉｃａｌ Ｃｏｎｆｅｒｅｎｃｅｓ， ＶＩＢ−４８３９７ （２００３）

しかしながら、非特許文献１に記載の方法は、処理時間が長く、実用できるものではなかった。

また、非特許文献２に記載の方法によると、幅を広くすると計算精度が低下し、また、時刻ｔが更新される度に全区間の分割点について積分計算を実行しなければならず、依然として実用レベルの精度確保と処理時間の短縮化を実現できるものではなかった。

このように従来方法では、精度と処理時間を両立することが困難であり、実用化において大きな課題となっているが、コストや効率を含めて実用的な計算方法は未だ開発されていない。

分数階微積分が適用、又は期待されている具体例として、電気炉やエアコンなどの温調機能を有する機器の温度制御などの機器制御技術がある。代表的な制御技術としては、現在ＰＩＤ（比例・積分・微分）制御が適用されており、これに分数階微積分を適用したＰＩ^λＤ^μ制御等が、より正確な制御法として研究されているが、分数階微積分値の計算に際して処理時間がかかりすぎて、未だ十分な精度で制御が行えないのが現状である。

また、挙動解析においては、例えば近年研究が進んでいる粘弾性の変形挙動解析について、粘弾性体のゴムやポリウレタンフォームなどでは、試料の荷重試験による周波数応答からその粘弾性を推定する研究が数多く報告されているが、より正確な予測、解析を目的に、分数階微分モデルを用いた粘弾性の推定手法（非特許文献３）が提案され、また、実際に分数階微分モデルを用いてポリウレタンフォームの粘弾性を推定すること（非特許文献４）などが報告されている。

しかしながら、この分数階微積分からなる要素の計算方法は、台形則やシンプソン則やＧｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖ法など様々な方法が提案されているが、精度を求めるほど処理時間が膨大となり、精度と処理時間を両立させる大きな課題となっており、コストや効率を含めて実用的な計算方法は未だ開発されておらず、分数階微積分モデルを用いた解析や製品設計の実用化できない原因となっていた。

このようなことから、高精度かつ短時間で処理可能な分数階微積分値の計算方法が開発されれば、様々な分野での実用化が進み、産業的に大きな貢献をもたらすため、コンピュータを用いた革新的な分数階微積分値の高速演算装置の開発が切望されていた。

従って、本発明の目的は、上記問題点に鑑み、高い精度を維持しつつ短い計算時間で分数階微積分値を計算することができる高速演算装置、高速演算プログラム及び高速演算プログラムを記録した記録媒体、機器制御システム、並びにシミュレーションシステムを提供することにある。

本発明者等は、上記問題点を解決するために、分数階微分をＲｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義（以下ＲＬ定義と称する）、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義（以下ＧＬ定義と称する）、及びＣａｐｕｔｏの定義（以下Ｃ定義と称する）による式を用いて分数階微分値及び／又は分数階積分値を求める際に、ｆ（τ）の区間［ａ，ｔ］の積分計算を行うにあたり、区間［ａ，ｔ］をその区間内の点ｃ（以下分割点と称する）で区間［ａ，ｃ］と区間［ｃ，ｔ］とにまず分割し、区間［ａ，ｃ］の積分を特定の条件の下で、（ｔ−τ）のベキをｔのベキとτのベキの積に変換することによって、従来のｆ（τ）に代えて、変数ｔと独立なτによるｆ（τ）の重み付き積分を用いて、ｆ（τ）の重み付き積分項と、変数ｔのベキに依存する重み付き係数との積の和として計算するように構成した高速演算装置を提供する。これにより、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させておき、その値を読み出して使用することによって、その都度、最初から再計算する必要がないため、さらに、全ての計算ステップに対する関数値ｆ（τ）を記憶部上に記憶しておく必要もなくなるため、極めて高速に分数階微積分値計算を行うことができる。また、分数階微分の定義に分数階積分が含まれていることから、分数階微分の場合と同様に分数階積分値も極めて高速に計算することができるのである。

即ち、本発明によれば、高速演算装置は、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求めるように構成されており、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えている。演算部は、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ_１進む毎に式（１）を満たしながら更新される分割点ｃ_１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する第１の分割設定手段と、

（ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ_１，ｔ］又は区間［ｔ，ｃ_１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］において、変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、分割点ｃ_１を上端とする幅Δ_１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ_１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ_１］の分数階微積分値Ａと区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えている。

変数ｔに関する積分区間が［ａ，ｔ］であって、
第１の分割設定手段は、式（２）の条件を満たす設定値Ｍ及びｍにおいて、ｔ−（Ｍ＋ｍ）ｈ＝ｃ_１となった時点を分割点ｃ_１として設定するものであり、
積分演算手段は、区間［ａ，ｃ_１］の変数ｔの更新に依存しない重み付き積分値の算出を、分割点ｃ_１が設定された時点からｃ_１＝ａとして開始し、以降、変数ｔがｍｈ（＝Δ_１）進みｃ_１＝ｔ−（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｃ_１，ｔ−Ｍｈ］の幅ｈ刻みのｆ（τ）データを用い、ｃ_１＝ｔ−Ｍｈへの更新を伴って、更新されたｃ_１を上端とする幅Δ_１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して記憶部に記憶させる第１の積分演算部を備えている

（ただし、ε_Ｔは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ_１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
ことが好ましい。

変数ｔに関する積分区間が［ｔ，ａ］であって、
第１の分割設定手段は、式（３）の条件を満たす設定値Ｍ及びｍにおいて、ｔ＋（Ｍ＋ｍ）ｈ＝ｃ_１となった時点を分割点ｃ_１として設定するものであり、
積分演算手段は、区間［ｃ_１，ａ］の変数ｔの更新に依存しない重み付き積分値の算出を、分割点ｃ_１が設定された時点からｃ_１＝ａとして開始し、以降、変数ｔがｍｈ（＝Δ_１）進みｃ_１＝ｔ＋（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｔ＋Ｍｈ，ｃ_１］の幅ｈ刻みのｆ（τ）データを用い、ｃ_１＝ｔ＋Ｍｈへの更新を伴って、更新されたｃ_１を上端とする幅Δ_１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して記憶部に記憶させる第１の積分演算部を備えている

（ただし、ε_Ｔは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ_１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
ことも好ましい。

第１の分割設定手段によって、式（４）の条件を満たす設定値Ｍ及びｍにおいて、変数ｔのｈ刻みの更新に伴い分割点ｃ_１がｃ_１＝ｔ−Ｍｈとして更新され、さらに、ｃ_１＝ａを初期値として、区間［ａ，ｃ_１］内のｃ_１を上端とする幅Δ’が変数ｔのｈ刻みの更新によってｈからｍｈまで変化するものであって、

（ただし、ε_Ｔは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ_１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
積分演算手段は、幅Δ’がｍｈ未満の場合には、ｃ_１を上端とする幅Δ’の小区間は幅ｈ刻みのｆ（τ）データを用い、変数ｔの更新に依存する積分値の計算を変数ｔの更新毎に繰り返し実行して記憶部に記憶させ、幅Δ’がｍｈとなった場合には、小区間の幅ｈ刻みのｆ（τ）データを用い、分割点ｃ_１を上端とする幅Δ_１（＝ｍｈ）の小区間の変数ｔの更新に依存しない重み付け積分値を算出する重み付き計算を変数ｔの更新毎に繰り返し実行して記憶部に記憶させる第１の積分演算部を備えており、
第２の分数階微積分値演算手段は、幅Δ’がｍｈとなった場合に分割点ｃ_１に対応する重み付き積分値を記憶部から読み出して、忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して算出する分数階微積分値と、変数ｔの更新に依存する幅Δ’がｍｈ未満の場合の小区間の積分値から算出する分数階微積分値とを合算して算出する第１の分数階微積分値演算部を備えていることも好ましい。

積分演算手段は、変数ｔを更新しながら区間［ａ，ｃ_１］の重み付き積分を算出する際に、幅Δ_１の小区間の個数Ｊ_ｓ−１（ｓ＝２）が式（５）で表されるＬ_ｓと等しくなった時点で、ｃ_２＝ｃ_１−（Ｌ_ｓ−ｍ_１）Δ_１である新たな分割点ｃ_２を設定して区間［ａ，ｃ_１］を中区間［ａ，ｃ_２］と中区間［ｃ_２，ｃ_１］とに分割する第２の分割設定手段を具備し、
中区間［ａ，ｃ_２］においては、ｃ_２の設定された時点並びに以降、変数ｔがｍ_１Δ_１進む度に、分割点ｃ_２から直近の分割点ｃ_１方向に連続するｍ_１個の幅Δ_１の小区間を纏めて幅Δ_２の小区間とし、記憶部に記憶されているｍ_１個の各小区間の重み付き積分値を重み付き積分値として算出して記憶部に記憶させる第２の積分演算部を備えており、

第２の分数階微積分値演算手段は、区間［ａ，ｃ_１］の分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ_２，ｃ_１］と、中区間［ａ，ｃ_２］の各小区間の重み付け積分値、をそれぞれ読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。

積分演算手段は、変数ｔを更新しながら区間［ｃ_１，ａ］の重み付き積分を算出する際に、幅Δ_１の小区間の個数Ｊ_ｓ−１（ｓ＝２）が式（６）で表されるＬ_ｓと等しくなった時点で、ｃ_２＝ｃ_１＋（Ｌ_ｓ−ｍ_１）Δ_１である新たな分割点ｃ_２を設定して区間［ｃ_１，ａ］を中区間［ｃ_１，ａ］と中区間［ｃ_１，ｃ_２］とに分割する第２の分割設定手段を具備し、
中区間［ｃ_２，ａ］においては、ｃ_２の設定された時点並びに以降、変数ｔがｍ_１Δ_１進む度に、分割点ｃ_２から直近の分割点ｃ_１方向に連続するｍ_１個の幅Δ_１の小区間を纏めて幅Δ_２の小区間とし、記憶部に記憶されているｍ_１個の各小区間の重み付き積分値を重み付き積分値として算出して記憶部に記憶させる第２の積分演算部を備えており、

第２の分数階微積分値演算手段は、区間［ｃ_１，ａ］の分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ_１，ｃ_２］と、中区間［ｃ_２，ａ］の各小区間の重み付け積分値、をそれぞれ読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。

積分演算手段は、変数ｔを更新しながら区間［ａ，ｃ_２］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃ_ｓ，ｃ_ｓ−１］の幅Δ_ｓ−１の小区間の個数Ｊ_ｓ−１（ｓ≧３）が式（７）で表されるＬ_ｓと等しくなる度に、ｃ_ｓ＝ｃ_ｓ−１−（Ｌ_ｓ−ｍ_１）Δ_ｓ−１である新たな分割点ｃ_ｓによってｃ_ｓ＋１＝ａを初期値とする中区間［ｃ_ｓ＋１，ｃ_ｓ］と［ｃ_ｓ，ｃ_ｓ−１］とを設定する第３の分割手段を具備し、ｃ_ｓが設定された時点及びそれ以降、変数ｔがｍ_１Δ_ｓ−１進む度に、分割点ｃ_ｓから直近の分割点ｃ_ｓ−１方向に連続するｍ_１個のΔ_ｓ−１の小区間を纏めて幅Δ_ｓの小区間とし、記憶部に記憶されているｍ_１個のΔ_ｓ−１の各小区間の重み付き積分値を用いてΔ_ｓの小区間の重み付き積分値として算出して記憶部に記憶させ、さらに、以降変数ｔの更新によってΔ_ｓの小区間の個数がＬ_ｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、

第２の分数階微積分値演算手段は、区間［ａ，ｃ_１］又は分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ_２，ｃ_１］及び中区間［ｃ_３，ｃ_２］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする中区間［ｃ_ｓ＋１，ｃ_ｓ］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。

積分演算手段は、変数ｔを更新しながら区間［ｃ_２，ａ］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃ_ｓ−１，ｃ_ｓ］の幅Δ_ｓ−１の小区間の個数Ｊ_ｓ−１（ｓ≧３）が式（８）で表されるＬ_ｓと等しくなる度に、ｃ_ｓ＝ｃ_ｓ−１＋（Ｌ_ｓ−ｍ_１）Δ_ｓ−１である新たな分割点ｃ_ｓによってｃ_ｓ＋１＝ａを初期値とする中区間［ｃ_ｓ，ｃ_ｓ＋１］と［ｃ_ｓ−１，ｃ_ｓ］とを設定する第３の分割手段を具備し、ｃ_ｓが設定された時点及びそれ以降、変数ｔがｍ_１Δ_ｓ−１進む度に、分割点ｃ_ｓから直近の分割点ｃ_ｓ−１方向に連続するｍ_１個のΔ_ｓ−１の小区間を纏めて幅Δ_ｓの小区間とし、記憶部に記憶されているｍ_１個のΔ_ｓ−１の各小区間の重み付き積分値を用いてΔ_ｓの小区間の重み付き積分値として算出して記憶部に記憶させ、さらに、以降変数ｔの更新によってΔ_ｓの小区間の個数がＬ_ｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、

第２の分数階微積分値演算手段は、区間［ｃ_１，ａ］の分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ_１，ｃ_２］及び中区間［ｃ_２，ｃ_３］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする中区間［ｃ_ｓ，ｃ_ｓ＋１］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。

ｍが２以上であることも好ましい。

ｍ_１が２であることことも好ましい。

忘却関数が式（９）であり、
積分演算手段は、式（９）の忘却関数又は式（９）の分子をあらかじめ計算し、得られたデータを記憶部に記憶させる第３の積分演算部を備えており、
第２の分数階微分値演算手段を実行する

（ただし、ｂは、個々の小区間を識別するための任意の点、αは分数階微積分の積分定義で規定される積分階数、ｇは分数階微積分の積分定義で規定される関数、ｋは展開次数番号である）
ことも好ましい。

第３の積分演算部は、式（９）の分子をあらかじめ計算し、得られたデータを記憶部に記憶させるものであり、
積分演算手段は、さらに、重み付き積分にあらかじめ式（９）中の係数１／（ｋ！Γ（α−ｋ））を乗じた重み付き積分値として算出して記憶部に記憶させる第４の積分演算部を備えており、
第２の分数階微積分値演算手段は、記憶部に記憶されている小区間の前述の係数を乗じた重み付け積分値を読み出し、全ての中区間の各小区間に対応する式（９）の忘却関数の分子を乗じた値を小区間毎に算出し、全小区間について合算することによって区間［ａ，ｃ_１］の分数階微積分値を算出する第３の分数階微積分値演算部を備えている
ことも好ましい。

区間［ａ，ｃ_１］の分数階微積分定義がＲｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義であり、区間［ｃ_１，ｔ］の分数階微積分定義がＣａｐｕｔｏの定義であることも好ましい。

本発明によれば、さらに、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微分要素の変数ｔにおける分数階微分値を求めるコンピュータの高速演算プログラムであって、
コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
コンピュータを、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ_１進む毎に式（１０）を満たしながら更新される分割点ｃ_１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、

（ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ_１，ｔ］又は区間［ｔ，ｃ_１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］において、変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、分割点ｃ_１を上端とする幅Δ_１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ_１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ_１］の分数階微積分値Ａと区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
して機能させる高速演算プログラムが提供される。

本発明によれば、さらにまた、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求めるコンピュータの高速演算プログラムであって、
コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
コンピュータを、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ_１進む毎に式（１１）を満たしながら更新される分割点ｃ_１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、

（ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ_１，ｔ］又は区間［ｔ，ｃ_１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］において、変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、分割点ｃ_１を上端とする幅Δ_１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ_１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ_１］の分数階微積分値Ａと区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
して機能させるための高速演算プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。

本発明によれば、さらに、設定データである入力データと制御対象のシステム出力とを入力するデータ入力部と、データ入力部から得られる入力データとシステム出力との差分を計算する差分演算部と、差分演算部から得られる差分と制御器の分数階微分要素を有する伝達関数とから制御器出力値を計算する制御器出力演算部と、制御器出力演算部から得られる制御器出力値を制御対象に出力する制御器出力用出力部とを備えており、制御器出力用出力部から得られる制御器出力値と制御対象の伝達関数とからシステム出力を出力して制御対象の被制御機器を作動させる制御を、システム出力と次の設定データとから繰り返し実行するように構成した機器制御システムであって、
制御器出力演算部は、制御器の伝達関数における分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
高速演算装置の演算部は、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ_１進む毎に式（１２）を満たしながら更新される分割点ｃ_１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、

（ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ_１，ｔ］又は区間［ｔ，ｃ_１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］において、変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、分割点ｃ_１を上端とする幅Δ_１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ_１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ_１］の分数階微積分値Ａと区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えている機器制御システムが提供される。

本発明によれば、さらに、物理現象が時間を変数とした分数階微積分要素を含む分数階微積分式を用い、分数階微積分式が積分式で表現されるようにモデル化し、分数階微積分要素から分数階微積分値を計算する第１の演算手段と、分数階微積分値を用いて運動方程式の解を計算する第２の演算手段と、物理現象の測定データと運動方程式の経時的な解とを用いて物理現象の経時的な特性値を計算する第３の演算手段とを備え、物理現象の挙動を解析するシミュレーションシステムであって、
第１の演算手段は、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
高速演算装置の演算部は、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ_１進む毎に式（１３）を満たしながら更新される分割点ｃ_１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、

（ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ_１，ｔ］又は区間［ｔ，ｃ_１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］において、変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、分割点ｃ_１を上端とする幅Δ_１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ_１］又は区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ_１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ_１］の分数階微積分値Ａと区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えているシミュレーションシステムが提供される。

本発明の高速演算装置による分数階微積分値の計算手法によれば、分数階微分の定義式中の積分計算で記憶部に記憶するのはｆ（τ）ではなく、変数ｔに依存しない重み付き積分Ｉ_ｓ，ｋｆ（τ）であるので、記憶部に記憶させるデータ数が低減されて記憶部容量を小さくでき、変数ｔが更新される毎に全積分区間の積分を再計算する必要がなく、その重み付き積分値を読み出して使用すればよいため、また、計算ステップも低減されるので計算時間を大幅に短縮した高速化を行うことができる。さらに、分数階微分の定義に分数階積分が含まれていることから、分数階微分の場合と同様に分数階積分値も極めて高速に計算することができる。

しかも、その計算誤差も小さいので、高精度かつ高速に分数階微積分値を求めることができる。従って、分数階微積分要素を含んでモデル化された様々な現象、事象の解析について、精度よく、かつ高速に結果が得られるので、例えば、ＰＩＤ制御を組み込んだ装置の制御をはじめ、様々な分野の装置の設計や解析などに活用できる。また、本発明の高速演算プログラム及びこの高速演算プログラムを記録したコンピュータ読み取り可能な記録媒体は、プロセス量のＰＩＤ制御や各種物理現象の経時的挙動のシミュレーション、それらを応用した製品設計などを効率的に実施することができる。

さらに、機器制御システムは、その分数階微積分値計算が上述の高速演算装置によって高速で行われるため、高い精度を確保しつつかつ高速な機器制御を行うことができる。また、シミュレーションシステムも、その分数階微積分値計算が上述の高速演算装置によって高速で行われるため、高い精度を確保しつつかつ高速なシミュレーションを行うことができる。

本発明の第１の実施形態として、高速演算装置のハードウェア構成を概略的に示すブロック図である。 図１の実施形態における高速演算装置の一部動作を説明するフローチャートである。 図１の実施形態における高速演算装置の一部動作を説明するフローチャートである。 図１の実施形態における高速演算装置の一部動作の別形態を説明するフローチャートである。 図１の実施形態においてソフトウェアによって構築される高速演算装置の主要部の構成を概略的に示すブロック図である。 図１の実施形態における高速演算装置において、分数階微分値を計算する場合に、区間［０，ｔ］内のｃ点で、区間［０，ｃ］と区間［ｃ，ｔ］とに分割することを説明するための説明図である。 分数階微分値を計算する場合に、微分変数の置き換えを説明するための説明図である。 分数階微分値を計算する場合の計算ステップ１を説明するための模式図である。 分数階微分値を計算する場合の計算ステップ２を説明するための模式図である。 分数階微分値を計算する場合の計算ステップ３を説明するための模式図である。 分数階微分値を計算する場合に、ｓ＝２としたときの計算ステップ４を説明するための模式図である。 分数階微分値を計算する場合の計算ステップ４を説明するための模式図である。 分割点ｃ_１＝ｔ−Ｍｈに固定する条件で計算するアルゴリズムにより分数階微分値を計算する方法の説明図である。 本発明の第２の実施形態として、ＰＩＤ制御システムのハードウェア構成を概略的に示すブロック図である。 図８の実施形態におけるＰＩＤ制御システムの一部動作を説明するフローチャートである。 図８の実施形態においてソフトウェアによって構築されるＰＩＤ制御システムの主要部の構成を概略的に示すブロック図である。 本発明の第３の実施形態において、粘弾性体の変形挙動シミュレーションを行う場合を説明するためのモデル図である。 本発明の第３の実施形態として、粘弾性解析を行うシミュレーション装置の一部動作を説明するフローチャートである。 図１２の実施形態においてソフトウェアによって構築されるシミュレーション装置の主要部の構成を概略的に示すブロック図である。 本発明の高速演算装置における計算方法と従来の計算方法との効果検証結果を対比するためのグラフである。 本発明の高速演算装置における計算方法と従来の計算方法との効果検証結果を対比するためのグラフである。

（Ｉ）第１の実施形態（高速演算装置）
図１は本発明の第１の実施形態として、分数階微積分を求める高速演算装置のハードウェア構成を概略的に示している。なお、分数階微積分の計算手法は、分数階積分値の計算にも適用できるので、以下分数階微分の場合について説明する。なお、分数階積分の計算の場合には、以下の分数階微分の場合の説明において、「分数階微分」を「分数階積分」と置き換えればよい。

図１に示すように、本実施形態における高速演算装置は、コンピュータ１０と、このコンピュータ１０に接続される外部装置１１とから構成されている。

コンピュータ１０は、バス１０ａを介して互いに接続された中央処理装置（ＣＰＵ）１０ｂと、リードオンリメモリ（ＲＯＭ）１０ｃと、ランダムアクセスメモリ（ＲＡＭ）１０ｄと、ハードディスク駆動装置（ＨＤＤ）１０ｅと、音声処理部１０ｆと、画像処理部１０ｇと、ブルーレイディスク／デジタルバーサタイルディスク／コンパクトディスク（ＢＲ／ＤＶＤ／ＣＤ）に関してデータ等の読み書きを行うディスク駆動装置１０ｈと、入出力インタフェース１０ｉとを主に備えたコンピュータ及びこれを作動させるプログラムから構成される。

音声処理部１０ｆはスピーカ１０ｊに接続されており、画像処理部１０ｇはディスプレイ１０ｋに接続されている。ディスク駆動装置１０ｈはＢＲ／ＤＶＤ／ＣＤが装着可能となっており、入出力インタフェース１０ｉにはキーボード１０ｌ、マウス１０ｍ及びプリンタ１０ｎ等が接続されている。さらに、この入出力インタフェース１０ｉには、外部装置１１が接続されている。

ＣＰＵ１０ｂは、ＲＯＭ１０ｃ又はＨＤＤ１０ｅに記憶されているオペレーションシステム（ＯＳ）やブートプログラム等の基本プログラムに従ってＲＡＭ１０ｄに記憶されているプログラムを実行して本実施形態の処理を行う。また、ＣＰＵ１０ｂは、ＲＯＭ１０ｃ、ＲＡＭ１０ｄ、ＨＤＤ１０ｅ、音声処理部１０ｆ、画像処理部１０ｇ、ディスク駆動装置１０ｈ、及び入出力インタフェース１０ｉの動作を制御する。

ＲＡＭ１０ｄは高速演算装置のメインメモリとして使用され、ＨＤＤ１０ｅやディスク駆動装置１０ｈから転送されたプログラムやデータを一時的に記憶する。また、ＲＡＭ１０ｄは、プログラム実行時の各種データが一時的に記憶されるワークエリアとしても使用される。

ＨＤＤ１０ｅには、プログラム及びデータがあらかじめ記憶されている。また、ＨＤＤ１０ｅは、ＲＡＭ１０ｄのワークエリアとしての容量が不足している場合に、補助的にその各種データを一時的に記憶する。

音声処理部１０ｆは、ＣＰＵ１０ｂの指示に従って必要な音声データを再生するための処理を行い、スピーカ１０ｊへその音声信号を出力する。

画像処理部１０ｇは、ＣＰＵ１０ｂの指示に従って画像データを生成する。生成された画像データは、ディスプレイ１０ｋに出力される。

ディスク駆動装置１０ｈは、必要な場合は、ＣＰＵ１０ｂの指示に従って、セットされたＢＲ／ＤＶＤ／ＣＤからプログラムやデータを読出し、ＲＡＭ１０ｄへ転送する。また、セットされたＢＲ／ＤＶＤ／ＣＤへプログラムやデータの書き込みをすることも可能である。

入出力インタフェース１０ｉは、キーボード１０ｌ、マウス１０ｍ及びプリンタ１０ｎとＣＰＵ１０ｂ又はＲＡＭ１０ｄとの間のデータのやり取りを制御する。さらに、この入出力インタフェース１０ｉは、外部装置１１とのデータの入出力を制御するように構成されている。

このようなハードウェア構成の高速演算装置において、ＣＰＵ１０ｂは、作動時は、まず、ＲＡＭ１０ｄ内にプログラム記憶領域、データ記憶領域及びワークエリアを確保し、ＨＤＤ１０ｅやディスク駆動装置１０ｈからプログラム及びデータを取り込んで、プログラム記憶領域及びデータ記憶領域に格納する。次いで、このプログラム記憶領域に格納されたプログラムに基づいて、図２ａ〜図２ｃに示す処理を実行する。ＣＰＵ１０ｂがプログラムを実行することによって、図３に概略的に示すごとき構成の主要部を備えた高速演算装置が構築される。

即ち、本実施形態の高速演算装置は、分数階微分を高速で実行可能な演算装置であり、積分下限ａを０とした例として、図３に示すように、分数階微分値を計算するために必要な各種の一次パラメータや、計算対象となる各変数の関数値データを入力する入力部３０と、入力された一次パラメータやデータから二次パラメータを算出する二次パラメータ演算部３１と、区間［ｃ_１，ｔ］（ただし、ｃ_１は分割点）の分数階微分値を算出する第１の分数階微分値演算部３２と、区間［０，ｃ_１］を構成する各小区間Δｓの入力データｆ（τ）に重み付けした変数ｔに依存しない（変数ｔと独立な）τによる重み付き積分（以下、単に重み付き積分と称す）した値を算出する重み付き積分演算部３３と、重み付き積分値を用いて区間［０，ｃ_１］の分数階微分値を算出する第２の分数階微分値演算部３４と、全区間［０，ｔ］の分数階微分値を合算により求める全分数階微分値算出部３５と、必要に応じて設けられ、計算結果を出力する出力部３６と、入力部３０及び各演算部とデータの保存及び／又は読み出し作業がなされる構成となっており、高速解析プログラム、入力データ及び各演算データを記憶する記憶部３７とを備えるように構築される。

ここで、図１のハードウェア構成と対比すると、入力部３０は、入出力インタフェース１０ｉを介してＣＰＵ１０ｂに接続されている外部装置１１及びキーボード１０ｌ、又はディスク駆動装置１０ｈ等から構成され、二次パラメータ演算部３１、第１の分数階微分値演算部３２、重み付き積分演算部３３、第２の分数階微分値演算部３４及び全分数階微分値算出部３５は、ＣＰＵ１０ｂから構成され、出力部３６は入出力インタフェース１０ｉを介してＣＰＵ１０ｂに接続されている外部装置１１及びプリンタ１０ｎ、又はディスク駆動装置１０ｈ等から構成され、記憶部３７は、ＲＡＭ１０ｄ及びＨＤＤ１０ｅ等から構成される。

図２ａ、図２ｂ及び図２ｃは本実施形態における高速演算装置の主要部の動作の一部を説明している。以下、これら図２ａ、図２ｂ及び図２ｃのフローチャートに基づいて説明を行うが、本発明の原理及び技術思想を適用しかつ本発明の効果を得ることができる構成であれば、これらフローチャートに限定されるものではない。

以下、積分区間が［ａ，ｔ］（ここでａ＜ｔ）において、変数ｔが時間であって、積分下限値ａが０であり、ｔの進む方向がａからｔの方向（正の方向という）である場合に、ＲＬ定義の分数階微分値を演算する手順について説明する。この演算は、図２ａに示す以下のステップＳ１〜Ｓ１１を実行することによって行われる。

（ステップＳ１）：分数階微分値計算用一次パラメータの入力（又は設定）
まず、一次パラメータとして、ｈ、ε、ｍ、ｍ_１、Ｍ、Ｋ、Ｓ_ｍａｘ及び計算終点条件が、入力部１０ａを介して高速演算装置に入力される。ただし、ｈは積分計算時の時間ｔの刻み幅、εは許容誤差、ｍは幅Δ_１の小区間に置き換えられる区間幅ｈの小区間（Δ_０）の個数（又は積分ステップの個数）、ｍ_１は小区間Δ_ｓに纏められる小区間Δ_ｓ−１の個数、Ｍは最初の分割点ｃ_１までの区間幅ｈの積分計算する回数（ステップ数）、Ｋはテイラー展開次数、Ｓ_ｍａｘは時間区間［ｃ_ｓ＋１，ｃ_ｓ］の最大ステップ数である。計算終点条件はパラメータを定数値とし、処理するプログラムに応じて、時刻、ステップ数、及びデータ個数等を適宜設定することができる。

（ステップＳ２）：読み込んだ一次パラメータから二次パラメータの計算（又は設定）
入力設定された一次パラメータから分数階微分値計算に必要な二次パラメータとして、Ｌ_０、Ｌ_ｓ、Δ_１、Δ_ｓ、ｍａｘｓｔｅｐが二次パラメータ演算部３１により式（１４）〜式（１８）を用いて計算され、記憶部３７に記憶される。ただし、Ｌ_０は［０，ｃ_１］区間の計算に移行するデータ個数、Ｌ_ｓは第ｓ番目の中区間を構成する幅Δｓの小区間の最大個数、Δ_１は第１番目の中区間を構成する小区間幅、Δ_ｓは第ｓ番目の中区間を構成する小区間幅、ｍａｘｓｔｅｐはタイムステップ数の上限値である。

（ステップＳ３）：初期条件の設定
計算開始時の初期条件として、パラメータｉ、Ｊ_ｓ、及びｃ_ｓに初期値０が与えられる。即ち、ｉ＝０、Ｊ_０＝０，Ｊ_１＝０，．．．，Ｊ_ｓ＝０，．．．，Ｊ_{Ｓｍａｘ-１}＝０、ｃ_１＝０，ｃ_２＝０，．．．，ｃ_ｓ＝０，．．．，ｃ_Ｓｍａｘ＝０。ただし、ｉはタイムステップ数、Ｊ_ｓ（ｓ＝０，１，２，．．．，Ｓ_ｍａｘ−１）は第ｓ番目の中区間のメモリ数、ｃ_ｓ（ｓ＝１，２，．．．，Ｓ_ｍａｘは第ｓ番目の中区間の上端である。

（ステップＳ４〜ステップＳ７）：区間［ｃ_１，ｔ］の分数階微分値の計算
第１の分数階微分値演算部３２により、区間［ｃ_１，ｔ］の分数階微分値が算出される。まず、計算区間［ｃ_１，ｔ］（ｓ＝０番目の中区間）のメモリ数Ｊ_０がＪ_０＝Ｌ_０となるまで（ステップＳ７においてＮｏと判断されている間）は、タイムステップｉとメモリ数Ｊ_０とを１つずつ歩進させながら（ステップＳ４）、時間幅ｈ刻みのデータｆ_ｉ＝ｆ（τ_ｉ）を入力部３０を介して例えば外部装置１１（図１）から読み込み（ステップＳ５）、台形則などの従来法で［ｃ_１，ｔ］の分数階微分値Ａを算出し、算出した分数階微分値Ａを記憶部３７に記憶させる（ステップ６）。この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ１に相当する。なお、ステップＳ５におけるデータｆ（τ_ｉ）の読み込みは、区間［０，ｔ］の全データをあらかじめ入力して記憶部３７に記憶させたものを読み出す方式や、信号制御等のようにリアルタイムで時間更新時のデータを外部装置１１等から順次取り込む方式等、目的に応じて適宜適用することができる。

（ステップＳ８）：区間［０，ｃ_１］の小区分の重み付き積分の計算
メモリ数Ｊ_０がＪ_０＝Ｌ_０となった時点で（ステップＳ７においてＹｅｓと判断されたとき）、区間［０，ｃ_１］の分数階微分値Ｂの計算に移行する。まず、ステップＳ８において、重み付き積分演算部３３により、区間［０，ｃ_１］の分数階微分値Ｂの計算に用いる、区間［０，ｃ_１］を構成する幅Δ_ｓの小区間（ｓ＝１，２，．．．，Ｓ_ｍａｘ−１）の重み付き積分値が算出される。重み付き積分値の計算は、まず、時刻ｔが更新されてＪ_０＝Ｌ_０となる度にｓ＝１である時刻ｔにおけるｃ_１を上端とするΔ_１幅の小区間についてのみ実行され、ｃ_１は重み付き積分の計算が実行される度にｔ側にΔ_１ずつ移動し、移動後の小区間［ｃ_１−Δ_１，ｃ_１］に固有の重み付き積分値として記憶部３７に記憶される。さらに、ｓ≧２以上に移行して区間［０，ｃ_１］内を幅Δ_ｓの小区間（Δ_ｓ＝ｍ_１Δ_ｓ−１）からなる中区間［ｃ_ｓ＋１，ｃ_ｓ］で構成し、前述のｃ_１をｃ_ｓ、前述のΔ_１をΔ_ｓとして同様の処理を実行する。こうして固有の重み付き積分値を有する各小区間［ｃ_ｓ−（ｊ＋１）Δ_ｓ，ｃ_ｓ−ｊΔ_ｓ］同士は、重複することなく時系列で連続して並んで区間［０，ｃ_１］を構成し、各小区間の重み付け積分値は、時刻ｔに依存しないため、一度計算して記憶させておけば、時刻ｔの更新毎に再計算する必要はないため、非常に高速な計算処理を行うことができる。このステップＳ８における計算の詳細については、サブルーチンフローとして後述する。

（ステップＳ９）：区間［０，ｃ_１］の分数階微分値Ｂの計算
記憶部３７に記憶されている区間［０，ｃ_１］を構成する小区間の重み付き積分値を用い、式（１９）又は式（２０）（ｓ＝１固定の場合）によって分数階微分値Ｂが計算される。式（１９）及び式（２０）の内容については後述する。応用例として、あらかじめ入力して記憶しておいた式（１９）又は式（２０）中の右辺中の（ｔ−ｃ_ｓ＋ｊΔ_ｓ）^{−ｑ−ｋ−１}のベキ計算値を読み出しながら式（１９）又は式（２０）の計算を実行してもよい。

（ステップＳ１０）：［０，ｔ］の分数階微分値の計算
全分数階微分値算出部３５により、分数階微分値Ａと分数階微分値Ｂとを合算して［０，ｔ］の分数階微分値（時刻ｔにおける分数階微分値）が算出される。

（ステップＳ１１）：計算終点判断
次いで計算終点判断が行われ、終点条件に達しない場合は、ステップＳ４へ戻り、タイムステップｉとメモリ数Ｊ_０とを１つずつ歩進（時間更新に相当）させながら計算終点条件までステップＳ４〜Ｓ１１を繰り返して実行することにより、目的の時刻ｔの分数階微分値を算出し、終点条件に達した時点（ステップＳ１１においてＹｅｓと判断したとき）にこの処理を終了する。なお、リアルタイム制御のように、分数階微分値計算のフローが制御フローのサブルーチンとして適用される場合には、制御フローの計算終点条件によって分数階微分値計算フローが終了される。

次に、ステップＳ８の処理内容を図２ｂに示すサブルーチンフローを用いて詳細に説明する。

まず、区間［０，ｃ_１］に移行した時点でｓ＝１が与えられ（ステップＳＳ１）、幅Δ_１の小区間［ｃ_１，ｃ_１＋Δ_１］のｍ個のｈ区間のデータｆ（τ）が記憶部３７から読み出され、その重み付き積分値Ｉ_１，ｋｆ（ｃ_１＋Δ_１）が積分上端をｍｈとした式（２１）によって算出され、又はｍ＝２の場合には式（２１）は、式（２２）によって算出される（ステップＳＳ２）。これら式（２１）及び式（２２）の内容については後述する。次いで、ｃ_１＝ｃ_１＋Δ_１として、分割点ｃ_１を現時刻ｔ側にΔ_１だけ移動させ、重み付き積分値Ｉ_１，ｋｆ（ｃ_１＋Δ_１）が重み付き積分値Ｉ_１，ｋｆ（ｃ_１）として記憶部３７に記憶される（ステップＳＳ３）。この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ２に相当する。

次いで、Ｊ_１をＪ_１＋１と歩進させ、合わせて区間［ｃ_１，ｔ］の幅ｈの区間数に相当するメモリ数Ｊ_０を重み付き積分として纏めたｍ個分を減じてＪ_０＝Ｊ_０−ｍとする（ステップＳＳ４）。このステップＳＳ４の処理は、時間更新されて図２ａのステップＳ４に戻ったときに、区間［ｃ_１，ｔ］の幅ｈの区間数に相当するメモリ数Ｊ_０がＭ個からステップＳ６の分数階微分値Ａの計算が再開され、ｍステップ進むとステップＳ８のサブルーチンに移行するように実行するものである。次いで、幅Δ_１の小区間数に相当するメモリ数Ｊ_１がＬ_ｓに達したかが判断され（ステップＳＳ５）、Ｎｏと判断された場合はサブルーチンを完了して図２ａのステップＳ９へ移行する。

さらに時間更新して図２ａのステップＳ４〜Ｓ１１の処理が繰り返し実行される流れの中で、その都度、図２ｂのサブルーチンにおけるステップＳＳ１〜ＳＳ５の処理が繰り返され、各小区間に固有の重み付き積分値が順次記憶部３７に記憶される。このサブルーチンにおける、この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ３に相当する。なお、図２ａのステップＳ９において、式（１９）によって分数階微分値Ｂが計算される場合においては、現時刻ｔを基準に、ｊ、ｓによって総和計算を行うため、幅Δ_１の小区間の重み付き積分値はｃ_１側から最新の重み付け積分値をｊ＝０番目データとして、過去に向かって順次１番目〜（Ｊ_１−１）番目まで識別ラベルを時間更新毎に貼り替える処理（代入計算）をして記憶部３７に記憶される。このラベルの貼り替え処理には任意の方法が適用できる。

図２ｂに示すサブルーチンのステップＳＳ５において、Ｊ_１＝Ｌ_ｓであると判断された時点で（Ｙｅｓと判断したとき）、（Ｊ_１−１）番目から（Ｊ_１−ｍ_１）番目の幅Δ_１の小区間を幅Δ_２の小区間に纏める処理に移行するために、ステップＳＳ６へ進む。

ステップＳＳ６においては、分点拡大処理を実行するかどうか判断し、実行する場合（Ｙｅｓと判断したとき）は、次のステップＳＳ７へ進む。ステップＳＳ７においては、ｓ＝ｓ＋１（＝２）として分点ｃ_２が設定され、次いで、中区間［ｃ_２，ｃ_１］の（Ｊ_１−１）番目から（Ｊ₁−ｍ_１）番目の幅Δ_１の小区間を纏めた区間［ｃ_２，ｃ_２＋Δ_２］（ここでΔ_２＝ｍ_１Δ_１）の重み付き積分値が、既に計算、記憶された上述の（Ｊ_１−１）番目から（Ｊ₁−ｍ_１）番目の小区間の重み付け積分値を用いて、ｓ＝２とした式（２３）によって算出され、ｍ_１＝２の場合には式（２４）を用いて算出される（ステップＳＳ８）。これら式（２３）及び式（２４）の内容については後述する。なお、以上はｓ＝２として説明したが、本ステップでは、第（ｓ−１）中区間を構成する幅Δ_ｓ−１の小区間のうち、（Ｊ_ｓ−１−１）番目から第（Ｊ_ｓ−１−ｍ_１）番目の小区間の積分データを用いて区間［ｃ_ｓ、ｃ_ｓ+Δ_ｓ］の重み付き積分値が式（２３）を用いて算出されるのである。

次いで、分点ｃ_２をΔ_２分だけ現時刻ｔ側に移動（ｃ_２＝ｃ_２＋Δ_２として幅Δ_２の小区間の上端をｃ_２とする処理）し、重み付き積分値Ｉ_２，ｋｆ（ｃ_２＋Δ_２）が重み付き積分値Ｉ_２，ｋｆ（ｃ_２）として記憶部３７に記憶される（ステップＳＳ９）。この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ４に相当する。

その後、区間［ｃ_３，ｃ_２］中の幅Δ_２の小区間の数に相当するメモリ数Ｊ_２をＪ_２＝Ｊ_２＋１と更新すると共に、幅Δ_２に纏められた幅Δ_１の小区間数ｍ_１個分をステップＳＳ４のＪ_１から減じる処理（Ｊ_１＝Ｊ_１−ｍ_１）が実行される（ステップＳＳ１０）。この処理により、次回のサブルーチン実行時に中区間［ｃ_２，ｃ_１］中の幅Δ_１の小区間数Ｊ_１は（Ｌ_ｓ−ｍ_１）個から開始される。次いで、幅Δ_２の小区間数に相当するメモリ数Ｊ_２がＬ_ｓに達したかを判断し（ステップＳＳ１１）、Ｙｅｓの判断であれば、ステップＳＳ７に戻ってｓ＝ｓ＋１としてｓ≧３の場合に移行し、ｓ＝２の場合と同様に、ステップＳＳ７〜ＳＳ１１の処理が繰り返して実行される。

即ち、ステップＳＳ７において分点ｃ_ｓが設定され、次いで第（ｓ−１）中区間を構成する幅Δ_ｓ−１の小区間のうち、（Ｊ_ｓ−１−１）番目から第（Ｊ_ｓ−１−ｍ_１）番目の小区間の積分データを用いて区間［ｃ_ｓ、ｃ_ｓ+Δ_ｓ］の重み付き積分値が式（２３）を用いて算出される（ステップＳＳ８）。次いで、分点ｃ_ｓをΔ_ｓ分だけ現時刻ｔ側に移動（ｃ_ｓ＝ｃ_ｓ＋Δ_ｓとして幅Δ_ｓの小区間の上端をｃ_ｓとする処理）し、重み付き積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ＋Δ_ｓ）が重み付き積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ）として記憶部３７に記憶され（ステップＳＳ９）、その後、区間［ｃ_ｓ＋１，ｃ_ｓ］中の幅Δ_ｓの小区間の数に相当するメモリ数Ｊ_ｓをＪ_ｓ＝Ｊ_ｓ＋１と更新すると共に、幅Δ_ｓに纏められた幅Δ_ｓ−１の小区間数ｍ_１個分をＪ_ｓ−１から減じる処理（Ｊ_ｓ−１＝Ｊ_ｓ−１−ｍ_１）が実行され、次いで、幅Δ_ｓの小区間数に相当するメモリ数Ｊ_ｓがＬ_ｓに達したかの判断（ステップＳＳ１１）処理が繰り返される。

一方、ステップＳＳ１１における判断がＮｏの場合は、このサブルーチンを完了して図２ａのステップＳ９の処理へ移行し、以後、ステップＳ８の実行の度にサブルーチンＳＳ１〜ＳＳ１１が実行される。

このように、図２ｂのサブルーチン処理は、ｓ≧１においてＪ_ｓがＬ_ｓに達した中区間［ｃ_ｓ，ｃ_ｓ−１］について、ｃ_ｓに直近の連続するｍ_１個（ｓ＝１の場合はｍ個）の幅Δ_ｓ−１の小区間を纏めて幅Δ_ｓの小区間の重み付き積分値とする演算を実行する。また、ｓ−１の場合と同様に、幅Δ_ｓの小区間の重み付き積分値（ｓ≧２）は、ｃ_ｓ側から最新の重み付け積分値をｊ＝０番目データとして、もとからあるデータは過去に向かって順次１番目〜（Ｊｓ−１）番目まで識別ラベルを時間更新毎に張り替える処理（代入計算）をして記憶部３７に記憶される。このラベルの張り替え処理には任意の方法が適用できる。

また、図２ａにおけるステップＳ８のサブルーチンフローとして、図２ｃに示す他のサブルーチンフローを適用することができる。図２ｃに示す他のサブルーチンは、図２ｂのステップＳＳ２〜ＳＳ３及びステップＳＳ８〜ＳＳ９に代えて、ステップＮＳ２〜ＮＳ３及びステップＮＳ８〜ＮＳ９の処理を行うものである。即ち、先にｃ_１及びｃ_ｓの移動を含め旧重み付き積分値データの並べ替えを実行した（ステップＮＳ２及びＮＳ８）後に、新Δ_１小区間の重み付き積分Ｉ_１，ｋｆ（ｃ_１）及び新Δ_ｓ小区間の重み付き積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ）（ここでｓ≧２）の計算を実行する（ステップＮＳ３及びＮＳ９）ものであり、その他のステップＮＳ１、ＮＳ４〜ＮＳ７、及びＮＳ１０〜ＮＳ１１は図２ｂのステップＳＳ１、ＳＳ４〜ＳＳ７、及びＳＳ１０〜ＳＳ１１とほぼ同様である。ただし、図２ｃのフローの方が図２ｂのフローに比べて数値計算の工程が少なくて済むため、より高速で処理することができる。

以上、本発明の第１の実施形態として、変数ｔを時間とし、積分下限ａをａ＝０とした高速演算装置について説明したが、本発明は、変数の変化に方向をもつものであれば時間以外の変数からなる分数階微分要素の分数階微分値の計算にも適用でき、積分下限も０に限定されない。また変数の方向も、変数を時間とした上記説明の例とは逆の現在から過去への流れ方向(負の方向という)の場合にも適用できる。時間以外の変数の例としては、分数階の空間微分モデル化した拡散現象の解析などに適用できる。また、後述する制御技術においては制御対象として、時間以外にも、温度、圧力、長さ（距離）、面積、体積、重量などがあり、これらの変化を計算し、適宜電気信号として制御器に入力し、対象物（製品、原料・反応物）の流量、タンク内のレベル、ポンプやコンプレッサーの回転数などを制御することができる。温度、圧力、長さ（距離）、面積、体積、重量のいずれかを変数とする場合は、上述の計算方法において、時間（時刻）の代わりに、これら温度などを用いることで全く同様に計算することができる。本発明の場合、変数が時間であれば、現時刻を基準として過去に遡って分割点を設けて区間を決め計算することにしているが、温度であれば、基準温度に対して低温側（あるいは高温側）に分割点を設けて区間を決め計算することになる。

（ＩＩ）高速演算装置における分数階微分値の計算方法
本実施形態の高速演算装置は、分数階微分要素を含むモデルにおける分数階微分要素を式（２５）のＲＬ定義、式（２６）のＧＬ定義、式（２７）のＣ定義のような積分形式の定義を用いて区間［ａ，ｔ］の分数階微分値を計算する。ただし、これらの式（２５）〜式（２７）において、正の非整数ｑに対して整数ｎ_ｑは、ｑ＜ｎ_ｑ＜ｑ＋１を満たすものである。また、積分方向が逆向き（変数の進みの方向が負）の場合には、式（２５）のＲＬ定義は、式（２８）、同様に式（２７）のＣ定義の場合には式（２９）となり、中区間及び小区間の並びを負の向きに並べ替えるなど、変数列（上述の例では時間列）を反転すれば、上述の変数の進みが正の場合と全く同じ手法で高速演算することができる。

本実施形態においては、特に、
（１）図４に示すように、区間［ａ，ｔ］内のｃ点（分割点）において、区間［ａ，ｃ］と［ｃ，ｔ］とに分割する、
（２）区間［ａ，ｃ］の積分計算において、入力データｆ（τ）に重み付けした変数ｔに依存しない（変数ｔと独立な）τによる重み付き積分（以下、単に重み付き積分と称す）Ｉ_ｓ，ｋｇ（ｂ）として記憶部に記憶する、
（３）この記憶した重み付き積分Ｉ_ｓ，ｋｇ（ｂ）を用いて積分計算する、
ことを基本思想とするものである。これにより、区間［ａ，ｃ］の積分計算で記憶部に記憶するのはｆ（τ）ではなく、変数ｔに依存しない重み付き積分Ｉ_ｓ，ｋｇ（ｂ）であるから、記憶部容量を小さくできると共に、過去の時刻における積分を再計算する必要がないので、計算時間を短縮でき、このときの計算時間は計算ステップ数ｎに対してｎｌｏｇｎで増加するので、従来のｎ^２に比べて劇的に短くできるようにしたものである。なお、図４において、区間［ｃ，ｔ］のｈ刻みの点τは、ｔ側からτ_０＝ｔ、τ_１＝ｔ−ｈ、τ_２＝ｔ−２ｈ、・・・、τ_Ｊ０−１＝ｔ−（Ｊ_０−１）ｈのように全て異なるが、簡易的に区別せずにτとして関数値を全てｆ（τ）と記載している。また、上述した重み付き積分Ｉ_ｓ，ｋｇ（ｂ）は、式（２５）〜式（２７）の各定義における区間［ａ，ｃ］の積分値Ｉ_ｃｆ（ｔ）を一般化した場合の式（３０）における関数ｇ（τ）の重み付け積分であり、τ＝ｂ−ｕ、ｂ＝ｃ_ｓ−ｊΔ_ｓとして式（３１）で定義される。ここで、ＲＬ定義の場合、α＝−ｑ、ｇ（τ）＝ｆ（τ）であり、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義（以下ＧＬ定義）及びＣａｐｕｔｏの定義（以下Ｃ定義）の場合は、α＝ｎ_ｑ−ｑ、ｇ（τ）＝ｆ^（ｎｑ）（τ）である。この一般化については後述する。

さらに、本実施形態においては、区間［ａ，ｃ］をさらに小区間に分割する分割点ｃ_ｓ（ここで、ｓ＝１，２，３・・・Ｓ）を設定したときに、（ｔ−ｃ_ｓ）が十分に大きい条件（分数階微分値を求める変数ｔの点から十分に離れた領域）では、（ｔ−ｃ_ｓ）の大きさ（変数が時間の場合には、現時刻ｔからどのくらい過去の時点か）に応じた範囲毎に既に計算されて記憶部に記憶されている重み付き積分Ｉ_ｓ，ｋｇ（ｂ）を一纏めとした積分値に置き換えて記憶部に記憶し、それを用いて計算しても計算精度を維持することができ、分割点間隔を拡大できるので、重み付き積分Ｉ_ｓ，ｋｇ（ｂ）の数を一層低減することによって計算時間の短縮化と精度の維持が実現できる。

（ＩＩ−１）高速演算化の原理と理論
本発明の分数階微分値の高速演算化の原理と理論について、積分区間が［ａ，ｔ］（ここでａ＜ｔ）において、変数ｔが時間であって、ｔの進む方向がａからｔの方向（正の方向）であるＲＬ定義を例にとり説明する。

式（２５）のＲＬ定義式において、積分する時間区間［ａ，ｔ］を分割点ｃで分割すると、式（２５）は式（３２）に変形することができる。

式（３２）において、第１項の積分の上限ｃと第２項の積分の下限ｃとは共に時間とともに変化するが、第２項の積分の下限ｃに対する微分は、ちょうど第１項の積分の上限ｃの微分を相殺するので、式（３２）では、２つの積分のそれぞれ上限及び下限として固定されているとして計算する。そして、区間［ａ，ｃ］の積分に対応する右辺の第１項の積分に注目すると、この第１項の積分は、前述の通り積分区間の上下限が固定されているので、時刻ｔについての微分を積分の中に入れることが可能であり、第１項の積分をＩ_ｃｆ（ｔ）と定義して式（３３）のように変形できる。

通常、台形則やシンプソン則による数値積分の古典的な手法では、区間［ａ，ｃ］を短い間隔で分割した区間の積分を用いて数値計算されるので、今、時刻ｂ（ｂ∈（ａ，ｃ］）を上限とした小区間幅Δとして、この小区間の積分Ｉ_Δｆ（ｂ）について考えると、式（３３）は式（３４）の形の積分の和に置き換えられる。ここで、ｂ−Δ≧０である。ただし、時刻ｂは、個々の小区間を識別するラベルのような役割を有しており、個々の小区間を識別できれば、任意の点（時刻）に設定できるものであって、上記のｂの定義はその一例であり、これに限定されない。

この小区間の積分の計算について、本発明では、従来の方法では記憶しておいたｆ（τ）を用いて台形則などの近似方法に代えて以下の方法を適用する。即ち、時刻ｂにおける（ｔ−τ）^−１−ｑ項について、式（３５）のようにＴａｙｌａｒ展開する。ここで、ｋは展開次数であり、式（３５）は｜τ−ｂ｜＜ｔ−ｂのときに収束する。

この式（３５）を式（３４）に代入し、τ＝ｂ−ｕとする。ここでτ＝ｂ−ｕは、図５に示すように、区間［ｂ−Δ，ｂ］内の時刻τに関し、区間上限ｂを基準に定義したものであり、ｂを基準として変数ｕによってτが変化するので、変数τを変数ｕに置き換え、ｕの積分区間を［０，Δ］とすれば、式（３６）が導かれる。こうした置き換えは、区間幅Δの特定区間内における時刻τの積分計算をプログラムする場合においても有用である。

こうして、ある特定区間（区間幅Δ）の上限が決まれば、区間幅を積分範囲とする変数ｕを導入した方がプログラム構成でも取り扱いやすくなる。

さらに、ｕ以外の項を積分外に移動して式（３６）を変形し、式（３７）を得る。ここで、Ｉ_Δ，ｋｆ（ｂ）は式（３８）で定義される重み付け積分である。

式（３７）は、式（３３）のＩ_ｃｆ（ｔ）が、式（３８）で定義されるｕの重み付き積分によって計算できることを示している。そして注目すべきは、式（３８）の重み付き積分が時刻ｔに依存しないので、ｂがｔよりも過去において、式（３８）の重み付き積分Ｉ_Δ，ｋｆ（ｂ）を計算して記憶しておけば、時刻ｔが更新された後にも、その重み付き積分Ｉ_Δ，ｋｆ（ｂ）を再計算することなく式（３３）のＩ_ｃｆ（ｔ）の計算に利用できることであり、これが本発明の重要な第１の原理である。

小区間の幅Δは、式（３７）における級数の収束条件として、式（３９）の条件によって制限される。

式（３７）は無限和であるが、数値計算においては近似値として有限の展開次数（ｋ＝Ｋ）まで、式（４０）を用いて計算する。

このときの真値と近似値との誤差は、式（４１）に示す通りである。ここで、‖ｆ（τ）‖は、｜ｆ（τ）｜の上限であり、Ｏ（ ）は、Ｌａｎｄａｕのオーダー記号である。

式（３８）の小区間の幅Δは、式（３９）の条件を満たせば、任意に設定できる。そして、振動波形の解析問題においては、波形の周期よりも広い間隔で積分すると精度が著しく低下する従来の方式と異なり、この小区間の幅Δは、式（３９）の条件においては振動波長よりも広く設定することさえでき、この点が従来法と全く異なる本発明の重要な第２の原理である。

本発明は、上述した第１の原理と第２の原理に基づいた計算方法の発明であり、これにより高い精度と処理時間の短縮を実現しようとするものである。

（ＩＩ−２）分割点間隔の拡大による更なる高速演算化の理論
次に、時間発展によってｔが増加する状況を考える。時間経過とともにある固定された時点ｂにおける時間間隔Δは、ｔ−ｂに比例して（近似の精度を落とさずに）増加することができる。これが本発明の重要な第３の原理である。

時刻τの単位長さあたりの幅Δの小区間数は、εを比例定数としてΔ＝ε（ｔ−ｂ）とおくと、１／ε（ｔ−ｂ）となる。そして、分割された全ての合計は、（１／ε）ｌｏｇ（ｔ−ａ）で表される時刻τの単位長さ当たりの数を統合して計算される。このように、式（３３）の区間［ａ，ｃ］の積分Ｉ_ｃｆ（ｔ）を計算するための重み付き積分Ｉ_Δ，ｋｆ（ｂ）の総数は、ｌｏｇ（ｔ−ａ）で対数的に増加する。

従って、本発明は、前述した第１の原理及び第２の原理に、さらに第３の原理を適用して、飛躍的に処理時間の短縮化を実現することができる。なお、本原理の導出過程の式（３５）においてＴａｙｌａｒ展開を適用したが、他の展開方法に適用しても同様に、本発明の各原理に帰着する。

（分割点間隔の拡大原理を適用した場合の区間［０，ｃ_１］の分数階微分値計算理論）
この第３の原理は、「現時刻ｔから離れるほど（過去に遠ざかるほど）小区間の幅Δを広く設定しても計算精度が維持される」というものであるので、現時刻ｔ、式（３２）の積分下限ａ＝０とした区間［０，ｔ］において、現時刻ｔからの距離によって幅の異なる区間での積分を計算することを考える。

分割点ｃで分割した区間［０，ｃ］区間をさらに中区間［ｃ_ｓ＋１，ｃ_ｓ］（ここでｓ＝１、２、３、・・Ｓ−１）に分割する。ここで、ｃ_１＝ｃであり、ｃ_Ｓをｃ_Ｓ＝０（ｓ＝Ｓ−１のとき）に固定する。なお、ｓ＝０のときには、ｃ_０＝ｔとして、区間［ｃ_１，ｔ］となり、Δ_０＝ｈに相当する。そして、ｔ−Ｊ_０ｈがｃ_１に相当する。

さらに、中区間［ｃ_ｓ＋１，ｃ_ｓ］内をΔ_ｓ＝２^ｓｈの幅に分割した小区間における積分を考える。この幅Δ_ｓはΔ_ｓ＝２Δ_ｓ−１に相当するが、後述するように２倍に限定されず、２以上の整数倍であればよい。ここで、各小区間のΔ_ｓの個数は１個から複数個まで取り得るので、分割点ｃ_ｓ側から添え字ｊとして０、１、２・・・、Ｊ_ｓ−１と番号付けて考える。

添え字ｓが増えることは現時刻ｔから遠ざかることに対応するので、現時刻ｔから遠い区間［ｃ_ｓ＋１，ｃ_ｓ］（中区間という）内のΔ_ｓの幅も広くなる。

収束条件は、式（３９）において、ΔをΔ_ｓと置き換え、ｂをｃ_ｓと置き換えると、式（４２）となり（ここで、εは精度よい結果を得るための基準値）、その条件下において積分Ｉ_ｃｆ（ｔ）は、式（３８）において、ΔをΔ_ｓと置き、ｂをｃ_ｓ−ｊΔ_ｓと置き換えた式（４５）を用いて式（３７）を区間［ａ，ｃ］において和をとった式（４３）によって計算できる。なお、ここで、ｃ_ｓは中区間の右端点、Δ_ｓは小区間の幅であり、式（４３）の右辺の式（４４）の項は、忘却関数である。また、式（４３）では前記時刻ｂは、ｂ＝ｃ_ｓ−ｊΔ_ｓとおいている。前述のｂ＝ｃ_ｓ−ｊΔ_ｓとおいた場合の時刻ｂは、変数ｔの更新に伴う過去の分割点ｃ_ｓの履歴点（時刻）であって、式（４３）は、ｂ＝ｃ_ｓ−ｊΔ_ｓを上端とする小区間の重み付け積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ−ｊΔ_ｓ）と式（４４）の忘却関数との積を計算して、合算することによって分数階微分値が算出されることを示している。また、式（４２）において、ｓ＝１のとき、式（１）及び式（１０）〜式（１３）となる。また、式（４２）の分母は、積分区間が［ａ，ｔ］の場合にはｔ−ｃ_ｓであり、変数の進みが逆(負の方向)である積分区間［ｔ，ａ］の場合にはｃ_ｓ−ｔである。
なお、式（４３）の計算を実行するプログラムは、各小区間の重み付け積分値の識別は、必ずしもｂ＝ｃ_ｓ−ｊΔ_ｓで行う必要はなく、例えば（ｓ，ｊ）の二次元マトリックスデータとして識別処理する手法の他、公知のコーディング手法を適用してもよい。

式（４３）の中括弧｛｝内が中区間［ｃ_ｓ＋１，ｃ_ｓ］の積分値を与えるので、［０，ｃ_１］が複数の中区間［ｃ_ｓ＋１，ｃ_ｓ］で分割されているときには、それらの中区間の積分値を加算すればよいのである。また、式（４３）はｋについて無限和であるため、数値計算においては有限の展開次数（ｋ＝Ｋ）までコンピュータに計算させればよく、近似値として式（４６）によって計算できる。

このときの真値と近似値の誤差は、｜ｆ（ｔ）｜の上限値をＬ、式（４２）の誤差εを用いて式（４７）となる。ここでＯ（ ）は、Ｌａｎｄｏｕのオーダー記号である。

なお、ｓ＝１に限定し、小区間Δ_１の個数の上限を設けない場合には、式（４６）は式（４８）となり、区間［０，ｃ_１］が複数のΔ_１幅からなる小区間で構成される場合には、各小区間の重み付き積分値を合算することで区間［０，ｃ_１］の分数階微分値が計算できることになる。

ここで、Ｉ_１，ｋｆ（ｃ_１−ｊΔ_１）は、ｃ_１−ｊΔ_１＝ｂとして式（４９）で計算できる。

これは、前述した本発明の第１の原理及び第２の原理に基づいた最も基本的な計算方法に相当する。

以上まとめると、本発明は、前述の第１の原理、第２の原理に基づいて、前述の変数ｔが幅Δ_１進む毎に特定条件を満たしながら更新される分割点ｃ_１を設定することにより、区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割し、区間［ｃ_１，ｔ］は従来法で分数階微分値を算出し、
区間［ａ，ｃ_１］においては、変数ｔがΔ_１更新される度に、幅Δ_１の小区間内のｆ（τ）を用いて、小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶し、各小区間について各小区間の重み付け積分値と、各小区間に対応する忘却関数の積を計算したのち、全ての小区間について合算して分数階微分値を算出し、
区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］の各分数階微分値を合算して、全区間［ａ，ｔ］の分数階微分値を算出することに特徴点を有している。
さらに、前述の小区間の重み付け積分値として、第３の原理に基づき区間［０，ｃ_１］をさらに中区間［ｃ_ｓ＋１，ｃ_ｓ］（ここでｓ＝１、２、３、・・Ｓ−１）に分割した場合には、［ｃ_ｓ＋１，ｃ_ｓ］中の小区間Δ_ｓの重み付け積分値は、中区間ｃ_ｓに隣接して連続する中区間［ｃ_ｓ，ｃ_ｓ−１］の小区間Δ_ｓ−１の重み付け積分値を結合して算出したものを適用できることに特徴点を有している。なお、前述の第１から第３の原理やそれらに基づく上記特徴点は、前述したとおり、分数階積分の計算にも適用できるので、分数階微積分値の計算方法の原理及び特徴点といえる。

（ＩＩ−３）ＲＬ定義による分数階微分式以外への拡張
以上は、ＲＬ定義に沿って説明した本発明の第１の原理〜第３の原理であるが、これらの原理は、分数階微分を積分形式で表記した他の定義を適用した場合においても適用できる。例えば、式（２５）〜式（２７）における、区間［ａ，ｃ］の積分値は、前述の通り、式（３０）として一般表記できる。

このとき、式（４９）の積分下端ａ＝０とすると、式（４０）は次の式（５０）の通り一般化され、本発明の第１の原理及び第２の原理が、積分項を含む他の分数階微分の定義にも適用できる。

ここで、式（５０）中のＩ_Δ，ｋｇ（ｂ）は、式（５１）である。

同様に、区間［０，ｃ］の分数階微分値を計算する式（４６）は、次の式（５２）の通り一般化される。

ここで、式（５２）中のＩ_ｓ，ｋｇ（ｃ_ｓ−ｊΔ_ｓ）は、ｂ＝ｃ_ｓ−ｊΔ_ｓとおいた式（３１）である。また、式（５２）中の忘却関数を、ｂ＝ｃ_ｓ−ｊΔ_ｓとおくと式（５３）となる。なお、変数の進みが負の場合には、式（５３）の分子（ｔ−ｂ）^{α−ｋ−１}は、（ｂ−ｔ）^{α−ｋ−１}となる。

さらに、ｓ＝１に限定し、小区間Δ_１の個数の上限を設けない場合である式（４８）及び式（４９）は式（５４）及び式（５５）で一般化される。

式（５２）〜式（５５）において、ＲＬ定義の場合、α＝−ｑ、ｇ（τ）＝ｆ（τ）とし、ＧＬ定義及びＣ定義の場合は、α＝ｎ_ｑ−ｑ、ｇ（τ）＝ｆ^（ｎｑ）（τ）とすれば、各定義における区間［０，ｃ］（又は［０、ｃ_ｓ］）の積分値が得られ、この積分値から分数階微分値が得られ、また、α＞０、ｇ（τ）＝ｆ（τ）とすれば分数階積分の内区間［０，ｃ］における値が計算できる。

なお、変数の進みが負の場合には、変数の進みが正の方向の場合と逆に、変数が進むと変数は小さくなる（例えば、変数ｔを起点とするとｔ、ｔ−ｈ、ｔ−２ｈ、・・・）。そして、分割点ｃ_１はｃ_１＝ｔ＋ｍｈ、分割点ｃｓ（ｓ≧２）はｃ_ｓ＝ｃ_ｓ−１＋（Ｌ_ｓ−ｍ_１）Δ_ｓ−１であり、分数階微分値を算出は、式（５２）Ｉ_ｓ、ｋｇ（ｃ_ｓ−ｊΔ_ｓ）をＩ_ｓ、ｋｇ（ｃ_ｓ＋ｊΔ_ｓ）に、忘却関数の分子の（ｔ−ｃ_ｓ＋ｊΔ_ｓ）^{α−ｋ−１}を（ｔ＋ｃ_ｓ−ｊΔ_ｓ）^{α−ｋ−１}のように変数列（上述の例では時間列）を反転して計算すればよい。

（初期値について）
一般に、分数階微分値計算においては、τ≦ａのｆ（τ）の振る舞いを無視しているために、特にａ＝０の場合の事象を取り扱う場合には、τが０以前の履歴がないので、分数階微分方程式の解の解釈に曖昧さが生じる。

本発明においては、τ≦ａで初期値ｆ（τ）＝０、かつτ＝ａにおいてｆ（ａ）のあらゆる微分の微分係数が０、としている。即ち、ｆ（ａ）＝０、ｆ^（１）（ａ）＝０，・・・，ｆ^{（ｎｑ−１）}（ａ）＝０（τ＝ａ、微分階数ｐ＞０）である。この初期値設定の条件であっても本発明の分数階微分値計算の精度は確保されている。

なお、式（２６）のＧＬ定義には、第２項に初期値項を有するが、式（２６）の第１項の計算における初期値と上述の第２項の初期値項は相殺されるため、式（４９）及び式（５１）には式（２６）の初期値項に起因する項は含まれない。そして、ｆ（ａ）＝０、ｆ（１）（ａ）＝０，・・・，Ｄ^ｎｑ−１ｆ（ａ）＝０のとき、ＧＬ定義及びＣ定義は同じ計算式で計算できる。

（ＩＩ−４）重み付き積分の計算アルゴリズム
次に、式（４６）の重み付き積分Ｉ_ｓ，ｋｆ（ｂ）の計算方法について、４つの計算ステップに分けてさらに詳細に説明する。

ここで、現時刻ｔの近傍の時刻τ＝ｂとし、さらに、式（５６）の条件（ここで、εは精度よい結果を得るための上限値）を満たす最小の０より大きいｍの整数倍であるＭを定義する。

式（５６）中のＭは、区間［ｃ_１，ｔ］におけるΔ_０に等しい時間幅ｈ区間の基準個数である。また、式（５６）中のｍは、幅Δ_１の小区間に置き換えられる区間幅ｈの小区間の個数（又は積分ステップの個数）であり、２以上の整数が好ましく、さらにｍ＝２又はＫの値に依存して、ｍは２以上であると本発明の第１の原理〜第３の原理に基づいて計算が最も効率よくできるので特に好ましい。ｍ＝２以外の例として、Ｋ＝２の場合はｍ＝４が計算効率がよい。以下各状態に沿って図６ａ〜図６ｅを用いて説明する。なお、図６ａ〜図６ｅにおいて、区間［ｃ，ｔ］又は区間［ｃ_１，ｔ］のｈ刻みの点τは、ｔ側からτ_０＝ｔ、τ_１＝ｔ−ｈ、τ_２＝ｔ−２ｈ、・・・、τ_Ｊ０−１＝ｔ−（Ｊ_０−１）ｈのように全て異なるが、簡易的に区別せずにτとして関数値を全てｆ（τ）と記載している。

（計算ステップ１）
計算ステップ１は、現時刻ｔが、ｍｈ／（ｔ−ｍｈ）≧ε（又は、Ｊ_０＜Ｍ＋ｍ、ｍは２以上の整数）を満たし、ｃ_ｓ＝ｃ_１＝０となる区間の計算ステップである（図６ａ参照）。

この区間は、本発明の分割点間隔拡大の手法（第１の原理〜第３の原理）を適用すると誤差が大きくなるので、従来法を用いて等間隔の時間幅ｈで数値計算を行う。従って、この計算ステップにおいては、ｆ（τ）の値は全て記憶部に記憶される。

なお、上述のｍｈはｍｈ＝Δ_１に相当する。以下ｍ＝２の条件にて説明するが、ｍは３以上の整数でも構わない。

（計算ステップ２）
計算ステップ１の状態から時間が経過して、ある現時刻ｔで、初めて２ｈ／（ｔ−２ｈ）＜ε又は、Ｊ_０がＪ_０＝Ｍ＋２（Ｊ_０＝Ｍ＋ｍにおいてｍ＝２としている）を満たした状態となったとき、この時点で、τ＝０、ｈ、２ｈのｆ（τ）を使って式（５７）で積分計算する。これは、式（３８）及び式（４５）において、ｂ＝２ｈ、Δ＝２ｈ（式（４５）においてはΔ_１＝２ｈ）に該当する（図６ｂ参照）。

そして、得られた重み付き積分値は、式（５８）のようにラベルを貼って記憶部に記憶させる（プログラミングにおける代入して保存する作業に相当）。ここで記号（：＝）は、右辺を左辺に置換する操作と定義している。

このとき、Ｊ_１＝１である。そして、区間［ｃ_１，ｔ］の時間幅ｈ（＝Δ_０）による分割点数Ｊ_０は２個減少してＭ個となる。

現時刻ｔが（Ｍ＋２）ｈを超えたとき、区間［ｃ_１，ｔ］の間隔Δ_０のｆ（τ）＝ｆ（ｔ−ｊｈ）（ここでｊ＝０，１・・・Ｊ_０）の値は高速演算装置の記憶部に記憶される。これらのデータは、式（３２）の右辺の第２項の計算に使用する。また、これらのデータは、時間の更新によって、順次、式（４６）中の重み付き積分Ｉ_１，ｋｆ（ｃ_１）の計算にも使用する。

重み付き積分Ｉ_ｓ，ｋｆ（ｃ_ｓ−ｊΔ_ｓ）の計算の道筋は２つあり、第一の筋道として、最初にｓ＝１における幅Δ_１の小区間の重み付き積分Ｉ_１，ｋｆ（ｃ_１−ｊΔ_１）が、区間［ｃ_１，ｔ］のｆ（τ）の値を用いて計算される計算ステップ３と、次に、第二の筋道として、ｓ≧２の場合に、第ｓ中区間中の間隔Δ_ｓにおける重み付き積分Ｉ_ｓ，ｋｆ（ｃ_ｓ−ｊΔ_ｓ）が、第（ｓ−１）中区間中のΔ_ｓ−１の積分値Ｉ_{ｓ−１，ｋ}ｆ（ｃ_ｓ−１−ｊΔ_ｓ−１）を用いて計算される計算ステップ４とである。

（計算ステップ３）
計算ステップ３は、時刻ｔが進んで、２ｈ／（ｔ−ｃ_１−２ｈ）＜ε（ただし、ｃ_１≧０）、又はＪ_０がＪ_０＝Ｍ＋２（Ｊ_０＝Ｍ＋ｍにおいてｍ＝２としている）の条件を満たした状態となったときであって、時刻ｔが進む以前に算出して記憶部に記憶されているｊ＝０，１，２，．．．，Ｊ_１−１個の重み付き積分Ｉ_１，ｋｆ（ｃ_１−ｊΔ_１）に加えて、τ＝ｃ_１，ｃ_１＋ｈ，ｃ_１＋２ｈ（ｊ_０＝Ｍ＋２，Ｍ＋１，Ｍ又はｊ_０＝Ｊ_０，Ｊ_０−１，Ｊ_０−２に対応）におけるｆ（τ）を用いて、ｓ＝１であるΔ_１の重み付き積分Ｉ_１，ｋｆ（ｃ_１＋２ｈ）を式（５９）によって計算し、高速演算装置の記憶部に記憶させる計算ステップである。この式（５９）は、ｂ＝ｃ_１＋２ｈ、Δ_１＝２ｈとしたときの式（４５）の計算に相当する（図６ｃ参照）。

次に、既に第ｓ＝１中区間に記憶されているＩ_１，ｋｆ（ｃ_１−ｊΔ_１）のラベルを式（６０）の第２式のように貼り替え、式（５９）で得られた値をＩ_１，ｋｆ（ｃ_１）の値として式（６０）の第３式のように代入する。ただし、式（６０）の第２式及び第３式の右辺のｃ_１は、第１式によってｃ_１の値を変更（貼り替え）する前のものである。なお、図６ｃにおいては、式（６０）の第２式の右辺の重み付け積分Ｉ_１，ｋｆ（ｃ_１−ｊΔ_１）は、ｂ＝ｃ_１−ｊΔ_１としてＩ_１，ｋｆ（ｂ）と表している。

ここで注目すべきは、重み付き積分及び時刻ｔが進む以前に算出して記憶部に記憶されている重み付き積分は、ラベルのみ貼り替えており再計算しないで、代入記憶のみ実行することである。また、計算ステップ３において時刻ｔがさらに進んで（更新されて）時間経過したときにも、重み付き積分の計算は、更新された［ｃ_１，ｃ_１＋２ｈ］の区間のみ式（５８）で計算し、時間更新前に計算された重み付き積分は、順次式（６０）に倣ってラベルを貼り替えて記憶部に記憶させていくことである。

式（６０）におけるＪ_１及びＪ_０（ｓ＝０，１におけるＪ_ｓに相当）のラベルの貼り替えは、ｓ＝０の幅Δ_０の小区間からｓ＝１の幅Δ_１の小区間への変換の数の変化として反映されている。

（計算ステップ４）
この計算ステップ４は、時間を更新して、現時刻ｔが、Δ_ｓ／（ｔ−ｃ_ｓ−Δ_ｓ）＜ε（又はＪ_ｓ−１がＪ_ｓ−１＝（ｍ_１−１）Ｍ／ｍ＋ｍ_１）の条件を満たした状態となったときであって、ｓ番目の中区間［ｃ_ｓ＋１，ｃ_ｓ］（ここでｓ＝２，３，・・・，Ｓ−１）について、重み付き積分Ｉ_Ｓ，ｋｆ（ｃ_ｓ＋Δ_ｓ）を、既に計算されて記憶部に記憶された（ｓ−１）番目の中区間の幅Δ_ｓ−１小区間のうち、ｃ_ｓから直近の連続するｍ_１個の幅Δ_ｓ−１小区間の重み付き積分値を用いて計算する計算ステップである（図６ｄ及び図６ｅ参照）。ここで、ｍ_１は、２以上の整数であって、後述する通り、ｍ_１＝２が最も計算効率がよいので、以下ｍ＝２、ｍ_１＝２の例で説明するが、ｍ及びｍ_１は３以上の整数でも構わない。なお、図６ｄはｓ＝２の場合であって、式（６０）の第２式の右辺の重み付け積分Ｉ_１，ｋｆ（ｃ_１−ｊΔ_１）は、ｂ＝ｃ_１−ｊΔ_１としてＩ_１，ｋｆ（ｂ）と表している。

重み付き積分Ｉ_ｓ，ｋｆ（ｃ_ｓ＋Δ_ｓ）は、式（６１）で計算する。なお、同式中の２Δ_ｓ−１の２がｍ_１に相当する。

式（６１）の右辺は、式（６４）の通り、次の式（６２）及び式（６３）に示すｓ−１番目の中区間を構成するｊ＝Ｊ_ｓ−１−１番目とＪ_ｓ−１−２番目の小区間Δ_ｓ−１の重み付き積分値を合算して求める。

Ｉ_Ｓ，ｋｆ（ｃ_ｓ＋Δ_ｓ）が算出された後、式（６５）の通りラベルを貼り替えて代入記憶する。次に、既に算出して記憶されている第ｓ中区間を構成する小区間のΔ_ｓ重み付け積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ−ｊΔ_ｓ）のラベルを式（６５）の第２式のように貼り替え、式（６３）で得られた値をＩ_ｓ，ｋｆ（ｃ_ｓ）の値として式（６５）の第３式のように代入する。ただし、式（６５）の第２式及び第３式の右辺のｃ_ｓは、第１式によってｃ_ｓの値を変更（貼り替え）する前のものである。なお、図６ｅにおいては、式（６５）でｓ−１としたときの第２式の右辺の重み付け積分Ｉ_{ｓ−１，ｋ}ｆ（ｃ_ｓ−１−ｊΔ_ｓ−１）は、ｂ＝ｃ_ｓ−１−ｊΔ_ｓ−１としてＩ_{ｓ−１，ｋ}ｆ（ｂ）と表している。

ここで、Δ_ｓ＝２Δ_ｓ−１である。

一方、ｍ_１が３以上の場合、式（６１）及び式（６４）は、それぞれ式（６６）及び式（６７）で表される。なお、式（６７）において、ｍ_１＝２とするとν＝２となり式（６３）となる。

なお、式（６７）は分数階微分定義の一般形として式（６８）で表される。

式（６５）のラベル貼り替え作業によって、ラベルＪ_ｓは１個増加し、Ｊ_ｓ−１はＪ_ｓ−１＝Ｍ／２（Ｊ_ｓ−１＝（ｍ_１−１）Ｍ／ｍにおいてｍ＝２、ｍ_１＝２の場合に相当）となるように２個減少し、また、ｍ_１が３以上の場合は、Ｊ_ｓ−１＝（ｍ_１−１）Ｍ／２となるようにｍ_１個減少する。

もしｓ＝Ｓの場合には、τ＝０を指定するラベルの添え字は、ＳからＳ＋１に置き換わり、ｃ_Ｓ＋１＝０、となる。

このラベル貼り替え作業は、計算ステップ３と同様に、ラベルのみ貼り替えており再計算しないで、代入記憶のみ実行すればよく、時刻ｔがさらに進んで（更新されて）時間経過したときにも、重み付き積分の計算は、更新された［ｃ_ｓ，ｃ_ｓ＋２Δ_ｓ−１］の区間のみ式（６３）で計算し、時間更新前に計算された重み付き積分は、順次、式（６５）に倣ってラベルを貼り替えて記憶部に記憶させていくのである。

以上のように、時間の進行に合わせて、区間［ｃ_１，ｔ］（ｓ＝０の区間に相当）の時間幅ｈ毎のｆ（τ）＝ｆ（ｔ−ｊｈ）のデータを使って、計算ステップ１から計算ステップ４によって区間［０，ｃ_１］中の小区間の重み付き積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ−ｊΔ_ｓ）（ここで、ｊ＝０，１，２，・・・Ｊ_ｓ−１、ｓ＝１，２，３，・・・，Ｓ−１）に置き換えながら記憶部に記憶させ、記憶された重み付き積分値を用いて式（４６）によって、区間［０，ｃ_１］の分数階微分値が算出される。

なお、計算ステップ３及び計算ステップ４において、計算した重み付き積分値は、ラベルを貼り換える手法で時間ｔの進行に対応した時系列データとして識別できるようにしているが、本手法に限定されるものではなく、任意のプログラム手法を適用できることはいうまでもない。

（全区間［０，ｔ］の分数階微分値の計算）
このように計算ステップ１〜４で計算された重み付き積分値を用いて、区間［０，ｃ_１］の分数階微分値は、式（４６）（一般化した表記では式（５２））で計算し、区間［ｃ_１，ｔ］の分数階微分値は計算ステップ１のｆ（τ）で従来の方法で計算して、最後に区間［ｃ_１，ｔ］と区間［０，ｃ_１］の分数階微分値を足し合わせて全区間［０，ｔ］の分数階微分値を得る。

なお、このアルゴリズムから明らかなように、現時刻ｔにおける分数階微分値の計算において、データｆ（τ）は区間［ｃ_１，ｔ］分のみ必要であって、τ≦ｃ_１−ｈのｆ（τ）は重み付き積分に置き換わるので記憶部から削除（又は重み付き積分で上書き）して、記憶部容量をさらに節約することができる。

以上のアルゴリズムにおいて、式（４６）中の右辺中の（ｔ−ｃ_ｓ＋ｊΔ_ｓ）^{−ｑ−ｋ−１}のベキ計算をｓとｊとの二元表に対応するデータ群としてあらかじめ与えておけば（計算データとして計算機に読み込んでおけば）、その分，記憶部容量は必要となるが、計算速度をさらに飛躍的に速くできる。ベキ計算をあらかじめｓとｊとの二元表に対応するデータ群としておくという方法は、単一のプログラムで同一時刻の分数階微分値を多数求める必要のある有限要素解析（ＦＥＭ）や、統計データを求めるため繰り返し計算の必要な場合に、特に有効である。また、ただ一度の計算のみの場合も速度が重要視される計算に対して有効である。さらには、式（５３）の忘却関数自体を先に計算して記憶しておくことも計算速度を向上するために有効である。

また、別のアルゴリズムとして、図７に示すように、現時刻ｔの更新に合わせて分割点ｃ_１＝ｔ−Ｍｈとなるように設定して計算するアルゴリズムとしてもよい。この場合ｓ＝１の最初の区間幅Δ’は、ｈからｍｈの間を変化し、ｍｈのときにΔ’＝Δ_１となる。上述のΔ’がｈから（ｍ−１）ｈの間は、区間［ｃ_１，ｔ］の分数階微分値計算と同様にｈ刻み毎のｆ（τ）を用いて時間に依存する重み付き積分計算として計算する（図７（ｂ）及び（ｄ））。これは、上述の最初の区間幅Δ_１の上端を積分下限として現時刻ｔまでをｈ刻み毎のｆ（τ）を用いて、時間に依存する重み付き積分計算することに等しい。上述のΔ’がｍｈの場合にはΔ_１となって、現時刻ｔの更新に依存しない重み付け積分として、上述したアルゴリズムと実質的に同様に計算される（図７（ｃ）及び（ｅ））。なお、図７ａ〜図７ｅは、図６ａ〜図６ｅと同様に、区間［ｃ，ｔ］又は区間［ｃ_１，ｔ］のｈ刻みの点τは、ｔ側からτ_０＝ｔ、τ_１＝ｔ−ｈ、τ_２＝ｔ−２ｈ、・・・、τ_Ｊ０−１＝ｔ−（Ｊ_０−１）ｈのように全て異なるが、簡易的に区別せずにτとして関数値を全てｆ（τ）と記載している。

また、別のアルゴリズムとして、ｂをｂ_ｓ，ｊ＝ｃ_ｓ−ｊΔ_ｓとおいた式（５２）の重み付き積分値Ｉ_ｓ，ｋｇ（ｃ_ｓ−ｊΔ_ｓ）に代えて、式（５２）中の１／（ｋ！Γ（α−ｋ））を式（５２）に乗じた式（６９）の重み付き積分値Ｉ’_ｓ，ｋｇ（ｃ_ｓ−ｊΔ_ｓ）として記憶し、式（７０）で近似計算をすることによって、さらに高速で演算することができる。同様にｓ＝１の場合に対応する式（５４）及び式（５５）はそれぞれ式（７１）及び式（７２）となる。この方法の場合には、記憶した重み付き積分値を纏めて新たな拡大した幅の小区間の重み付き積分値を算出するステップにおけるｍ_１＝２の場合の式（６２）は、式（４９）の一般化した形式として式（７３）となり、ｍ_１が３以上の場合は、式（７４）から算出できる。なお、式（７４）において、ｍ_１＝２とするとν＝２となり式（７３）となる。

さらに、より計算精度を向上できる方法として、区間［０，ｃ_１］の分数階微分値計算においては、ＲＬ定義を適用し、かつ区間［ｃ_１，ｔ］の分数階微分値計算をＣ定義とすることがより好ましい。これは、ＲＬ定義においては、その積分核は積分区間を分割することによって現時刻ｔ近傍以外では（ｔ−τ）^−ｑ−１とできるため、過去の影響（特に誤差の影響）を小さくでき、一方、Ｃ定義は、現時刻ｔでのｎ_ｑ次までの微分係数が既知であればｔでの分数階微分係数定義式から直接求められるため、両者の特徴を活かして全区間［０，ｔ］の分数階微分値の精度を高めることができるのである。

この方法における微分定義式は式（７５）となり、ｃが一定の場合にもｔ−ｃが一定の場合にも成り立つ。そして、式（７５）の第１項の計算は、これまでに説明したＲＬ定義での分数階微分値計算のアルゴリズムを適用すればよい。

なお、図７のアルゴリズムの場合には、現時刻ｔに依存しない重み付き積分の計算領域をＲＬ定義、現時刻ｔに依存する積分計算の領域にＣ定義を適用すればよい。

ちなみに、式（７５）は以下の通り導出される。まず、式（２５）のＲＬ定義式を分点ｃで分割した式（３２）をａ＝０として、式（７６）の通り変形すると、ｃが定数の場合、式（７６）の右辺第１項は式（７７）となる。一方、式（７６）の第２項は式（７８）の通りＣ定義の微分式を関係付けられ、式（７７）及び式（７８）をそれぞれ式（７６）の第１項及び第２項に代入すると式（７５）が得られる。ｃがｔ−ｃで一定の場合にも導出の途中の変形が異なるが、結果として式（７５）となる。

（ＩＩ−５）精度について
上記の計算アルゴリズムにおける精度を確保するためには、幅Δ_ｓの小区間は、前述の式（４２）によって制限されるが、さらにこの制限条件について詳細説明する。

前述した計算ステップ２及び計算ステップ３から理解されるように、第１番目（ｊ＝０番目）の幅Δ_１の小区間は、式（４２）から式（７９）と書ける。

第１番目の小区間の幅Δ_１はΔ_１／（ｔ−ｃ_１）＜εを満足する。ｃ_ｓ（ここで、ｓ＝１，２，・・・、Ｓ_ｍａｘ−１）の値は、計算時刻ｔ＝ｎｈ（ここでｎはタイムステップ数）によって変化する。しかしながら、ｓ番目の中区間を構成する幅Δ_ｓの小区間の最少個数は、ｍ＝ｍ_１＝２の場合Ｍ／ｍであり、一方、小区間の幅Δ_ｓはΔ_ｓ＝２^ｓｈである。このようにｔ−ｃ_ｓの最小値は、ｍ＝ｍ_１＝２として、式（８０）の通り見積もられる。

Δ_ｓ＝２^ｓｈと式（８０）との比から明らかなように、式（４２）の左辺Δ_ｓ／（ｔ−ｃ_ｓ）の比は、常にεよりも小さい２／Ｍ以下であり、任意の整数ｓにおいて常に式（４２）を満足するので、分割点の間隔を拡張しても精度が確保できるのである。

（ＩＩ−６）最小の計算量となるｍ及びｍ_１条件
前述の通り、本発明の分数階微分値の計算方法によれば、従来法に比べてデータポイント数と分数階微分値計算に必要な総データ数を減らせるため、計算装置に記憶されるメモリ数も少なくて済むため処理時間を短縮できる。ここまではアルゴリズムを理解し易いように簡略的にｍ＝ｍ_１＝２の場合で説明したが、次に記憶部容量と計算量が最小となるｍとｍ_１について説明する。

前述の計算ステップ２及び３において、ｍ≧２の整数ｍ、上述のｍの倍数である整数Ｍにおいて、現時刻ｔ（又はＪ_０）がＪ_０＝Ｍ＋ｍ，ｍ≧２の条件を満たしたとき、ｃ_１＋ｊｈ，ｊ＝０，１，２，・・・，ｍにおける各ｆ（τ）を使ってＩ_１，ｋｆ（ｃ_１＋ｍｈ）が計算される。

また、前述の計算ステップ４においては、ｓ−１番目（ここで、ｓ＝２，・・・、Ｓ_ｍａｘ−１）の区間のＪ_ｓ−１がＪ_ｓ−１＝ｍ_１＋（ｍ_１−１）Ｍ／ｍ，ｍ_１≧２のとき、連続するｍ_１個分の小区間の重み付き積分Ｉ_{ｓ−１，ｋ}ｆ（ｃ_ｓ＋ｊΔ_ｓ−１）（ただし、ｊ＝０，１，２，・・・，ｍ_１、ｋ＝０，１，２，・・・，Ｋ）を用いて、ｓ番目の中区間の新たな小区間の重み付き積分値Ｉ_ｓ，ｋｆ（ｃ_ｓ＋Δ_ｓ）として計算される。

これらの手順では、式（４６）の積分の精度を固定するための式（８１）の条件が課される。この条件は、ｔ−ｃ_ｓが最少となる式（８１）の条件の場合であっても、満足されなければならない。ここで、式（８２）を満たすｓ番目の中区間の小区間の幅Δ_ｓは、式（８３）である。

時刻ｔ＝ｃ_ｓにおける計算に必要なデータポイント数は、以下の通りである。
（Ａ）ｒ番目の区間の最小幅は、（ｍ_１ ^ｒ−ｍ_１ ^ｒ−１）Ｍｈ（ｒ＝１，２，・・・，Ｓ−１）であり、小区間Δｒの幅はΔ_ｒ＝ｍｍ_１ ^ｒ−１ｈである、
（Ｂ）このように、中区間ｒ番目の幅Δ_ｒの小区間の最小の個数は、（ｍ_１−１）Ｍ／ｍであって、ｒには依存しない、
（Ｃ）必要なΔ_ｒの小区間数は、ｒ番目の区間（ｒ＞０）において、（ｍ_１−１）Ｍ／ｍ＋ｍ_１であり、一方、ｓ＝０番目の区間（区間［ｃ_１，ｔ］）では、データの組合せを考慮するとＭ＋ｍ＋１である。
（Ｄ）総データ量をＮ_ｍｅｍとして、展開次数がＫ（Ｋ次）の近似値は、式（８４）で与えられる。

時刻ｔ＝ｃ_ｓに至る分数階微分値の計算の総計算ステップ数ｎは、式（８２）よりｎ＝ｍ_１ ^ｓ−１Ｍであることからｓ−１は式（８５）となる。

式（８５）を式（８４）の（ｓ−１）に代入すると、式（８６）が得られる。

式（８６）中、ε_０＝ｍ／Ｍは固定された値である。ここで、ｍ_１／ｌｏｇ_２ｍ_１と（ｍ_１−１）／ｌｏｇ_２ｍ_１とは、ｍ_１≧２で単調増加する。ｍを固定した場合の総データ量が最小となるｍ_１の数は、ｍ_１＝２において総データ量が最小となる。ｍ_１＝２を式（８６）に代入すると、式（８７）となる。

総データ量が最小となるｍの数は、ｍを実変数として、式（８７）を微分して得た式（８８）の値が０となる条件から、式（８９）に示す値となる。

このように、総データ量Ｎ_ｍｅｍは、ｍ_ｍｉｎに最も近い整数ｍのときに最小となる。記憶部に記憶されるデータ量は、次の通り見積もられる。

今、ｍ＝２^ν＋１、またｍ＝２の場合にＭ＝Ｍ_２とおく。ここで、νは、分割点ｃ_１において、ｍ個の小区間でｆ（τ）を積分した場合、２個の小区間でｆ（τ）を積分するのに比べて減少する中区間の数を意味する。

ここで、Ｍ／ｍ＝Ｍ_２／２と仮定することは精度の点から自然である。このことは、異なるｍ間での積分精度が等しいことを保証している。このことから、Ｍ＝Ｍ_２・２^νとなる。

ｍ＝２の場合に、総データ数（Ｎ_ｍｅｍ）_２は、ｍ_１＝２の式（８４）から式（９０）となる。

一方、任意のｍにおいて分数階微分値計算に必要な総データ数（Ｎ_ｍｅｍ）_ｍは、式（９１）となる。

式（９０）と式（９１）との差は、Ｍ／ｍ＝Ｍ_２／２を用いて、式（９２）で表される。

式（９２）の差の値は、Ｋ（展開次数）、ｍ及びＭ_２／２に依存し、タイムステップには依存しない。

（総データ数の低減）
精度ε_Ｔは、ｍ_１の値によらず、式（９３）の条件で拘束される。

展開次数として、式（９３）を満たす最小の整数Ｋ_εＴとすると、式（９４）の右辺を式（９１）に代入すると、式（９５）の通り、与えられた精度を満たす最小のデータ数が求まる。

式（８６）は、ｍ_１＝２において式（９６）の通りに書き換えられ、総データ数はｎ／Ｍ≫１となるｎに対して対数的に増加する。

これによれば、例えば、Ｋ＝２、ｍ＝２、Ｍ＝４０の条件の場合、分数階微分値計算に必要な総データ数（従来の方法を１として）は、１２５０番目のタイムステップ（ｓ＝５）時点においては約１／３．４、１０２４０番目のタイムステップ（ｓ＝８）時点においては約１／１８、さらに、１６００００番目のタイムステップ（ｓ＝１２）時点においては約１／２００となり、タイムステップが多いほど従来法に比べて総データ量が少なくなるため、記憶部容量が節約されると共に、処理時間が劇的に短縮されるのである。

以上、分数階微分値の計算原理、アルゴリズムを説明してきたが、分数階積分値の計算への適用について説明すると、分数階積分_ａＩ^αｆ（ｔ）は、式（９７）で定義され、忘却関数とｆ(τ)の積の積分として定義される。式（９７）の積分区間を［ａ，ｃ］と［ｃ，ｔ］とに分割したときの区間［ａ，ｃ］の積分は式（３０）と等しい。従って、上記で述べた分数階微分値の高速計算の原理、理論が適用できるので、分数階積分値の高速計算にも有効である。また変数の進みが負の方向の場合（変数を時間とした場合には未来から現在への流れ）にも、中区間及び小区間の並びを負の向きに並べ替えるなど、変数列（上述の例では時間列）を反転すれば、上述した通常の変数の進み方向（変数の進みが正）の場合の分数階微積分を全く同じ手法を適用できる。変数の進みが逆の方向の場合、変数をｔとすると、分数階積分の定義式は式（９８）である。

（ＩＩＩ）第２の実施形態（機器制御システム）
本発明の第２の実施形態として、第１の実施形態における高速演算装置を用いて分数階微分値の高速計算を行う機器制御システムについて説明する。この機器制御システムは、その分数階微分値計算が第１の実施形態の場合と同様に高速で行われるため、高い精度を確保しつつかつ高速な機器制御システムが実現される。

本実施形態の機器制御システムのハードウェア構成は、第１の実施形態の場合と同様に、図１に示す、コンピュータ１０と、このコンピュータ１０に接続される外部装置１１とから構成されている。

このようなハードウェア構成の機器制御システムにおいて、ＣＰＵ１０ｂは、作動時は、まず、ＲＡＭ１０ｄ内にプログラム記憶領域、データ記憶領域及びワークエリアを確保し、ＨＤＤ１０ｅからプログラム及びデータを取り込んで、プログラム記憶領域及びデータ記憶領域に格納する。次いで、このプログラム記憶領域に格納されたプログラムに基づいて、図９、及び図２ａ〜図２ｃに示す処理を実行する。ＣＰＵ１０ｂがプログラムを実行することによって、図１０に概略的に示すごとき構成の主要部を備えた機器制御システムが構築される。

即ち、本実施形態の機器制御システムは、設定データである入力データと制御対象のシステム出力（実測値）とを読み込んで、その差分を計算し、その差分と制御器の分数階微分項を有する伝達関数とから制御器出力値を計算して制御対象に出力し、この制御器出力値と制御対象の伝達関数とからシステム出力を出力して制御対象である被制御機器を作動させると共に、システム出力と次の設定値とから上述の制御フローを繰り返しながら制御するように構成されている。その際、上述の伝達関数は、分数階微分要素の分数階微分式が積分式で定義され、この積分式の分数階微分値が、第１の実施形態における高速演算装置と同様に演算される。このように、制御器の分数階微分項を有する伝達関数に基づいた出力によって機器を制御するものであり、その出力計算に本発明の高速演算装置における分数階微分値計算手法が適用される。以下、フィードバック制御の一つであるＰＩＤ（比例・積分・微分）制御に適用される本実施形態の機器制御システムについて説明する。

（ＰＩＤ制御）
一般に、ＰＩＤ制御は、設定値と実測値との差である偏差に対して、この偏差に比例して入力値を変化させる動作（比例動作又はＰ動作）と、この偏差の積分に比例して入力値を変化させる動作（積分動作又はＩ動作）と、この偏差の微分に比例して入力値を変化させる動作（微分動作又はＤ動作）とを組み合わせて制御する技術である。そして、制御対象が分数階の動的なシステムの場合には、積分動作と微分動作とをそれぞれＩ^λ、Ｄ^μとしたλ階の積分器とμ階の微分器とを有する制御器とすることにより、精度のよい制御が可能となることが知られている。この代表的なものには、Ａ．Ｏｕｓｔａｌｏｕｐ等のＣＲＯＮＥ Ｃｏｎｔｒｏｌｌｅｒ等があり、盛んに研究されている。

このようなＰＩＤ制御は、エアコンでは温調制御に、また化学プラントでは、液面制御、圧力制御、又は反応物供給量制御などに適用されているが、積分動作では、周期的変動があり、遅れが生じやすいことが指摘されている。

ＰＩＤ制御において精度良く制御するには制御時間が重要であり、ＰＩＤの性能向上の一因として、制御時間（分数階微分処理の計算時間）の改善が有効である。

（ＰＩＤ制御システムの構成）
図８は、本実施形態の機器制御システムとして、ＰＩＤ制御システムの一例をブロック図で示している。

ＰＩＤ制御システムでは、まず、設定データである入力データＷ（ｓ）と制御対象のシステム出力Ｙ（ｓ）（実測値に相当）とを読み込んでその差分Ｅ（ｓ）を計算し、その差分Ｅ（ｓ）と制御器の伝達関数Ｇ_Ｃ（ｓ）とから制御器出力値Ｕ（ｓ）を計算して制御対象に出力し、制御対象は制御器出力値Ｕ（ｓ）と制御対象の伝達関数Ｇ（ｓ）とからシステム出力Ｙ（ｓ）を出力して、以降、システム出力Ｙ（ｓ）と次の設定値Ｗ（ｓ）とから上述の制御フローを繰り返しながら制御対象が制御される。なお、ここで、ｓはラプラス演算子であり、上記システムの説明は、時間ｔの関数をラプラス変換したｓ領域で行っている。

制御対象の伝達関数Ｇ（ｓ）は、実際には各制御対象ごとに異なるが、分数階微分系制御においては、Ｇ（ｓ）は式（１０１）のようなモデルに基づいて、式（１０２）のような分数階微分モデルであると仮定するが、実際の制御においては、式（１０２）の計算は不可能なので、結果として観測されるシステム出力Ｙ（ｓ）を上述の制御フローに適用している。よって、制御器において制御器出力値Ｕ（ｓ）を計算する部分に分数階微分値計算を実施する。

制御器の伝達関数Ｇ_Ｃ（ｓ）は、Ｐ動作の比例ゲインＫ_Ｐ、Ｉ動作の積分ゲインＫ_Ｉ、及びＤ動作の微分ゲインＫ_Ｄを用いて表記すると式（９９）となり、時刻ｔの微分方程式として、式（１００）となる。なお以降、実施形態の説明においては、微分又は積分の端点を表す下付き添字は省略する。

この式（１００）の分数階微分項の計算には、第１の実施形態における高速演算装置の分数階微分値計算手法を使用する。これにより、高速かつ精度よく分数階微分制御が可能となる。

本実施形態のＰＩＤ制御システムの一部動作を説明するフロー及びソフトウェアによって構築されるＰＩＤ制御システムの主要部の構成が図９及び図１０に示されている。

図１０に示すように、ＰＩＤ制御システムは、制御パラメータが入力される制御パラメータ入力部１００と、制御対象において目標とする出力の設定値Ｗ（ｓ）及びその実測値Ｙ（ｓ）が入力されるＷ（ｓ）及びＹ（ｓ）入力部１０１と、これら設定値Ｗ（ｓ）及び実測値Ｙ（ｓ）の差分Ｅ（ｓ）を計算する差分Ｅ(ｓ)演算部１０２と、この差分Ｅ（ｓ）と制御器の伝達関数Ｇ_Ｃ（ｓ）とから制御器出力値Ｕ（ｓ）を計算する制御器出力値Ｕ（ｓ）演算部１０３と、この制御器出力値Ｕ（ｓ）を出力するＵ（ｓ）出力部１０４と、記憶部１０５とを備えたコンピュータ１０６、並びに制御対象１０７から構成される。差分Ｅ(ｓ)演算部１０２、制御器出力値Ｕ（ｓ）演算部１０３及びＵ（ｓ）出力部１０４は、記憶部１０５に記憶された入力値又は記憶値を用いて計算処理及び出力処理を実行する。ここで、出力されたＵ（ｓ）は、制御対象１０７に入力され、制御対象１０７はその出力Ｕ（ｓ）に応じて作動して応答実測値Ｙ（ｓ）をフィードバックする。

このうち、制御器出力値Ｕ（ｓ）を計算する制御器出力値Ｕ（ｓ）演算部１０３に、第１の実施形態で述べた高速演算装置が適用されており、本発明の分数階微分値計算アルゴリズムに従って記憶部１０５に記憶された差分の履歴であるデータや重み付き積分値を読み出して、本発明の分数階微分値計算が実行される。

このＰＩＤ制御システムの動作を図９のフローに従って説明すると、まず、制御パラメータの入力を行って制御条件の設定を行う（ステップＳ２１）。これは、第１の実施形態の高速演算装置における一次パラメータの入力（ステップＳ１）、二次パラメータの計算（ステップＳ２）及び初期条件の設定（ステップＳ３）の処理とＰＩＤ制御システムとしての必要な制御条件の設定処理とを行うものである。次いで、設定値Ｗ（ｓ）と実測値Ｙ（ｓ）とが読み込まれる（ステップＳ２２）。次いで、これら設定値Ｗ（ｓ）と実測値Ｙ（ｓ）との差分Ｅ（ｓ）が算出され、記憶装置に記憶される（ステップＳ２３）。

その後、記憶装置に記憶された差分Ｅ（ｓ）の履歴データと、制御器の伝達関数Ｇ_Ｃ（ｓ）との積からなる分数階項を有する式が、第１の実施形態の高速演算装置の場合と同様に計算され、制御器出力値Ｕ（ｓ）が算出される（ステップＳ２４）。算出された制御器出力値Ｕ（ｓ）が制御対象１０７へ出力され（ステップＳ２５）、その後、制御終了の判断が行われる（ステップＳ２６）。終了条件に達しない場合は、ステップＳ２２へ戻り、ステップＳ２２〜Ｓ２６の処理が繰り返して実行される。終了条件に達した場合は（ステップＳ２６においてＹｅｓと判断したとき）にこの処理を終了する。

本発明の分数階微分値計算方法を適用したＰＩＤ制御は、電気炉やエアコンなどの温調機能を有する機器の温度制御、ハードディスクドライブ（ＨＤＤ）のヘッダの位置を合わせる制御等の各種アクチュエータの駆動制御を含む既存のＰＩＤ制御技術に広く適用できる。本実施形態の機器制御システムは、こうしたＰＩＤ制御の他にも、画像処理としてはフラクタル画像処理の他、データ圧縮方法の一つであるフラクタル圧縮の処理などにも適用することができる。

（ＩＶ）第３の実施形態（シミュレーションシステム）
本発明の第３の実施形態として、第１の実施形態における高速演算装置を用いて分数階微分値の高速計算を行うシミュレーションシステムについて説明する。このシミュレーションシステムは、その分数階微分値計算が第１の実施形態の場合と同様に高速で行われるため、高い精度を確保しつつかつ高速なシミュレーションが実現される。

本実施形態のシミュレーションシステムは、分数階微分要素を含む運動方程式でモデル化できる物理現象、例えば振動現象や、ブラウン運動及び拡散現象等、任意の物理現象の解析シミュレーションに応用でき、高い精度を確保しつつ、かつ高速な処理を実現できる。以下、本実施形態のシミュレーションシステムとして、分数階微分でモデル化された粘弾性体の外部応力による形状変形挙動シミュレーションについて説明するが、本発明はこれに限られるものではない。

本実施形態のシミュレーションシステムのハードウェア構成は、第１の実施形態の場合と同様に、図１に示す、コンピュータ１０と、このコンピュータ１０に接続される外部装置１１とから構成されている。

このようなハードウェア構成のシミュレーションシステムにおいて、ＣＰＵ１０ｂは、作動時は、まず、ＲＡＭ１０ｄ内にプログラム記憶領域、データ記憶領域及びワークエリアを確保し、ＨＤＤ１０ｅからプログラム及びデータを取り込んで、プログラム記憶領域及びデータ記憶領域に格納する。次いで、このプログラム記憶領域に格納されたプログラムに基づいて、図１２、及び図２ａ〜図２ｃに示す処理を実行する。ＣＰＵ１０ｂがプログラムを実行することによって、図１３に概略的に示すごとき構成の主要部を備えたシミュレーションシステムが構築される。

粘弾性材料の変形解析を行うシミュレーションとしては、公知の種々のシミュレーション法が適用できる。例えば、有限要素解析法（ＦＥＭ）を用いた分数階微分モデルの粘弾性解析が盛んに研究されているため、ＦＥＭ解析手法を例に説明する。その分数階微分要素の計算に本発明の分数階微分値計算方法が適用される。

（解析モデル）
図１１は粘弾性体の振動シミュレーションを行う場合を説明するためのモデル図である。一例として、粘弾性体に荷重印加される１自由度系の振動子モデルを挙げると、運動方程式は、式（１０３）となる。ただし、式（１０３）において、ｍは質量、ｘ（ｔ）は時間の関数である変位、ｋはバネ定数、ｃは粘弾性係数、Ｐ（ｔ）は時間の関数である外力、Ｄは微分演算子（ｄ／ｄｔ）である。なお以降、実施形態の説明においては、微分又は積分の端点を表す下付き添字は省略する。

一般に、３次元の等温・等方性の粘弾性体の動的挙動を支配する運動方程式は，有限要素により定式化すると、要素単位では式（１０４）で表される。

この式（１０４）において、Ｍは質量マトリックスであり、ｄ（ｔ）は節点変位、Ｆ（ｔ）は要素内力ベクトル、Ｐ（ｔ）は等価節点外力ベクトルである。これらＭ、Ｆ（ｔ）、及びＰ（ｔ）は、形状関数をＮ、変位とひずみに関係するマトリックスをＬ、密度をρとすると、体積力ベクトルｆ及び表面力ベクトルＴ（ｔ）を用いて、式（１０５）で定義される。なお、式（１０５）において、Ωは物体の体積Ｖに関する二重積分、Ｓは物体の表面の面積に関する三重積分、上付き添字Ｔは転置を意味する。

図１１のモデルを３自由度に拡張した場合、例えば、日本機械学会［Ｎｏ．００−６］Ｄｙｎａｍｉｃｓ ａｎｄ Ｄｅｓｉｇｎ Ｃｏｎｆｅｒｅｎｃｅ２０００ ＣＤ−ＲＯＭ論文集［２０００．９．５−８東京］２０７「分数階微分モデル粘弾性体の基礎理論」に開示されているように式（１０５）のＦは、式（１０６）で表され、運動方程式は式（１０７）となる。これは、偏差変形が分数階Ｋｅｌｖｉｎ−Ｖｏｉｇｔ則に，体積変形がＨｏｏｋｅ則に従うモデルである。ただし、式（１０６）及び式（１０７）において、Ｋ’は剛性マトリックス、Ｃ_１は粘弾性減衰マトリックス、τは一般化緩和時間、ｑは分数階の微分階数である。また、式（１０７）において、Ｂは変位歪みマトリックスであり、Ｄ_０及びＤ_１は緩和関数マトリックスであって、Ｄ_０は線形の応力歪みマトリックスに対応し、Ｄ_１は粘弾性の効果を表すマトリックスに対応する。

さらに、式（１０７）の左辺に時刻ｔにおける未知の変位ベクトルｄ（ｔ）、右辺に既知のベクトルがくるように変形すると、式（１０８）が得られる。

ここで、Ｋは式（１０９）で表される一般化剛性行列であり、Ｐ’（ｔ）は式（１１０）で表される時刻ｔにおける既知のベクトル（以下、一般化等価節点力ベクトルという）である。

ここで、Ｓ（ｔ_ｉ−１）は、式（１０８）の左辺中の時刻ｔ_ｉ−１における既知量を纏めた項であり、Ｎｅｗｍａｒｋ−β法の数値積分を適用すると、ｔ_０を初期値として式（１１１）で表され、節点変位ベクトルｄ（ｔ）の分数階微分の記憶積分項を含んでいる。ただし、式（１１１）中のα、βは、Ｎｅｗｍａｒｋ−β法のパラメータである。

ＦＥＭ手法による解析では、現時刻のＰ’（ｔ）と剛性行列Ｋとから現時刻の節点変位ベクトルｄ（ｔ）、節点速度ベクトル、節点加速度ベクトルを計算し、その変位挙動や歪み、応力などの挙動を解析する。

（シミュレーションシステムの構成）
図１３に示すように、シミュレーションシステムは、粘弾性体の形状メッシュデータ、材料データ及び外力データをシステムの記憶部１３６や入力ファイルから入力するデータ入力部１３０と、入力されたデータを用いて質量マトリックスＭ、粘弾性減衰マトリックスＣ_１、剛性行列Ｋ’を作成して記憶部１３６に記憶させるマトリックス演算部１３１と、節点変位ベクトルｄ（ｔ）の履歴を用いて分数階微分値を計算し、形状メッシュデータに対応した一般化等価節点外力ベクトルＰ’（ｔ）を算出して記憶部１３６に記憶させる演算部１３２と、一般化剛性行列Ｋを算出する一般化剛性行列Ｋ演算部１３３と、記憶部１３６に記憶されている一般化剛性行列Ｋと一般化等価節点外力ベクトルＰ’（ｔ）とを用いて剛性方程式を解いて節点変位ベクトルｄ（ｔ）を算出して記憶部１３６に記憶させる変位ベクトル演算部１３４と、入力結果や演算結果をディスプレイ等に出力する出力部１３５と、記憶部１３６とを備えたコンピュータから構成される。データ入力部１３０、マトリックス演算部１３１、一般化等価節点外力ベクトルＰ’（ｔ）演算部１３２、一般化剛性行列Ｋ演算部１３３、変位ベクトル演算部１３４及び出力部１３５と、記憶部１３６とは、入力データや各種算出データを記憶、及び／又は読み出し作業がなされる構成となっている。さらに、節点変位ベクトルｄ（ｔ）、節点変位速度ベクトル、節点変位加速度ベクトルを用いて算出される物性値の挙動を計算する場合には、物性予測値演算部を含んでもよい。なお、上述した一般化等価節点外力ベクトルＰ’（ｔ）演算部１３２において、節点変位ベクトルｄ（ｔ）の分数階微分値Ｄ^ｑｄ（ｔ）を求める際に、第１の実施形態で述べた高速演算装置による分数階微分値計算が実行される。

（解析処理フロー）
以下、本実施形態における有限要素法（ＦＥＭ）による粘弾性材料の外力による形状変形挙動のシミュレーションを図１２のフローチャートを用いて説明する。線形の分数階微分モデルにおいては、節点変位ベクトルｄ（ｔ）の分数階微分値Ｄ ^ｑ ｄ（ｔ）の計算に、第１の実施形態で述べた高速演算装置による分数階微分値計算手法が適用される。

まず、粘弾性体の形状データ、材料データ、等価節点外力ベクトルＰ（ｔ）データ、初期条件、一般化剛性行列Ｋの拘束条件（幾何学的拘束条件）、一般化等価節点ベクトルＰ’（ｔ）の拘束条件（力学的拘束条件）、一般化緩和時間τ等の解析条件を読み込み、記憶部１３６に記憶させる処理等を含む解析初期条件の設定を行う（ステップＳ３１）。

次いで、時間繰り返しインデックスｉをｉ＝０に初期設定する（ステップＳ３２）。

次いで、記憶部１３６に記憶されている粘弾性体の形状データ、材料データ、等価節点外力ベクトルＰ（ｔ）データを読み出して、質量マトリックスＭ、粘弾性減衰マトリックスＣ_１、剛性行列Ｋ’を式（１０５）中のＭの定義式及び式（１０７）中のＣ_１及び式（１０９）中のＫ’の定義式によって算出して記憶部１３６に記憶させる（ステップＳ３３）。

次いで、時間繰り返しインデックスｉをｉ＝ｉ＋１として、時間を更新した後（ステップＳ３４）、節点変位ベクトルｄ（ｔ）分数階微分値Ｄ^ｑｄ（ｔ）を算出して記憶部１３６に記憶させる（ステップＳ３５）。

その後、記憶された節点変位ベクトルｄ（ｔ）の分数階微分値Ｄ^ｑｄ（ｔ）を記憶部１３６から読み出してＳ（ｔ_ｉ−１）を算出し、記憶部１３６に記憶させる（ステップＳ３６）。

次いで、記憶部１３６から等価節点外力ベクトルＰ（ｔ）とＳ（ｔ_ｉ−１）とを読み出して、式（１１０）から一般化等価節点外力ベクトルＰ’（ｔ）を算出して記憶部に記憶させる（ステップＳ３７）。

次いで、記憶部１３６から質量マトリックスＭ、粘弾性減衰マトリックスＣ_１、剛性行列Ｋ及び一般化緩和時間τを読み出して、式（１０９）により一般化剛性行列Ｋ’を算出し、記憶部１３６に記憶させる（ステップＳ３８）。

次いで、一般化剛性行列Ｋ’に幾何学的拘束の処理と、一般化等価節点ベクトルＰ’（ｔ）に力学的拘束の処理を行った後、記憶部１３６に記憶させる（ステップＳ３９）。

次いで、記憶部１３６に記憶されている幾何学的拘束処理された一般化剛性行列Ｋ’と、力学的拘束処理された一般化等価節点外力ベクトルＰ’（ｔ）とを読み出し、式（１０８）によって、時刻ｔの節点変位ベクトルｄ（ｔ_ｉ）、節点速度ベクトル、節点加速度ベクトルを算出し、記憶部１３６に記憶させる（ステップＳ４０）。

その後、計算終了かどうかの判断を行い（ステップＳ４１）、時間繰り返しインデックスｉがｉ≧ｎではない場合（Ｎｏの場合）は、ステップＳ３４へ戻って時間を更新し、ステップＳ３４〜Ｓ４１の同様のループで、各時刻毎の節点変位ベクトルｄ（ｔ_ｉ）、節点速度ベクトル、節点加速度ベクトルを算出し、各時刻の履歴データとして記憶部１３６に記憶させる。この変位ベクトルｄ（ｔ）の各時刻の履歴データは、変位ベクトルの分数階微分値Ｄ^ｑｄ（ｔ_ｉ−１）において、式（３１）の分数階微分定義式におけるｆ（τ）に相当し、この履歴データについて本発明の分数階微分値の計算が適用される。これにより分数階微分値Ｄ^ｑｄ（ｔ）の分数階微分値が精度よく高速で得られ、その結果、シミュレーションの精度を保ちつつ従来に比較して短時間で解析を行うことが可能となる。

ステップＳ４１において、時間繰り返しインデックスｉがｉ≧ｎと判断した場合（Ｙｅｓの場合）は、算出されて記憶部１３６に記憶されている各時刻のデータ、即ち変位ベクトルｄ（ｔ）、変位速度ベクトル、変位加速度ベクトルが、記憶部１３６から読み出され、公知のデータ加工手段によって、グラフィック表示などの形式で出力される（ステップＳ４２）。また、算出して記憶部に記憶されている節点変位ベクトルｄ（ｔ）、節点速度ベクトル、節点加速度ベクトルを用いて、節点歪みベクトルや要素内応力ベクトルを算出することもできる。

なお、図１２においては、分数階微分値Ｄ^ｑｄ（ｔ_ｉ−１）の算出とＰ’（ｔ_ｉ−１）の算出とを分けて実行するように記載されているが、プログラムにおいては、Ｐ’（ｔ_ｉ−１）を算出する計算ステップの中でサブルーチンとして分数階微分値Ｄ^ｑｄ（ｔ_ｉ−１）の算出を実行させてもよい。

本発明の高速演算装置の処理時間効果を検証するために、式（１１２）で与えられる分数階微分モデルの伝達関数と、式（１１３)で与えられる分数階微分モデルの制御システムとからなる閉ループ系ＰＤ^μ制御（α＜μ＜β）における出力の計算時間について、本発明を適用した実施例１と、従来法を適用した比較例１及び２とを比較した。

（実施例１）
前述のＰＤ^μ制御は、出力をｙ（ｔ）、入力をｗ（ｔ）として、式（１１４）の微分方程式でモデル化できる。式（１１４）の微分方程式は、式（１１５）を用いて式（１１４）を（１−β）階微分して式（１１６）へと変形し、式（１１６）を離散化した式（１１７）（添え字「０」はｎ−１ステップ目の値を示す）から現時刻ｔの出力ｙ（ｔ）の値として出力ｙ_ｎを数値計算する。

具体的には、前述の離散化により、式（１１７）右辺の各分数階微分項は式（１１８）及び式（１１９）の分数階積分として計算できるので、式（１１７）中のＡのうち式（１１８）及び式（１１９）中のｊ＝０からｎ−１までの和の部分をＡ_０ｙ_ｎ（実際には出力ｙ_ｎにあたる項は含まれない）、ｊ＝ｎ項目の部分をＡ_ｎｙ_ｎとすると、式（１１７）は式（１２０）となり、式（１２０）の両辺を［１＋ｈＡ_ｎ／２］で除して出力ｙ_ｎを数値計算した。式（１１６）の係数ａ_０及びａ_１は、ａ_０＝ａ_１＝１とし、ゲイン値Ｋ、Ｔ_ｄはそれぞれＫ＝２０．５、Ｔ_ｄ＝３．７３４３とした。また、μ＝１．１５である。

分数階微分値の計算は、式（９４）の条件を満足するパラメータとして、ｍ＝ｍ_１＝２、Ｍ＝４０、Ｋ＝２、ε_０＝１／２０とし、ＲＬ定義に基づく式（４６）によって、区間［０，ｃ_１］の分数階微分値を算出し、区間［ｃ_１，ｔ］の分数階微分値は、式（３２）の第２項のＲＬ微分の値を台形則によって求め、区間［０，ｃ_１］の分数階微分値と区間［ｃ_１，ｔ］の分数階微分値とを合算して全区間［０，ｔ］の分数階微分値を求めるプログラムとした。なお、区間［ｃ_１，ｔ］と［０，ｃ_１］の積分中のベキ式の値は、あらかじめ計算したデータを記憶部に記憶しておいて計算を実行した。

（比較例１）
比較例１として、実施例１においてＲＬ定義による分数階微分値の計算を「台形則」で計算を行った。「台形則」とは、台形公式はニュートン・コーツの公式と呼ばれる数値積分の公式のひとつであって、被積分関数を一次関数で近似し、台形の面積の公式に帰着させて積分の近似値を求める方法であり、式（１２１）の要領で計算するものである。このＲＬ微分中のベキ式の値は、あらかじめ計算したデータを記憶部に記憶しておき、これを読み出して計算を実行した。

（比較例２）
比較例２として、実施例１において従来のＧＬ定義による分数階微分値の計算を行った。「ＧＬ定義」による計算法は、式（１２２）を用いて計算する。ＧＬ定義による計算において、式（１２２）の係数である式（１２３）は、あらかじめ計算して記憶部に記憶させておき、これを読み出して式（１２２）の計算を実行している。

実施例１並びに比較例１及び２について、計算時間（ＣＰＵ時間）を計算ステップ数に対してプロットした結果を図１４に示す。図１４より、本発明の高速計算装置では、計算の数千ステップ以上における台形則、ＧＬ定義の場合に比べてＣＰＵ時間が低く抑えられていることが分かる。このことから、本発明の高速演算装置を用いた制御は、従来法よりも極めて短時間で出力計算できるので、追従性の高い制御が実現できることが分かる。計算の精度は、一次の精度でステップ数が等しければ、本発明の実施例１、比較例１である台形則、及び比較例２であるＧＬ法の精度はほぼ同等である。

本発明の高速演算装置の処理時間効果を検証するために、シミュレーションシステムについて、式（１２４）の分数階微分を含む粘弾性体の１自由度系の振動子モデルの分数階微分値計算の処理時間について、本発明を適用した実施例２と、従来法を適用した比較例３及び４とを比較した。

ただし、式（１２４）において、ｍ_ｓは質量、ｘ（ｔ）は変位、ｔは時間、ｃは粘弾性係数、ｑは分数階数（０より大きい実数）、ｋ_ｓはバネ定数、ｆ（ｔ）は質量に与えられる外力を表す。

（実施例２）
まず、時刻ｔ＞０において、式（１２４）にｆ（ｔ）＝ｓｉｎ（ｔ）を入力した場合の応答を計算する。その他の式（１０９）のパラメータは、ｍ＝ｃ＝ｋ＝１，ｑ＝０．５とおいた。本発明の高速計算装置を使用して、式（１２４）の分数階微分値の計算を行った。この実施例２の分数階微分値の計算は、式（９４）の条件を満足するパラメータとして、ｍ＝ｍ_１＝２、Ｍ＝４０、Ｋ＝２、ε_０＝１／２０とし、ＲＬ定義に基づく式（４６）によって、区間［０，ｃ_１］の分数階微分値を算出し、区間［ｃ_１，ｔ］の分数階微分値は、式（３２）の第２項のＲＬ微分の値を台形則によって求め、区間［０，ｃ_１］の分数階微分値と区間［ｃ_１，ｔ］の分数階微分値とを合算して全区間［０，ｔ］の分数階微分値を求めるプログラムとした。なお、区間［ｃ_１，ｔ］と［０，ｃ_１］の積分中のベキ式の値は、あらかじめ計算したデータを記憶部に記憶しておいて計算を実行した。

（比較例３）
実施例２において、分数階微分値の計算を比較例１と同じ台形則で計算した。

（比較例４）
実施例２において、分数階微分値の計算を比較例２と同じＧＬ法で計算した。

実施例２並びに比較例３及び４について、計算時間（ＣＰＵ時間）を計算ステップ数に対してプロットした結果を図１５に示す。図１５より、本発明の高速計算方法では、計算の数千ステップ以上における台形則、ＧＬ定義の場合に比べてＣＰＵ時間が低く抑えられていることが分かる。各ステップ数までに必要なメモリ数は、台形則ではステップ数に等しい。ＧＬ法ではｘ（ｔ）のための記憶部と係数のための記憶部が必要で、ステップ数の２倍である。従って、台形法とＧＬ法はＮ^２アルゴリズムの結果となっている。一方、本発明によれば、必要なメモリ数は、Ｍ＝４０、ｍ＝ｍ_１＝２、Ｋ＝２の場合、ステップ数＝４０×２^ｎに対して４３＋３×２２ｎであり、ＮｌｏｇＮアルゴリズムである。計算の精度は、一次の精度でステップ数が等しければ、本発明の実施例２、比較例３である台形則、及び比較例４であるＧＬ法の精度はほぼ同等である。

このように本発明により高速計算アルゴリズムを使用することによって、これまで実現が難しいとされてきた分数階微分方程式の数値計算が飛躍的に短時間で計算でき、分数階微分要素を含む現象の解析などを実用的な時間で行うことができる。

以上述べた実施形態は全て本発明を例示的に示すものであって限定的に示すものではなく、本発明は他の種々の変形態様及び変更態様で実施することができる。従って本発明の範囲は特許請求の範囲及びその均等範囲によってのみ規定されるものである。

本発明の分数階微分値の計算を行う高速演算装置は、分数階微分要素を含む様々な現象の経時的な挙動や特性解析（シミュレーション）や信号制御、画像処理、経済推移予測などに幅広く適用でき、この装置を用いれば、精度よく短時間で処理できるので、商品開発の効率アップや、生産性、製品の品質向上、さらに経済変動の予測情報の提供などにも利用することができる。

１０ コンピュータ
１０ａ バス
１０ｂ ＣＰＵ
１０ｃ ＲＯＭ
１０ｄ ＲＡＭ
１０ｅ ＨＤＤ
１０ｆ 音声処理部
１０ｇ 画像処理部
１０ｈ ディスク駆動装置
１０ｉ 入出力インタフェース
１０ｊ スピーカ
１０ｋ ディスプレイ
１０ｌ キーボード
１０ｍ マウス
１０ｎ プリンタ
１１ 外部装置
３０ 入力部
３１ 二次パラメータ演算部
３２ 第１の分数階微分値演算部
３３ 重み付き積分演算部
３４ 第２の分数階微分値演算部
３５ 全分数階微分値算出部
３６、１３５ 出力部
３７、１０５、１３６ 記憶部
１００ 制御パラメータ入力部
１０１ Ｗ（ｓ）及びＹ（ｓ）入力部
１０２ 差分Ｅ(ｓ)演算部
１０３ 制御システム出力値Ｕ（ｓ）演算部
１０４ Ｕ（ｓ）出力部
１０６ コンピュータ
１０７ 制御対象
１３０ データ入力部
１３１ マトリックス演算部
１３２ 一般化等価節点外力ベクトルＰ’（ｔ）演算部
１３３ 一般剛性行列Ｋ演算部
１３４ 変位ベクトル演算部

Claims (18)

1. 分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求める高速演算装置であって、
   パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
   前記演算部は、
   前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ_１進む毎に式（１）を満たしながら更新される分割点ｃ_１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する第１の分割設定手段と、
   

   

   （ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
   前記変数ｔが幅ｈ更新される度に、前記区間［ｃ_１，ｔ］又は前記区間［ｔ，ｃ_１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
   前記区間［ａ，ｃ_１］又は前記区間［ｃ_１，ａ］において、前記変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、前記分割点ｃ_１を上端とする幅Δ_１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
   前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ_１］又は該区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
   前記区間［ｃ_１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ_１］の分数階微積分値Ａと前記区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
   を備えていることを特徴とする高速演算装置。
2. 前記変数ｔに関する積分区間が［ａ，ｔ］であって、
   前記第１の分割設定手段は、式（２）の条件を満たす設定値Ｍ及びｍにおいて、ｔ−（Ｍ＋ｍ）ｈ＝ｃ_１となった時点を分割点ｃ_１として設定するものであり、
   前記積分演算手段は、前記区間［ａ，ｃ_１］の前記変数ｔの更新に依存しない重み付き積分値の算出を、前記分割点ｃ_１が設定された時点からｃ_１＝ａとして開始し、以降、前記変数ｔがｍｈ（＝Δ_１）進みｃ_１＝ｔ−（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｃ_１，ｔ−Ｍｈ］の幅ｈ刻みのｆ（τ）データを用い、ｃ_１＝ｔ−Ｍｈへの更新を伴って、前記更新されたｃ_１を上端とする幅Δ_１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して前記記憶部に記憶させる第１の積分演算部を備えている、
   

   

   （ただし、ε_Ｔは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ_１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
   ことを特徴とする請求項１に記載の高速演算装置。
3. 前記変数ｔに関する積分区間が［ｔ，ａ］であって、
   前記第１の分割設定手段は、式（３）の条件を満たす設定値Ｍ及びｍにおいて、ｔ＋（Ｍ＋ｍ）ｈ＝ｃ_１となった時点を分割点ｃ_１として設定するものであり、
   前記積分演算手段は、前記区間［ｃ_１，ａ］の前記変数ｔの更新に依存しない重み付き積分値の算出を、前記分割点ｃ_１が設定された時点からｃ_１＝ａとして開始し、以降、前記変数ｔがｍｈ（＝Δ_１）進みｃ_１＝ｔ＋（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｔ＋Ｍｈ，ｃ_１］の幅ｈ刻みのｆ（τ）データを用い、ｃ_１＝ｔ＋Ｍｈへの更新を伴って、前記更新されたｃ_１を上端とする幅Δ_１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して前記記憶部に記憶させる第１の積分演算部を備えている、
   

   

   （ただし、ε_Ｔは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ_１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
   ことを特徴とする請求項１に記載の高速演算装置。
4. 前記第１の分割設定手段によって、式（４）の条件を満たす設定値Ｍ及びｍにおいて、前記変数ｔのｈ刻みの更新に伴い前記分割点ｃ_１がｃ_１＝ｔ−Ｍｈとして更新され、さらに、ｃ_１＝ａを初期値として、区間［ａ，ｃ_１］内のｃ_１を上端とする幅Δ’が前記変数ｔのｈ刻みの更新によってｈからｍｈまで変化するものであって、
   

   

   （ただし、ε_Ｔは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ_１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
   前記積分演算手段は、前記幅Δ’がｍｈ未満の場合には、ｃ_１を上端とする前記幅Δ’の小区間は幅ｈ刻みのｆ（τ）データを用い、前記変数ｔの更新に依存する積分値の計算を前記変数ｔの更新毎に繰り返し実行して前記記憶部に記憶させ、前記幅Δ’がｍｈとなった場合には、前記小区間の幅ｈ刻みのｆ（τ）データを用い、分割点ｃ_１を上端とする幅Δ_１（＝ｍｈ）の小区間の前記変数ｔの更新に依存しない重み付け積分値を算出する重み付き計算を前記変数ｔの更新毎に繰り返し実行して前記記憶部に記憶させる第１の積分演算部を備えており、
   前記第２の分数階微積分値演算手段は、前記幅Δ’がｍｈとなった場合に分割点ｃ_１に対応する重み付き積分値を前記記憶部から読み出して、前記忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して算出する分数階微積分値と、前記変数ｔの更新に依存する前記幅Δ’がｍｈ未満の場合の小区間の積分値から算出する分数階微積分値とを合算して算出する第１の分数階微積分値演算部を備えている
   ことを特徴とする請求項１に記載の高速演算装置。
5. 前記積分演算手段は、前記変数ｔを更新しながら前記区間［ａ，ｃ_１］の重み付き積分を算出する際に、前記幅Δ_１の小区間の個数Ｊ_ｓ−１（ｓ＝２）が式（５）で表されるＬ_ｓと等しくなった時点で、ｃ_２＝ｃ_１−（Ｌ_ｓ−ｍ_１）Δ_１である新たな分割点ｃ_２を設定して前記区間［ａ，ｃ_１］を中区間［ａ，ｃ_２］と中区間［ｃ_２，ｃ_１］とに分割する第２の分割設定手段を具備し、
   前記中区間［ａ，ｃ_２］においては、前記ｃ_２の設定された時点並びに以降、前記変数ｔがｍ_１Δ_１進む度に、前記分割点ｃ_２から直近の前記分割点ｃ_１方向に連続するｍ_１個の幅Δ_１の小区間を纏めて幅Δ_２の小区間とし、前記記憶部に記憶されている前記ｍ_１個の各小区間の重み付き積分値を重み付き積分値として算出して前記記憶部に記憶させる第２の積分演算部を備えており、
   

   

   前記第２の分数階微積分値演算手段は、前記区間［ａ，ｃ_１］の分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ_２，ｃ_１］と、前記中区間［ａ，ｃ_２］の各小区間の重み付け積分値、をそれぞれ読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
   ことを特徴とする請求項２に記載の高速演算装置。
6. 前記積分演算手段は、前記変数ｔを更新しながら前記区間［ｃ_１，ａ］の重み付き積分を算出する際に、前記幅Δ_１の小区間の個数Ｊ_ｓ−１（ｓ＝２）が式（６）で表されるＬ_ｓと等しくなった時点で、ｃ_２＝ｃ_１＋（Ｌ_ｓ−ｍ_１）Δ_１である新たな分割点ｃ_２を設定して前記区間［ｃ_１，ａ］を中区間［ｃ_１，ａ］と中区間［ｃ_１，ｃ_２］とに分割する第２の分割設定手段を具備し、
   前記中区間［ｃ_２，ａ］においては、前記ｃ_２の設定された時点並びに以降、前記変数ｔがｍ_１Δ_１進む度に、前記分割点ｃ_２から直近の前記分割点ｃ_１方向に連続するｍ_１個の幅Δ_１の小区間を纏めて幅Δ_２の小区間とし、前記記憶部に記憶されている前記ｍ_１個の各小区間の重み付き積分値を重み付き積分値として算出して前記記憶部に記憶させる第２の積分演算部を備えており、
   

   

   前記第２の分数階微積分値演算手段は、前記区間［ｃ_１，ａ］の分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ_１，ｃ_２］と、前記中区間［ｃ_２，ａ］の各小区間の重み付け積分値、をそれぞれ読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
   ことを特徴とする請求項３に記載の高速演算装置。
7. 前記積分演算手段は、前記変数ｔを更新しながら前記区間［ａ，ｃ_２］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃ_ｓ，ｃ_ｓ−１］の幅Δ_ｓ−１の小区間の個数Ｊ_ｓ−１（ｓ≧３）が式（７）で表されるＬ_ｓと等しくなる度に、ｃ_ｓ＝ｃ_ｓ−１−（Ｌ_ｓ−ｍ_１）Δ_ｓ−１である新たな分割点ｃ_ｓによってｃ_ｓ＋１＝ａを初期値とする前記中区間［ｃ_ｓ＋１，ｃ_ｓ］と［ｃ_ｓ，ｃ_ｓ−１］とを設定する第３の分割手段を具備し、前記ｃ_ｓが設定された時点及びそれ以降、前記変数ｔがｍ_１Δ_ｓ−１進む度に、前記分割点ｃ_ｓから直近の前記分割点ｃ_ｓ−１方向に連続するｍ_１個のΔ_ｓ−１の小区間を纏めて幅Δ_ｓの小区間とし、前記記憶部に記憶されている前記ｍ_１個のΔ_ｓ−１の各小区間の重み付き積分値を用いて前記Δ_ｓの小区間の重み付き積分値として算出して前記記憶部に記憶させ、さらに、以降前記変数ｔの更新によってΔ_ｓの小区間の個数がＬ_ｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、
   

   

   前記第２の分数階微積分値演算手段は、前記区間［ａ，ｃ_１］又は分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ_２，ｃ_１］及び前記中区間［ｃ_３，ｃ_２］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする前記中区間［ｃ_ｓ＋１，ｃ_ｓ］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する前記忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
   ことを特徴とする請求項４に記載の高速演算装置。
8. 前記積分演算手段は、前記変数ｔを更新しながら前記区間［ｃ_２，ａ］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃ_ｓ−１，ｃ_ｓ］の幅Δ_ｓ−１の小区間の個数Ｊ_ｓ−１（ｓ≧３）が式（８）で表されるＬ_ｓと等しくなる度に、ｃ_ｓ＝ｃ_ｓ−１＋（Ｌ_ｓ−ｍ_１）Δ_ｓ−１である新たな分割点ｃ_ｓによってｃ_ｓ＋１＝ａを初期値とする前記中区間［ｃ_ｓ，ｃ_ｓ＋１］と［ｃ_ｓ−１，ｃ_ｓ］とを設定する第３の分割手段を具備し、前記ｃ_ｓが設定された時点及びそれ以降、前記変数ｔがｍ_１Δ_ｓ−１進む度に、前記分割点ｃ_ｓから直近の前記分割点ｃ_ｓ−１方向に連続するｍ_１個のΔ_ｓ−１の小区間を纏めて幅Δ_ｓの小区間とし、前記記憶部に記憶されている前記ｍ_１個のΔ_ｓ−１の各小区間の重み付き積分値を用いて前記Δ_ｓの小区間の重み付き積分値として算出して前記記憶部に記憶させ、さらに、以降前記変数ｔの更新によってΔ_ｓの小区間の個数がＬ_ｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、
   

   

   前記第２の分数階微積分値演算手段は、前記区間［ｃ_１，ａ］の分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ_１，ｃ_２］及び前記中区間［ｃ_２，ｃ_３］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする前記中区間［ｃ_ｓ，ｃ_ｓ＋１］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する前記忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
   ことを特徴とする請求項６に記載の高速演算装置。
9. 前記ｍが２以上であることを特徴とする請求項２から８のいずれか１項に記載の高速演算装置。
10. 前記ｍ_１が２であることを特徴とする請求項５から８のいずれか１項に記載の高速演算装置。
11. 前記忘却関数が式（９）であり、
    前記積分演算手段は、前記式（９）の忘却関数又は前記式（９）の分子をあらかじめ計算し、得られたデータを前記記憶部に記憶させる第３の積分演算部を備えており、
    前記第２の分数階微分値演算手段を実行する
    

    

    （ただし、ｂは、個々の小区間を識別するための任意の点、αは分数階微積分の積分定義で規定される積分階数、ｇは分数階微積分の積分定義で規定される関数、ｋは展開次数番号である）
    ことを特徴とする請求項１、２、４、５及び７のいずれか１項に記載の高速演算装置。
12. 前記第３の積分演算部は、前記式（９）の分子をあらかじめ計算し、得られたデータを前記記憶部に記憶させるものであり、
    前記積分演算手段は、さらに、前記重み付き積分にあらかじめ前記式（９）中の係数１／（ｋ！Γ（α−ｋ））を乗じた重み付き積分値として算出して前記記憶部に記憶させる第４の積分演算部を備えており、
    前記第２の分数階微積分値演算手段は、前記記憶部に記憶されている小区間の前記係数を乗じた重み付け積分値を読み出し、全ての中区間の各小区間に対応する式（９）の前記忘却関数の分子を乗じた値を前記小区間毎に算出し、全小区間について合算することによって区間［ａ，ｃ_１］の分数階微積分値を算出する第３の分数階微積分値演算部を備えている
    ことを特徴とする請求項１１に記載の高速演算装置。
13. 前記積分式で定義される分数階微積分要素の分数階微積分式が、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義、Ｃａｐｕｔｏの定義及び分数階積分の定義式のいずれかであることを特徴とする請求項１から１２のいずれか１項に記載の高速演算装置。
14. 前記区間［ａ，ｃ_１］の分数階微積分定義がＲｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義であり、区間［ｃ_１，ｔ］の分数階微積分定義がＣａｐｕｔｏの定義であることを特徴とする請求項１、２、４、５及び７のいずれか１項に記載の高速演算装置。
15. 分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微分要素の前記変数ｔにおける分数階微分値を求めるコンピュータの高速演算プログラムであって、
    前記コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
    前記コンピュータを、
    前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ_１進む毎に式（１０）を満たしながら更新される分割点ｃ_１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、
    

    

    （ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
    前記変数ｔが幅ｈ更新される度に、前記区間［ｃ_１，ｔ］又は前記区間［ｔ，ｃ_１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
    前記区間［ａ，ｃ_１］又は前記区間［ｃ_１，ａ］において、前記変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、前記分割点ｃ_１を上端とする幅Δ_１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
    前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ_１］又は該区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
    前記区間［ｃ_１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ_１］の分数階微積分値Ａと前記区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
    して機能させる高速演算プログラム。
16. 分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求めるコンピュータの高速演算プログラムであって、
    前記コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
    前記コンピュータを、
    前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ_１進む毎に式（１１）を満たしながら更新される分割点ｃ_１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、
    

    

    （ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
    前記変数ｔが幅ｈ更新される度に、前記区間［ｃ_１，ｔ］又は前記区間［ｔ，ｃ_１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
    前記区間［ａ，ｃ_１］又は前記区間［ｃ_１，ａ］において、前記変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、前記分割点ｃ_１を上端とする幅Δ_１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
    前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ_１］又は該区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
    前記区間［ｃ_１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ_１］の分数階微積分値Ａと前記区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
    して機能させるための高速演算プログラムを記録したコンピュータ読み取り可能な記録媒体。
17. 設定データである入力データと制御対象のシステム出力とを入力するデータ入力部と、該データ入力部から得られる入力データとシステム出力との差分を計算する差分演算部と、該差分演算部から得られる差分と制御器の分数階微分要素を有する伝達関数とから制御器出力値を計算する制御器出力演算部と、該制御器出力演算部から得られる制御器出力値を制御対象に出力する制御器出力用出力部とを備えており、前記制御器出力用出力部から得られる制御器出力値と前記制御対象の伝達関数とからシステム出力を出力して前記制御対象の被制御機器を作動させる制御を、前記システム出力と次の設定データとから繰り返し実行するように構成した機器制御システムであって、
    前記制御器出力演算部は、前記制御器の前記伝達関数における分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
    前記高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
    前記高速演算装置の前記演算部は、
    前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ_１進む毎に式（１２）を満たしながら更新される分割点ｃ_１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、
    

    

    （ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
    前記変数ｔが幅ｈ更新される度に、前記区間［ｃ_１，ｔ］又は前記区間［ｔ，ｃ_１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
    前記区間［ａ，ｃ_１］又は前記区間［ｃ_１，ａ］において、前記変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、前記分割点ｃ_１を上端とする幅Δ_１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
    前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ_１］又は該区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
    前記区間［ｃ_１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ_１］の分数階微積分値Ａと前記区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
    を備えていることを特徴とする機器制御システム。
18. 物理現象が時間を変数とした分数階微積分要素を含む分数階微積分式を用い、該分数階微積分式が積分式で表現されるようにモデル化し、前記分数階微積分要素から分数階微積分値を計算する第１の演算手段と、前記分数階微積分値を用いて運動方程式の解を計算する第２の演算手段と、物理現象の測定データと前記運動方程式の経時的な解とを用いて前記物理現象の経時的な特性値を計算する第３の演算手段とを備え、物理現象の挙動を解析するシミュレーションシステムであって、
    前記第１の演算手段は、前記分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
    前記高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
    前記高速演算装置の前記演算部は、
    前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ_１進む毎に式（１３）を満たしながら更新される分割点ｃ_１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ_１］と区間［ｃ_１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ_１］と区間［ｃ_１，ａ］とに分割する分割設定手段と、
    

    

    （ただし、Δ_１は分割点ｃ_１を上端とする変数の小区間の幅（Δ_１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
    前記変数ｔが幅ｈ更新される度に、前記区間［ｃ_１，ｔ］又は前記区間［ｔ，ｃ_１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
    前記区間［ａ，ｃ_１］又は前記区間［ｃ_１，ａ］において、前記変数ｔがΔ_１更新される度に、分割点ｃ_１を端点とする幅Δ_１の小区間内のｆ（τ）を用いて、前記分割点ｃ_１を上端とする幅Δ_１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
    前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ_１］又は該区間［ｃ_１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
    前記区間［ｃ_１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ_１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ_１］の分数階微積分値Ａと前記区間［ｃ_１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
    を備えていることを特徴とするシミュレーションシステム。

Images