(57)【要約】分数階微積分演算を行う高速演算装置であって、関数ｆ（τ）の積分区間［ａ，ｔ］を、区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割し、区間［ａ，ｃ１］においては、計算してメモリに記憶しておいたｔに依存しない重み付き積分値を用いて微積分演算を行い、区間［ｃ１，ｔ］について計算した微積分演算の結果と合算することにより、微積分演算結果を得る。ｔが更新されるのにしたがって、新たな分割点ｃ２を設定し、区間［ｃ１，ｔ］における計算結果に基づいて区間［ａ，ｃ２］の重み付き積分値を算出し、メモリに記憶する。
2013024810.tif 000002
【特許請求の範囲】
【請求項１】
分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求める高速演算装置であって、
パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
前記演算部は、
前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ１進む毎に式（１）を満たしながら更新される分割点ｃ１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する第１の分割設定手段と、
【数１】
2013024810.tif 000127
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
前記変数ｔが幅ｈ更新される度に、前記区間［ｃ１，ｔ］又は前記区間［ｔ，ｃ１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
前記区間［ａ，ｃ１］又は前記区間［ｃ１，ａ］において、前記変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、前記分割点ｃ１を上端とする幅Δ１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ１］又は該区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
前記区間［ｃ１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ１］の分数階微積分値Ａと前記区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えていることを特徴とする高速演算装置。
【請求項２】
前記変数ｔに関する積分区間が［ａ，ｔ］であって、
前記第１の分割設定手段は、式（２）の条件を満たす設定値Ｍ及びｍにおいて、ｔ−（Ｍ＋ｍ）ｈ＝ｃ１となった時点を分割点ｃ１として設定するものであり、
前記積分演算手段は、前記区間［ａ，ｃ１］の前記変数ｔの更新に依存しない重み付き積分値の算出を、前記分割点ｃ１が設定された時点からｃ１＝ａとして開始し、以降、前記変数ｔがｍｈ（＝Δ１）進みｃ１＝ｔ−（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｃ１，ｔ−Ｍｈ］の幅ｈ刻みのｆ（τ）データを用い、ｃ１＝ｔ−Ｍｈへの更新を伴って、前記更新されたｃ１を上端とする幅Δ１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して前記記憶部に記憶させる第１の積分演算部を備えている、
【数２】
2013024810.tif 000128
（ただし、εＴは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
ことを特徴とする請求項１に記載の高速演算装置。
【請求項３】
前記変数ｔに関する積分区間が［ｔ，ａ］であって、
前記第１の分割設定手段は、式（３）の条件を満たす設定値Ｍ及びｍにおいて、ｔ＋（Ｍ＋ｍ）ｈ＝ｃ１となった時点を分割点ｃ１として設定するものであり、
前記積分演算手段は、前記区間［ｃ１，ａ］の前記変数ｔの更新に依存しない重み付き積分値の算出を、前記分割点ｃ１が設定された時点からｃ１＝ａとして開始し、以降、前記変数ｔがｍｈ（＝Δ１）進みｃ１＝ｔ＋（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｔ＋Ｍｈ，ｃ１］の幅ｈ刻みのｆ（τ）データを用い、ｃ１＝ｔ＋Ｍｈへの更新を伴って、前記更新されたｃ１を上端とする幅Δ１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して前記記憶部に記憶させる第１の積分演算部を備えている、
【数３】
2013024810.tif 000129
（ただし、εＴは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
ことを特徴とする請求項１に記載の高速演算装置。
【請求項４】
前記第１の分割設定手段によって、式（４）の条件を満たす設定値Ｍ及びｍにおいて、前記変数ｔのｈ刻みの更新に伴い前記分割点ｃ１がｃ１＝ｔ−Ｍｈとして更新され、さらに、ｃ１＝ａを初期値として、区間［ａ，ｃ１］内のｃ１を上端とする幅Δ’が前記変数ｔのｈ刻みの更新によってｈからｍｈまで変化するものであって、
【数４】
2013024810.tif 000130
（ただし、εＴは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
前記積分演算手段は、前記幅Δ’がｍｈ未満の場合には、ｃ１を上端とする前記幅Δ’の小区間は幅ｈ刻みのｆ（τ）データを用い、前記変数ｔの更新に依存する積分値の計算を前記変数ｔの更新毎に繰り返し実行して前記記憶部に記憶させ、前記幅Δ’がｍｈとなった場合には、前記小区間の幅ｈ刻みのｆ（τ）データを用い、分割点ｃ１を上端とする幅Δ１（＝ｍｈ）の小区間の前記変数ｔの更新に依存しない重み付け積分値を算出する重み付き計算を前記変数ｔの更新毎に繰り返し実行して前記記憶部に記憶させる第１の積分演算部を備えており、
前記第２の分数階微積分値演算手段は、前記幅Δ’がｍｈとなった場合に分割点ｃ１に対応する重み付き積分値を前記記憶部から読み出して、前記忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して算出する分数階微積分値と、前記変数ｔの更新に依存する前記幅Δ’がｍｈ未満の場合の小区間の積分値から算出する分数階微積分値とを合算して算出する第１の分数階微積分値演算部を備えている
ことを特徴とする請求項１に記載の高速演算装置。
【請求項５】
前記積分演算手段は、前記変数ｔを更新しながら前記区間［ａ，ｃ１］の重み付き積分を算出する際に、前記幅Δ１の小区間の個数Ｊｓ−１（ｓ＝２）が式（５）で表されるＬｓと等しくなった時点で、ｃ２＝ｃ１−（Ｌｓ−ｍ１）Δ１である新たな分割点ｃ２を設定して前記区間［ａ，ｃ１］を中区間［ａ，ｃ２］と中区間［ｃ２，ｃ１］とに分割する第２の分割設定手段を具備し、
前記中区間［ａ，ｃ２］においては、前記ｃ２の設定された時点並びに以降、前記変数ｔがｍ１Δ１進む度に、前記分割点ｃ２から直近の前記分割点ｃ１方向に連続するｍ１個の幅Δ１の小区間を纏めて幅Δ２の小区間とし、前記記憶部に記憶されている前記ｍ１個の各小区間の重み付き積分値を重み付き積分値として算出して前記記憶部に記憶させる第２の積分演算部を備えており、
【数５】
2013024810.tif 000131
前記第２の分数階微積分値演算手段は、前記区間［ａ，ｃ１］の分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ２，ｃ１］と、前記中区間［ａ，ｃ２］の各小区間の重み付け積分値、をそれぞれ読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
ことを特徴とする請求項２に記載の高速演算装置。
【請求項６】
前記積分演算手段は、前記変数ｔを更新しながら前記区間［ｃ１，ａ］の重み付き積分を算出する際に、前記幅Δ１の小区間の個数Ｊｓ−１（ｓ＝２）が式（６）で表されるＬｓと等しくなった時点で、ｃ２＝ｃ１＋（Ｌｓ−ｍ１）Δ１である新たな分割点ｃ２を設定して前記区間［ｃ１，ａ］を中区間［ｃ１，ａ］と中区間［ｃ１，ｃ２］とに分割する第２の分割設定手段を具備し、
前記中区間［ｃ２，ａ］においては、前記ｃ２の設定された時点並びに以降、前記変数ｔがｍ１Δ１進む度に、前記分割点ｃ２から直近の前記分割点ｃ１方向に連続するｍ１個の幅Δ１の小区間を纏めて幅Δ２の小区間とし、前記記憶部に記憶されている前記ｍ１個の各小区間の重み付き積分値を重み付き積分値として算出して前記記憶部に記憶させる第２の積分演算部を備えており、
【数６】
2013024810.tif 000132
前記第２の分数階微積分値演算手段は、前記区間［ｃ１，ａ］の分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ１，ｃ２］と、前記中区間［ｃ２，ａ］の各小区間の重み付け積分値、をそれぞれ読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
ことを特徴とする請求項３に記載の高速演算装置。
【請求項７】
前記積分演算手段は、前記変数ｔを更新しながら前記区間［ａ，ｃ２］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃｓ，ｃｓ−１］の幅Δｓ−１の小区間の個数Ｊｓ−１（ｓ≧３）が式（７）で表されるＬｓと等しくなる度に、ｃｓ＝ｃｓ−１−（Ｌｓ−ｍ１）Δｓ−１である新たな分割点ｃｓによってｃｓ＋１＝ａを初期値とする前記中区間［ｃｓ＋１，ｃｓ］と［ｃｓ，ｃｓ−１］とを設定する第３の分割手段を具備し、前記ｃｓが設定された時点及びそれ以降、前記変数ｔがｍ１Δｓ−１進む度に、前記分割点ｃｓから直近の前記分割点ｃｓ−１方向に連続するｍ１個のΔｓ−１の小区間を纏めて幅Δｓの小区間とし、前記記憶部に記憶されている前記ｍ１個のΔｓ−１の各小区間の重み付き積分値を用いて前記Δｓの小区間の重み付き積分値として算出して前記記憶部に記憶させ、さらに、以降前記変数ｔの更新によってΔｓの小区間の個数がＬｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、
【数７】
2013024810.tif 000133
前記第２の分数階微積分値演算手段は、前記区間［ａ，ｃ１］又は分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ２，ｃ１］及び前記中区間［ｃ３，ｃ２］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする前記中区間［ｃｓ＋１，ｃｓ］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する前記忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
ことを特徴とする請求項４に記載の高速演算装置。
【請求項８】
前記積分演算手段は、前記変数ｔを更新しながら前記区間［ｃ２，ａ］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃｓ−１，ｃｓ］の幅Δｓ−１の小区間の個数Ｊｓ−１（ｓ≧３）が式（８）で表されるＬｓと等しくなる度に、ｃｓ＝ｃｓ−１＋（Ｌｓ−ｍ１）Δｓ−１である新たな分割点ｃｓによってｃｓ＋１＝ａを初期値とする前記中区間［ｃｓ，ｃｓ＋１］と［ｃｓ−１，ｃｓ］とを設定する第３の分割手段を具備し、前記ｃｓが設定された時点及びそれ以降、前記変数ｔがｍ１Δｓ−１進む度に、前記分割点ｃｓから直近の前記分割点ｃｓ−１方向に連続するｍ１個のΔｓ−１の小区間を纏めて幅Δｓの小区間とし、前記記憶部に記憶されている前記ｍ１個のΔｓ−１の各小区間の重み付き積分値を用いて前記Δｓの小区間の重み付き積分値として算出して前記記憶部に記憶させ、さらに、以降前記変数ｔの更新によってΔｓの小区間の個数がＬｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、
【数８】
2013024810.tif 000134
前記第２の分数階微積分値演算手段は、前記区間［ｃ１，ａ］の分数階微積分値Ｂを、前記記憶部に記憶されている前記中区間［ｃ１，ｃ２］及び前記中区間［ｃ２，ｃ３］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする前記中区間［ｃｓ，ｃｓ＋１］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する前記忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えている
ことを特徴とする請求項６に記載の高速演算装置。
【請求項９】
前記ｍが２以上であることを特徴とする請求項２から８のいずれか１項に記載の高速演算装置。
【請求項１０】
前記ｍ１が２であることを特徴とする請求項５から８のいずれか１項に記載の高速演算装置。
【請求項１１】
前記忘却関数が式（９）であり、
前記積分演算手段は、前記式（９）の忘却関数又は前記式（９）の分子をあらかじめ計算し、得られたデータを前記記憶部に記憶させる第３の積分演算部を備えており、
前記第２の分数階微分値演算手段を実行する
【数９】
2013024810.tif 000135
（ただし、ｂは、個々の小区間を識別するための任意の点、αは分数階微積分の積分定義で規定される積分階数、ｇは分数階微積分の積分定義で規定される関数、ｋは展開次数番号である）
ことを特徴とする請求項１、２、４、５及び７のいずれか１項に記載の高速演算装置。
【請求項１２】
前記第３の積分演算部は、前記式（９）の分子をあらかじめ計算し、得られたデータを前記記憶部に記憶させるものであり、
前記積分演算手段は、さらに、前記重み付き積分にあらかじめ前記式（９）中の係数１／（ｋ！Γ（α−ｋ））を乗じた重み付き積分値として算出して前記記憶部に記憶させる第４の積分演算部を備えており、
前記第２の分数階微積分値演算手段は、前記記憶部に記憶されている小区間の前記係数を乗じた重み付け積分値を読み出し、全ての中区間の各小区間に対応する式（９）の前記忘却関数の分子を乗じた値を前記小区間毎に算出し、全小区間について合算することによって区間［ａ，ｃ１］の分数階微積分値を算出する第３の分数階微積分値演算部を備えている
ことを特徴とする請求項１１に記載の高速演算装置。
【請求項１３】
前記積分式で定義される分数階微積分要素の分数階微積分式が、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義、Ｃａｐｕｔｏの定義及び分数階積分の定義式のいずれかであることを特徴とする請求項１から１２のいずれか１項に記載の高速演算装置。
【請求項１４】
前記区間［ａ，ｃ１］の分数階微積分定義がＲｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義であり、区間［ｃ１，ｔ］の分数階微積分定義がＣａｐｕｔｏの定義であることを特徴とする請求項１、２、４、５及び７のいずれか１項に記載の高速演算装置。
【請求項１５】
分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微分要素の前記変数ｔにおける分数階微分値を求めるコンピュータの高速演算プログラムであって、
前記コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
前記コンピュータを、
前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ１進む毎に式（１０）を満たしながら更新される分割点ｃ１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１０】
2013024810.tif 000136
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
前記変数ｔが幅ｈ更新される度に、前記区間［ｃ１，ｔ］又は前記区間［ｔ，ｃ１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
前記区間［ａ，ｃ１］又は前記区間［ｃ１，ａ］において、前記変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、前記分割点ｃ１を上端とする幅Δ１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ１］又は該区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
前記区間［ｃ１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ１］の分数階微積分値Ａと前記区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
して機能させる高速演算プログラム。
【請求項１６】
分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求めるコンピュータの高速演算プログラムであって、
前記コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
前記コンピュータを、
前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ１進む毎に式（１１）を満たしながら更新される分割点ｃ１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１１】
2013024810.tif 000137
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
前記変数ｔが幅ｈ更新される度に、前記区間［ｃ１，ｔ］又は前記区間［ｔ，ｃ１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
前記区間［ａ，ｃ１］又は前記区間［ｃ１，ａ］において、前記変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、前記分割点ｃ１を上端とする幅Δ１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ１］又は該区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
前記区間［ｃ１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ１］の分数階微積分値Ａと前記区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
して機能させるための高速演算プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１７】
設定データである入力データと制御対象のシステム出力とを入力するデータ入力部と、該データ入力部から得られる入力データとシステム出力との差分を計算する差分演算部と、該差分演算部から得られる差分と制御器の分数階微分要素を有する伝達関数とから制御器出力値を計算する制御器出力演算部と、該制御器出力演算部から得られる制御器出力値を制御対象に出力する制御器出力用出力部とを備えており、前記制御器出力用出力部から得られる制御器出力値と前記制御対象の伝達関数とからシステム出力を出力して前記制御対象の被制御機器を作動させる制御を、前記システム出力と次の設定データとから繰り返し実行するように構成した機器制御システムであって、
前記制御器出力演算部は、前記制御器の前記伝達関数における分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
前記高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
前記高速演算装置の前記演算部は、
前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ１進む毎に式（１２）を満たしながら更新される分割点ｃ１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１２】
2013024810.tif 000138
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
前記変数ｔが幅ｈ更新される度に、前記区間［ｃ１，ｔ］又は前記区間［ｔ，ｃ１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
前記区間［ａ，ｃ１］又は前記区間［ｃ１，ａ］において、前記変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、前記分割点ｃ１を上端とする幅Δ１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ１］又は該区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
前記区間［ｃ１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ１］の分数階微積分値Ａと前記区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えていることを特徴とする機器制御システム。
【請求項１８】
物理現象が時間を変数とした分数階微積分要素を含む分数階微積分式を用い、該分数階微積分式が積分式で表現されるようにモデル化し、前記分数階微積分要素から分数階微積分値を計算する第１の演算手段と、前記分数階微積分値を用いて運動方程式の解を計算する第２の演算手段と、物理現象の測定データと前記運動方程式の経時的な解とを用いて前記物理現象の経時的な特性値を計算する第３の演算手段とを備え、物理現象の挙動を解析するシミュレーションシステムであって、
前記第１の演算手段は、前記分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより前記分数階微積分要素の前記変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
前記高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、前記変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、前記入力部及び前記記憶部に電気的に接続されている演算部とを備えており、
前記高速演算装置の前記演算部は、
前記積分区間［ａ，ｔ］又は［ｔ，ａ］内において、前記変数ｔがΔ１進む毎に式（１３）を満たしながら更新される分割点ｃ１を設定することにより、前記積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は前記積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１３】
2013024810.tif 000139
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
前記変数ｔが幅ｈ更新される度に、前記区間［ｃ１，ｔ］又は前記区間［ｔ，ｃ１］において、前記入力部を介して入力された前記データｆ（τ）を用い前記変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
前記区間［ａ，ｃ１］又は前記区間［ｃ１，ａ］において、前記変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、前記分割点ｃ１を上端とする幅Δ１の小区間における、前記変数ｔの更新に依存しない重み付き積分値を計算して前記記憶部に記憶させる積分演算手段と、
前記記憶部に記憶されている前記重み付け積分値を読み出し、前記各小区間に対応する忘却関数を乗じた値を前記小区間毎に算出し、全小区間について合算して該区間［ａ，ｃ１］又は該区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
前記区間［ｃ１，ｔ］の分数階微積分値Ａと前記区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は前記区間［ｔ，ｃ１］の分数階微積分値Ａと前記区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えていることを特徴とするシミュレーションシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分数階微積分値を計算する高速演算装置、高速演算プログラム及び高速演算プログラムを記録した記録媒体、この高速演算装置を備えた機器制御システム、並びにこの高速演算装置を備えたシミュレーションシステムに関する。
【背景技術】
【０００２】
工学設計や現象解析、装置の制御などの作業において、コンピュータを用いた数値計算やシミュレーション、さらにはそれを応用した製品設計が多くの分野で盛んに実践されている。特に、近年、コンピュータの飛躍的な性能向上により、従来の単純化したモデル式から、より精度に優れたモデル式でのシミュレーションや制御方法が研究されている。その一手法として、従来の整数階の微積分要素でのモデルに代えて、非整数階である分数階微積分要素を含むモデルの適用が有効であり、例えば、振動現象などの挙動解析や、信号処理、画像処理など多くの分野で研究がされている。【０００３】
分数階微分は、Ｎｅｗｔｏｎ力学で使用される整数階微分を一般化したものであり、現在では工学分野のみならず多くの自然科学分野でも使用されつつある。しかし、分数階微分は、整数階微分に比べて解析が困難なことから、一部では１以上の導関数を含む整数階微分モデルに変換して解析する試みも行われている。【０００４】
分数階微分は、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義、及びＣａｐｕｔｏの定義で代表されるように、積分式で記述されるため、コンピュータによる分数階微分値の計算は、積分値から求める手法が用いられている。また、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義やＣａｐｕｔｏの定義のように分数階微分は、分数階積分の形で定義されることから明らかなように、分数階微分の計算手法は、分数階積分値の計算にも適用できる。【０００５】
積分計算としては、従来から台形則やシンプソン則など公知の方法が使用されている。これら従来の方法は、積分計算に際して積分間隔を細かく分割して計算するため、高い精度の計算値を得るには変数軸の分割点を多くする必要がある。しかし、分割点の数が多くなるほど計算時間が増加し、計算精度を求めるほど処理時間が膨大となる。例えば、Ｒｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義における閉区間［ａ，ｔ］についての積分において、変数の点（変数が時間の場合には時刻）τについての積分は、変数ｔに対する相対値（ｔ−τ）で表された変数ｔに依存の積分であるため、従来の分数階微分値や分数階積分値の計算では、全ての分割点のｆ（τ）をメモリに記憶しておき、変数ｔが進む毎にメモリのｆ（τ）を用いて更新前の積分計算を行う必要があることから、第ｎステップまでの積分計算ではｎ２オーダーの計算が必要となる。これは、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義、及びＣａｐｕｔｏの定義など、積分形式を含む場合にも共通する課題である。【０００６】
このような課題を解消する方法として、整数階微分の計算法であるＮｅｗｍａｒｋ−β法を分数階微分に拡張して応用する方法が提案されている（非特許文献１）。【０００７】
前述の課題を解消するために他の方法として、変数ｔを時間として現時刻から過去の時間領域の分割幅が広くなるようにベキ時間を適用して計算ステップ数を減じることにより、計算時間の短縮化を試みる提案がされている（非特許文献２）。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】那須野洋，清水信行：「非線形分数階微分方程式のための数値積分アルゴリズム」，日本機械学会論文集（Ｃ編）７２巻７２２号（２００６−１０）論文Ｎｏ．０６−１３７
【非特許文献２】那須野洋，清水信行：「ベキ時間を用いた分数階微分方程式の数値積分法」，日本機械学会論文集（Ｃ編）７２巻７２４号（２００６−１２）論文Ｎｏ．０６−６０５
【非特許文献３】那須野洋，清水信行，安野拓也：「分数階微分で記述される粘弾性体の幾何学的非線形静的・動的モデル」，日本機械学会論文集（Ｃ編）７２巻７１６号（２００６−４）論文Ｎｏ．０５−１３０
【非特許文献４】Ｒ．Ｄｅｎｇ，  Ｐ．Ｄａｖｉｅｓ，  Ａ．Ｋ．Ｂａｊａｊ：「Ａｐｐｌｉｃａｔｉｏｎ  ｏｆ  Ｆｒａｃｔｉｏｎａｌ  Ｄｅｒｉｖａｔｉｖｅｓ  ｔｏ  Ｍｏｄｅｌｉｎｇ  ｔｈｅ  Ｑｕａｓｉ−Ｓｔａｔｉｃ  Ｒｅｓｐｏｎｃｅ  ｏｆ  Ｐｏｌｙｕｒｅｔｈａｎｅ  Ｆｏａｍ」，  Ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ＡＳＭＥ  ２００３  Ｄｅｓｉｇｎ  Ｅｎｇｉｎｅｅｒｉｎｇ  Ｔｅｃｈｎｉｃａｌ  Ｃｏｎｆｅｒｅｎｃｅｓ，  ＶＩＢ−４８３９７  （２００３）
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、非特許文献１に記載の方法は、処理時間が長く、実用できるものではなかった。【００１０】
また、非特許文献２に記載の方法によると、幅を広くすると計算精度が低下し、また、時刻ｔが更新される度に全区間の分割点について積分計算を実行しなければならず、依然として実用レベルの精度確保と処理時間の短縮化を実現できるものではなかった。【００１１】
このように従来方法では、精度と処理時間を両立することが困難であり、実用化において大きな課題となっているが、コストや効率を含めて実用的な計算方法は未だ開発されていない。【００１２】
分数階微積分が適用、又は期待されている具体例として、電気炉やエアコンなどの温調機能を有する機器の温度制御などの機器制御技術がある。代表的な制御技術としては、現在ＰＩＤ（比例・積分・微分）制御が適用されており、これに分数階微積分を適用したＰＩλＤμ制御等が、より正確な制御法として研究されているが、分数階微積分値の計算に際して処理時間がかかりすぎて、未だ十分な精度で制御が行えないのが現状である。【００１３】
また、挙動解析においては、例えば近年研究が進んでいる粘弾性の変形挙動解析について、粘弾性体のゴムやポリウレタンフォームなどでは、試料の荷重試験による周波数応答からその粘弾性を推定する研究が数多く報告されているが、より正確な予測、解析を目的に、分数階微分モデルを用いた粘弾性の推定手法（非特許文献３）が提案され、また、実際に分数階微分モデルを用いてポリウレタンフォームの粘弾性を推定すること（非特許文献４）などが報告されている。【００１４】
しかしながら、この分数階微積分からなる要素の計算方法は、台形則やシンプソン則やＧｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖ法など様々な方法が提案されているが、精度を求めるほど処理時間が膨大となり、精度と処理時間を両立させる大きな課題となっており、コストや効率を含めて実用的な計算方法は未だ開発されておらず、分数階微積分モデルを用いた解析や製品設計の実用化できない原因となっていた。【００１５】
このようなことから、高精度かつ短時間で処理可能な分数階微積分値の計算方法が開発されれば、様々な分野での実用化が進み、産業的に大きな貢献をもたらすため、コンピュータを用いた革新的な分数階微積分値の高速演算装置の開発が切望されていた。【００１６】
従って、本発明の目的は、上記問題点に鑑み、高い精度を維持しつつ短い計算時間で分数階微積分値を計算することができる高速演算装置、高速演算プログラム及び高速演算プログラムを記録した記録媒体、機器制御システム、並びにシミュレーションシステムを提供することにある。
【課題を解決するための手段】
【００１７】
本発明者等は、上記問題点を解決するために、分数階微分をＲｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義（以下ＲＬ定義と称する）、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義（以下ＧＬ定義と称する）、及びＣａｐｕｔｏの定義（以下Ｃ定義と称する）による式を用いて分数階微分値及び／又は分数階積分値を求める際に、ｆ（τ）の区間［ａ，ｔ］の積分計算を行うにあたり、区間［ａ，ｔ］をその区間内の点ｃ（以下分割点と称する）で区間［ａ，ｃ］と区間［ｃ，ｔ］とにまず分割し、区間［ａ，ｃ］の積分を特定の条件の下で、（ｔ−τ）のベキをｔのベキとτのベキの積に変換することによって、従来のｆ（τ）に代えて、変数ｔと独立なτによるｆ（τ）の重み付き積分を用いて、ｆ（τ）の重み付き積分項と、変数ｔのベキに依存する重み付き係数との積の和として計算するように構成した高速演算装置を提供する。これにより、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させておき、その値を読み出して使用することによって、その都度、最初から再計算する必要がないため、さらに、全ての計算ステップに対する関数値ｆ（τ）を記憶部上に記憶しておく必要もなくなるため、極めて高速に分数階微積分値計算を行うことができる。また、分数階微分の定義に分数階積分が含まれていることから、分数階微分の場合と同様に分数階積分値も極めて高速に計算することができるのである。【００１８】
即ち、本発明によれば、高速演算装置は、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求めるように構成されており、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えている。演算部は、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ１進む毎に式（１）を満たしながら更新される分割点ｃ１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する第１の分割設定手段と、
【数１】
2013024810.tif 000003
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ１，ｔ］又は区間［ｔ，ｃ１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ１］又は区間［ｃ１，ａ］において、変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、分割点ｃ１を上端とする幅Δ１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ１］又は区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ１］の分数階微積分値Ａと区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えている。【００１９】
変数ｔに関する積分区間が［ａ，ｔ］であって、
第１の分割設定手段は、式（２）の条件を満たす設定値Ｍ及びｍにおいて、ｔ−（Ｍ＋ｍ）ｈ＝ｃ１となった時点を分割点ｃ１として設定するものであり、
積分演算手段は、区間［ａ，ｃ１］の変数ｔの更新に依存しない重み付き積分値の算出を、分割点ｃ１が設定された時点からｃ１＝ａとして開始し、以降、変数ｔがｍｈ（＝Δ１）進みｃ１＝ｔ−（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｃ１，ｔ−Ｍｈ］の幅ｈ刻みのｆ（τ）データを用い、ｃ１＝ｔ−Ｍｈへの更新を伴って、更新されたｃ１を上端とする幅Δ１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して記憶部に記憶させる第１の積分演算部を備えている
【数２】
2013024810.tif 000004
（ただし、εＴは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
ことが好ましい。【００２０】
変数ｔに関する積分区間が［ｔ，ａ］であって、
第１の分割設定手段は、式（３）の条件を満たす設定値Ｍ及びｍにおいて、ｔ＋（Ｍ＋ｍ）ｈ＝ｃ１となった時点を分割点ｃ１として設定するものであり、
積分演算手段は、区間［ｃ１，ａ］の変数ｔの更新に依存しない重み付き積分値の算出を、分割点ｃ１が設定された時点からｃ１＝ａとして開始し、以降、変数ｔがｍｈ（＝Δ１）進みｃ１＝ｔ＋（Ｍ＋ｍ）ｈの条件を満たす度に、区間［ｔ＋Ｍｈ，ｃ１］の幅ｈ刻みのｆ（τ）データを用い、ｃ１＝ｔ＋Ｍｈへの更新を伴って、更新されたｃ１を上端とする幅Δ１の小区間の重み付け積分値を算出する重み付き計算ステップを実行して記憶部に記憶させる第１の積分演算部を備えている
【数３】
2013024810.tif 000005
（ただし、εＴは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
ことも好ましい。【００２１】
第１の分割設定手段によって、式（４）の条件を満たす設定値Ｍ及びｍにおいて、変数ｔのｈ刻みの更新に伴い分割点ｃ１がｃ１＝ｔ−Ｍｈとして更新され、さらに、ｃ１＝ａを初期値として、区間［ａ，ｃ１］内のｃ１を上端とする幅Δ’が変数ｔのｈ刻みの更新によってｈからｍｈまで変化するものであって、
【数４】
2013024810.tif 000006
（ただし、εＴは計算精度、Ｍは時間刻み幅ｈ区間の積分計算のステップ数、ｍはΔ１に纏められる幅ｈ刻みの区間の個数であり２以上の整数、Ｋは展開次数）、
積分演算手段は、幅Δ’がｍｈ未満の場合には、ｃ１を上端とする幅Δ’の小区間は幅ｈ刻みのｆ（τ）データを用い、変数ｔの更新に依存する積分値の計算を変数ｔの更新毎に繰り返し実行して記憶部に記憶させ、幅Δ’がｍｈとなった場合には、小区間の幅ｈ刻みのｆ（τ）データを用い、分割点ｃ１を上端とする幅Δ１（＝ｍｈ）の小区間の変数ｔの更新に依存しない重み付け積分値を算出する重み付き計算を変数ｔの更新毎に繰り返し実行して記憶部に記憶させる第１の積分演算部を備えており、
第２の分数階微積分値演算手段は、幅Δ’がｍｈとなった場合に分割点ｃ１に対応する重み付き積分値を記憶部から読み出して、忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して算出する分数階微積分値と、変数ｔの更新に依存する幅Δ’がｍｈ未満の場合の小区間の積分値から算出する分数階微積分値とを合算して算出する第１の分数階微積分値演算部を備えていることも好ましい。【００２２】
積分演算手段は、変数ｔを更新しながら区間［ａ，ｃ１］の重み付き積分を算出する際に、幅Δ１の小区間の個数Ｊｓ−１（ｓ＝２）が式（５）で表されるＬｓと等しくなった時点で、ｃ２＝ｃ１−（Ｌｓ−ｍ１）Δ１である新たな分割点ｃ２を設定して区間［ａ，ｃ１］を中区間［ａ，ｃ２］と中区間［ｃ２，ｃ１］とに分割する第２の分割設定手段を具備し、
中区間［ａ，ｃ２］においては、ｃ２の設定された時点並びに以降、変数ｔがｍ１Δ１進む度に、分割点ｃ２から直近の分割点ｃ１方向に連続するｍ１個の幅Δ１の小区間を纏めて幅Δ２の小区間とし、記憶部に記憶されているｍ１個の各小区間の重み付き積分値を重み付き積分値として算出して記憶部に記憶させる第２の積分演算部を備えており、
【数５】
2013024810.tif 000007
第２の分数階微積分値演算手段は、区間［ａ，ｃ１］の分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ２，ｃ１］と、中区間［ａ，ｃ２］の各小区間の重み付け積分値、をそれぞれ読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。【００２３】
積分演算手段は、変数ｔを更新しながら区間［ｃ１，ａ］の重み付き積分を算出する際に、幅Δ１の小区間の個数Ｊｓ−１（ｓ＝２）が式（６）で表されるＬｓと等しくなった時点で、ｃ２＝ｃ１＋（Ｌｓ−ｍ１）Δ１である新たな分割点ｃ２を設定して区間［ｃ１，ａ］を中区間［ｃ１，ａ］と中区間［ｃ１，ｃ２］とに分割する第２の分割設定手段を具備し、
中区間［ｃ２，ａ］においては、ｃ２の設定された時点並びに以降、変数ｔがｍ１Δ１進む度に、分割点ｃ２から直近の分割点ｃ１方向に連続するｍ１個の幅Δ１の小区間を纏めて幅Δ２の小区間とし、記憶部に記憶されているｍ１個の各小区間の重み付き積分値を重み付き積分値として算出して記憶部に記憶させる第２の積分演算部を備えており、
【数６】
2013024810.tif 000008
第２の分数階微積分値演算手段は、区間［ｃ１，ａ］の分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ１，ｃ２］と、中区間［ｃ２，ａ］の各小区間の重み付け積分値、をそれぞれ読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。【００２４】
積分演算手段は、変数ｔを更新しながら区間［ａ，ｃ２］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃｓ，ｃｓ−１］の幅Δｓ−１の小区間の個数Ｊｓ−１（ｓ≧３）が式（７）で表されるＬｓと等しくなる度に、ｃｓ＝ｃｓ−１−（Ｌｓ−ｍ１）Δｓ−１である新たな分割点ｃｓによってｃｓ＋１＝ａを初期値とする中区間［ｃｓ＋１，ｃｓ］と［ｃｓ，ｃｓ−１］とを設定する第３の分割手段を具備し、ｃｓが設定された時点及びそれ以降、変数ｔがｍ１Δｓ−１進む度に、分割点ｃｓから直近の分割点ｃｓ−１方向に連続するｍ１個のΔｓ−１の小区間を纏めて幅Δｓの小区間とし、記憶部に記憶されているｍ１個のΔｓ−１の各小区間の重み付き積分値を用いてΔｓの小区間の重み付き積分値として算出して記憶部に記憶させ、さらに、以降変数ｔの更新によってΔｓの小区間の個数がＬｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、
【数７】
2013024810.tif 000009
第２の分数階微積分値演算手段は、区間［ａ，ｃ１］又は分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ２，ｃ１］及び中区間［ｃ３，ｃ２］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする中区間［ｃｓ＋１，ｃｓ］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。【００２５】
積分演算手段は、変数ｔを更新しながら区間［ｃ２，ａ］の重み付き積分を算出する際に、さらにｓ≧３（ｓ＝３，４，．．．，Ｓ）として、中区間［ｃｓ−１，ｃｓ］の幅Δｓ−１の小区間の個数Ｊｓ−１（ｓ≧３）が式（８）で表されるＬｓと等しくなる度に、ｃｓ＝ｃｓ−１＋（Ｌｓ−ｍ１）Δｓ−１である新たな分割点ｃｓによってｃｓ＋１＝ａを初期値とする中区間［ｃｓ，ｃｓ＋１］と［ｃｓ−１，ｃｓ］とを設定する第３の分割手段を具備し、ｃｓが設定された時点及びそれ以降、変数ｔがｍ１Δｓ−１進む度に、分割点ｃｓから直近の分割点ｃｓ−１方向に連続するｍ１個のΔｓ−１の小区間を纏めて幅Δｓの小区間とし、記憶部に記憶されているｍ１個のΔｓ−１の各小区間の重み付き積分値を用いてΔｓの小区間の重み付き積分値として算出して記憶部に記憶させ、さらに、以降変数ｔの更新によってΔｓの小区間の個数がＬｓとなる度にｓ＝ｓ＋１として、以上の処理を繰り返す第２の積分演算部を備えており、
【数８】
2013024810.tif 000010
第２の分数階微積分値演算手段は、区間［ｃ１，ａ］の分数階微積分値Ｂを、記憶部に記憶されている中区間［ｃ１，ｃ２］及び中区間［ｃ２，ｃ３］並びにｓ≧３（ｓ＝３，４，．．．，Ｓ）とする中区間［ｃｓ，ｃｓ＋１］の各小区間の重み付け積分値をそれぞれ読み出し、全ての中区間の各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算することによって算出する第２の分数階微積分値演算部を備えていることも好ましい。【００２６】
ｍが２以上であることも好ましい。【００２７】
ｍ１が２であることことも好ましい。【００２８】
忘却関数が式（９）であり、
積分演算手段は、式（９）の忘却関数又は式（９）の分子をあらかじめ計算し、得られたデータを記憶部に記憶させる第３の積分演算部を備えており、
第２の分数階微分値演算手段を実行する
【数９】
2013024810.tif 000011
（ただし、ｂは、個々の小区間を識別するための任意の点、αは分数階微積分の積分定義で規定される積分階数、ｇは分数階微積分の積分定義で規定される関数、ｋは展開次数番号である）
ことも好ましい。【００２９】
第３の積分演算部は、式（９）の分子をあらかじめ計算し、得られたデータを記憶部に記憶させるものであり、
積分演算手段は、さらに、重み付き積分にあらかじめ式（９）中の係数１／（ｋ！Γ（α−ｋ））を乗じた重み付き積分値として算出して記憶部に記憶させる第４の積分演算部を備えており、
第２の分数階微積分値演算手段は、記憶部に記憶されている小区間の前述の係数を乗じた重み付け積分値を読み出し、全ての中区間の各小区間に対応する式（９）の忘却関数の分子を乗じた値を小区間毎に算出し、全小区間について合算することによって区間［ａ，ｃ１］の分数階微積分値を算出する第３の分数階微積分値演算部を備えている
ことも好ましい。【００３０】
区間［ａ，ｃ１］の分数階微積分定義がＲｉｅｍａｎｎ−Ｌｉｏｕｖｉｌｌｅの定義であり、区間［ｃ１，ｔ］の分数階微積分定義がＣａｐｕｔｏの定義であることも好ましい。【００３１】
本発明によれば、さらに、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微分要素の変数ｔにおける分数階微分値を求めるコンピュータの高速演算プログラムであって、
コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
コンピュータを、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ１進む毎に式（１０）を満たしながら更新される分割点ｃ１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１０】
2013024810.tif 000012
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ１，ｔ］又は区間［ｔ，ｃ１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ１］又は区間［ｃ１，ａ］において、変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、分割点ｃ１を上端とする幅Δ１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ１］又は区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ１］の分数階微積分値Ａと区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
して機能させる高速演算プログラムが提供される。【００３２】
本発明によれば、さらにまた、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求めるコンピュータの高速演算プログラムであって、
コンピュータは、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
コンピュータを、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ１進む毎に式（１１）を満たしながら更新される分割点ｃ１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１１】
2013024810.tif 000013
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ１，ｔ］又は区間［ｔ，ｃ１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ１］又は区間［ｃ１，ａ］において、変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、分割点ｃ１を上端とする幅Δ１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ１］又は区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ１］の分数階微積分値Ａと区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
して機能させるための高速演算プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。【００３３】
本発明によれば、さらに、設定データである入力データと制御対象のシステム出力とを入力するデータ入力部と、データ入力部から得られる入力データとシステム出力との差分を計算する差分演算部と、差分演算部から得られる差分と制御器の分数階微分要素を有する伝達関数とから制御器出力値を計算する制御器出力演算部と、制御器出力演算部から得られる制御器出力値を制御対象に出力する制御器出力用出力部とを備えており、制御器出力用出力部から得られる制御器出力値と制御対象の伝達関数とからシステム出力を出力して制御対象の被制御機器を作動させる制御を、システム出力と次の設定データとから繰り返し実行するように構成した機器制御システムであって、
制御器出力演算部は、制御器の伝達関数における分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
高速演算装置の演算部は、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ１進む毎に式（１２）を満たしながら更新される分割点ｃ１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１２】
2013024810.tif 000014
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ１，ｔ］又は区間［ｔ，ｃ１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ１］又は区間［ｃ１，ａ］において、変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、分割点ｃ１を上端とする幅Δ１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ１］又は区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ１］の分数階微積分値Ａと区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えている機器制御システムが提供される。【００３４】
本発明によれば、さらに、物理現象が時間を変数とした分数階微積分要素を含む分数階微積分式を用い、分数階微積分式が積分式で表現されるようにモデル化し、分数階微積分要素から分数階微積分値を計算する第１の演算手段と、分数階微積分値を用いて運動方程式の解を計算する第２の演算手段と、物理現象の測定データと運動方程式の経時的な解とを用いて物理現象の経時的な特性値を計算する第３の演算手段とを備え、物理現象の挙動を解析するシミュレーションシステムであって、
第１の演算手段は、分数階微積分要素の分数階微積分式を積分式で定義し、変数ｔに関する積分区間［ａ，ｔ］又は逆向きの［ｔ，ａ］について、刻み幅ｈの点τ（区間［ａ，ｔ］においてはａ≦τ≦ｔ、区間［ｔ，ａ］においてはｔ≦τ≦ａ）におけるデータｆ（τ）を用いて積分計算することにより分数階微積分要素の変数ｔにおける分数階微積分値を求める高速演算装置を備えており、
高速演算装置は、パラメータ及びデータｆ（τ）を入力するための入力部と、変数ｔの更新に依存しない重み付き積分値を少なくとも記憶する記憶部と、入力部及び記憶部に電気的に接続されている演算部とを備えており、
高速演算装置の演算部は、
積分区間［ａ，ｔ］又は［ｔ，ａ］内において、変数ｔがΔ１進む毎に式（１３）を満たしながら更新される分割点ｃ１を設定することにより、積分区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割、又は積分区間［ｔ，ａ］を区間［ｔ，ｃ１］と区間［ｃ１，ａ］とに分割する分割設定手段と、
【数１３】
2013024810.tif 000015
（ただし、Δ１は分割点ｃ１を上端とする変数の小区間の幅（Δ１＝ｍｈ、ｍは２以上の整数）、εは計算精度が確保されるための上限値（０＜ε＜１））、
変数ｔが幅ｈ更新される度に、区間［ｃ１，ｔ］又は区間［ｔ，ｃ１］において、入力部を介して入力されたデータｆ（τ）を用い変数ｔに依存する積分値を算出して分数階微積分値Ａを求める第１の分数階微積分値演算手段と、
区間［ａ，ｃ１］又は区間［ｃ１，ａ］において、変数ｔがΔ１更新される度に、分割点ｃ１を端点とする幅Δ１の小区間内のｆ（τ）を用いて、分割点ｃ１を上端とする幅Δ１の小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶部に記憶させる積分演算手段と、
記憶部に記憶されている重み付け積分値を読み出し、各小区間に対応する忘却関数を乗じた値を小区間毎に算出し、全小区間について合算して区間［ａ，ｃ１］又は区間［ｃ１，ａ］の分数階微積分値Ｂを算出する第２の分数階微積分値演算手段と、
区間［ｃ１，ｔ］の分数階微積分値Ａと区間［ａ，ｃ１］の重み付け分数階微積分値Ｂとを合算して区間［ａ，ｔ］の分数階微積分値を算出、又は区間［ｔ，ｃ１］の分数階微積分値Ａと区間［ｃ１，ａ］の重み付け分数階微積分値Ｂとを合算して区間［ｔ，ａ］の分数階微積分値を算出する全分数階微積分値算出手段と
を備えているシミュレーションシステムが提供される。
【発明の効果】
【００３５】
本発明の高速演算装置による分数階微積分値の計算手法によれば、分数階微分の定義式中の積分計算で記憶部に記憶するのはｆ（τ）ではなく、変数ｔに依存しない重み付き積分Ｉｓ，ｋｆ（τ）であるので、記憶部に記憶させるデータ数が低減されて記憶部容量を小さくでき、変数ｔが更新される毎に全積分区間の積分を再計算する必要がなく、その重み付き積分値を読み出して使用すればよいため、また、計算ステップも低減されるので計算時間を大幅に短縮した高速化を行うことができる。さらに、分数階微分の定義に分数階積分が含まれていることから、分数階微分の場合と同様に分数階積分値も極めて高速に計算することができる。【００３６】
しかも、その計算誤差も小さいので、高精度かつ高速に分数階微積分値を求めることができる。従って、分数階微積分要素を含んでモデル化された様々な現象、事象の解析について、精度よく、かつ高速に結果が得られるので、例えば、ＰＩＤ制御を組み込んだ装置の制御をはじめ、様々な分野の装置の設計や解析などに活用できる。また、本発明の高速演算プログラム及びこの高速演算プログラムを記録したコンピュータ読み取り可能な記録媒体は、プロセス量のＰＩＤ制御や各種物理現象の経時的挙動のシミュレーション、それらを応用した製品設計などを効率的に実施することができる。【００３７】
さらに、機器制御システムは、その分数階微積分値計算が上述の高速演算装置によって高速で行われるため、高い精度を確保しつつかつ高速な機器制御を行うことができる。また、シミュレーションシステムも、その分数階微積分値計算が上述の高速演算装置によって高速で行われるため、高い精度を確保しつつかつ高速なシミュレーションを行うことができる。
【図面の簡単な説明】
【００３８】
【図１】本発明の第１の実施形態として、高速演算装置のハードウェア構成を概略的に示すブロック図である。
【図２ａ】図１の実施形態における高速演算装置の一部動作を説明するフローチャートである。
【図２ｂ】図１の実施形態における高速演算装置の一部動作を説明するフローチャートである。
【図２ｃ】図１の実施形態における高速演算装置の一部動作の別形態を説明するフローチャートである。
【図３】図１の実施形態においてソフトウェアによって構築される高速演算装置の主要部の構成を概略的に示すブロック図である。
【図４】図１の実施形態における高速演算装置において、分数階微分値を計算する場合に、区間［０，ｔ］内のｃ点で、区間［０，ｃ］と区間［ｃ，ｔ］とに分割することを説明するための説明図である。
【図５】分数階微分値を計算する場合に、微分変数の置き換えを説明するための説明図である。
【図６ａ】分数階微分値を計算する場合の計算ステップ１を説明するための模式図である。
【図６ｂ】分数階微分値を計算する場合の計算ステップ２を説明するための模式図である。
【図６ｃ】分数階微分値を計算する場合の計算ステップ３を説明するための模式図である。
【図６ｄ】分数階微分値を計算する場合に、ｓ＝２としたときの計算ステップ４を説明するための模式図である。
【図６ｅ】分数階微分値を計算する場合の計算ステップ４を説明するための模式図である。
【図７】分割点ｃ１＝ｔ−Ｍｈに固定する条件で計算するアルゴリズムにより分数階微分値を計算する方法の説明図である。
【図８】本発明の第２の実施形態として、ＰＩＤ制御システムのハードウェア構成を概略的に示すブロック図である。
【図９】図８の実施形態におけるＰＩＤ制御システムの一部動作を説明するフローチャートである。
【図１０】図８の実施形態においてソフトウェアによって構築されるＰＩＤ制御システムの主要部の構成を概略的に示すブロック図である。
【図１１】本発明の第３の実施形態において、粘弾性体の変形挙動シミュレーションを行う場合を説明するためのモデル図である。
【図１２】本発明の第３の実施形態として、粘弾性解析を行うシミュレーション装置の一部動作を説明するフローチャートである。
【図１３】図１２の実施形態においてソフトウェアによって構築されるシミュレーション装置の主要部の構成を概略的に示すブロック図である。
【図１４】本発明の高速演算装置における計算方法と従来の計算方法との効果検証結果を対比するためのグラフである。
【図１５】本発明の高速演算装置における計算方法と従来の計算方法との効果検証結果を対比するためのグラフである。
【発明を実施するための形態】
【００３９】
（Ｉ）第１の実施形態（高速演算装置）
図１は本発明の第１の実施形態として、分数階微積分を求める高速演算装置のハードウェア構成を概略的に示している。なお、分数階微積分の計算手法は、分数階積分値の計算にも適用できるので、以下分数階微分の場合について説明する。なお、分数階積分の計算の場合には、以下の分数階微分の場合の説明において、「分数階微分」を「分数階積分」と置き換えればよい。
【００４０】
図１に示すように、本実施形態における高速演算装置は、コンピュータ１０と、このコンピュータ１０に接続される外部装置１１とから構成されている。
【００４１】
コンピュータ１０は、バス１０ａを介して互いに接続された中央処理装置（ＣＰＵ）１０ｂと、リードオンリメモリ（ＲＯＭ）１０ｃと、ランダムアクセスメモリ（ＲＡＭ）１０ｄと、ハードディスク駆動装置（ＨＤＤ）１０ｅと、音声処理部１０ｆと、画像処理部１０ｇと、ブルーレイディスク／デジタルバーサタイルディスク／コンパクトディスク（ＢＲ／ＤＶＤ／ＣＤ）に関してデータ等の読み書きを行うディスク駆動装置１０ｈと、入出力インタフェース１０ｉとを主に備えたコンピュータ及びこれを作動させるプログラムから構成される。
【００４２】
音声処理部１０ｆはスピーカ１０ｊに接続されており、画像処理部１０ｇはディスプレイ１０ｋに接続されている。ディスク駆動装置１０ｈはＢＲ／ＤＶＤ／ＣＤが装着可能となっており、入出力インタフェース１０ｉにはキーボード１０ｌ、マウス１０ｍ及びプリンタ１０ｎ等が接続されている。さらに、この入出力インタフェース１０ｉには、外部装置１１が接続されている。
【００４３】
ＣＰＵ１０ｂは、ＲＯＭ１０ｃ又はＨＤＤ１０ｅに記憶されているオペレーションシステム（ＯＳ）やブートプログラム等の基本プログラムに従ってＲＡＭ１０ｄに記憶されているプログラムを実行して本実施形態の処理を行う。また、ＣＰＵ１０ｂは、ＲＯＭ１０ｃ、ＲＡＭ１０ｄ、ＨＤＤ１０ｅ、音声処理部１０ｆ、画像処理部１０ｇ、ディスク駆動装置１０ｈ、及び入出力インタフェース１０ｉの動作を制御する。
【００４４】
ＲＡＭ１０ｄは高速演算装置のメインメモリとして使用され、ＨＤＤ１０ｅやディスク駆動装置１０ｈから転送されたプログラムやデータを一時的に記憶する。また、ＲＡＭ１０ｄは、プログラム実行時の各種データが一時的に記憶されるワークエリアとしても使用される。
【００４５】
ＨＤＤ１０ｅには、プログラム及びデータがあらかじめ記憶されている。また、ＨＤＤ１０ｅは、ＲＡＭ１０ｄのワークエリアとしての容量が不足している場合に、補助的にその各種データを一時的に記憶する。
【００４６】
音声処理部１０ｆは、ＣＰＵ１０ｂの指示に従って必要な音声データを再生するための処理を行い、スピーカ１０ｊへその音声信号を出力する。
【００４７】
画像処理部１０ｇは、ＣＰＵ１０ｂの指示に従って画像データを生成する。生成された画像データは、ディスプレイ１０ｋに出力される。
【００４８】
ディスク駆動装置１０ｈは、必要な場合は、ＣＰＵ１０ｂの指示に従って、セットされたＢＲ／ＤＶＤ／ＣＤからプログラムやデータを読出し、ＲＡＭ１０ｄへ転送する。また、セットされたＢＲ／ＤＶＤ／ＣＤへプログラムやデータの書き込みをすることも可能である。
【００４９】
入出力インタフェース１０ｉは、キーボード１０ｌ、マウス１０ｍ及びプリンタ１０ｎとＣＰＵ１０ｂ又はＲＡＭ１０ｄとの間のデータのやり取りを制御する。さらに、この入出力インタフェース１０ｉは、外部装置１１とのデータの入出力を制御するように構成されている。
【００５０】
このようなハードウェア構成の高速演算装置において、ＣＰＵ１０ｂは、作動時は、まず、ＲＡＭ１０ｄ内にプログラム記憶領域、データ記憶領域及びワークエリアを確保し、ＨＤＤ１０ｅやディスク駆動装置１０ｈからプログラム及びデータを取り込んで、プログラム記憶領域及びデータ記憶領域に格納する。次いで、このプログラム記憶領域に格納されたプログラムに基づいて、図２ａ〜図２ｃに示す処理を実行する。ＣＰＵ１０ｂがプログラムを実行することによって、図３に概略的に示すごとき構成の主要部を備えた高速演算装置が構築される。
【００５１】
即ち、本実施形態の高速演算装置は、分数階微分を高速で実行可能な演算装置であり、積分下限ａを０とした例として、図３に示すように、分数階微分値を計算するために必要な各種の一次パラメータや、計算対象となる各変数の関数値データを入力する入力部３０と、入力された一次パラメータやデータから二次パラメータを算出する二次パラメータ演算部３１と、区間［ｃ１，ｔ］（ただし、ｃ１は分割点）の分数階微分値を算出する第１の分数階微分値演算部３２と、区間［０，ｃ１］を構成する各小区間Δｓの入力データｆ（τ）に重み付けした変数ｔに依存しない（変数ｔと独立な）τによる重み付き積分（以下、単に重み付き積分と称す）した値を算出する重み付き積分演算部３３と、重み付き積分値を用いて区間［０，ｃ１］の分数階微分値を算出する第２の分数階微分値演算部３４と、全区間［０，ｔ］の分数階微分値を合算により求める全分数階微分値算出部３５と、必要に応じて設けられ、計算結果を出力する出力部３６と、入力部３０及び各演算部とデータの保存及び／又は読み出し作業がなされる構成となっており、高速解析プログラム、入力データ及び各演算データを記憶する記憶部３７とを備えるように構築される。
【００５２】
ここで、図１のハードウェア構成と対比すると、入力部３０は、入出力インタフェース１０ｉを介してＣＰＵ１０ｂに接続されている外部装置１１及びキーボード１０ｌ、又はディスク駆動装置１０ｈ等から構成され、二次パラメータ演算部３１、第１の分数階微分値演算部３２、重み付き積分演算部３３、第２の分数階微分値演算部３４及び全分数階微分値算出部３５は、ＣＰＵ１０ｂから構成され、出力部３６は入出力インタフェース１０ｉを介してＣＰＵ１０ｂに接続されている外部装置１１及びプリンタ１０ｎ、又はディスク駆動装置１０ｈ等から構成され、記憶部３７は、ＲＡＭ１０ｄ及びＨＤＤ１０ｅ等から構成される。
【００５３】
図２ａ、図２ｂ及び図２ｃは本実施形態における高速演算装置の主要部の動作の一部を説明している。以下、これら図２ａ、図２ｂ及び図２ｃのフローチャートに基づいて説明を行うが、本発明の原理及び技術思想を適用しかつ本発明の効果を得ることができる構成であれば、これらフローチャートに限定されるものではない。
【００５４】
以下、積分区間が［ａ，ｔ］（ここでａ＜ｔ）において、変数ｔが時間であって、積分下限値ａが０であり、ｔの進む方向がａからｔの方向（正の方向という）である場合に、ＲＬ定義の分数階微分値を演算する手順について説明する。この演算は、図２ａに示す以下のステップＳ１〜Ｓ１１を実行することによって行われる。
【００５５】
（ステップＳ１）：分数階微分値計算用一次パラメータの入力（又は設定）
まず、一次パラメータとして、ｈ、ε、ｍ、ｍ１、Ｍ、Ｋ、Ｓｍａｘ及び計算終点条件が、入力部１０ａを介して高速演算装置に入力される。ただし、ｈは積分計算時の時間ｔの刻み幅、εは許容誤差、ｍは幅Δ１の小区間に置き換えられる区間幅ｈの小区間（Δ０）の個数（又は積分ステップの個数）、ｍ１は小区間Δｓに纏められる小区間Δｓ−１の個数、Ｍは最初の分割点ｃ１までの区間幅ｈの積分計算する回数（ステップ数）、Ｋはテイラー展開次数、Ｓｍａｘは時間区間［ｃｓ＋１，ｃｓ］の最大ステップ数である。計算終点条件はパラメータを定数値とし、処理するプログラムに応じて、時刻、ステップ数、及びデータ個数等を適宜設定することができる。
【００５６】
（ステップＳ２）：読み込んだ一次パラメータから二次パラメータの計算（又は設定）
入力設定された一次パラメータから分数階微分値計算に必要な二次パラメータとして、Ｌ０、Ｌｓ、Δ１、Δｓ、ｍａｘｓｔｅｐが二次パラメータ演算部３１により式（１４）〜式（１８）を用いて計算され、記憶部３７に記憶される。ただし、Ｌ０は［０，ｃ１］区間の計算に移行するデータ個数、Ｌｓは第ｓ番目の中区間を構成する幅Δｓの小区間の最大個数、Δ１は第１番目の中区間を構成する小区間幅、Δｓは第ｓ番目の中区間を構成する小区間幅、ｍａｘｓｔｅｐはタイムステップ数の上限値である。
【数１４】
2013024810.tif 000016
【数１５】
2013024810.tif 000017
【数１６】
2013024810.tif 000018
【数１７】
2013024810.tif 000019
【数１８】
2013024810.tif 000020
【００５７】
（ステップＳ３）：初期条件の設定
計算開始時の初期条件として、パラメータｉ、Ｊｓ、及びｃｓに初期値０が与えられる。即ち、ｉ＝０、Ｊ０＝０，Ｊ１＝０，．．．，Ｊｓ＝０，．．．，ＪＳｍａｘ-１＝０、ｃ１＝０，ｃ２＝０，．．．，ｃｓ＝０，．．．，ｃＳｍａｘ＝０。ただし、ｉはタイムステップ数、Ｊｓ（ｓ＝０，１，２，．．．，Ｓｍａｘ−１）は第ｓ番目の中区間のメモリ数、ｃｓ（ｓ＝１，２，．．．，Ｓｍａｘは第ｓ番目の中区間の上端である。
【００５８】
（ステップＳ４〜ステップＳ７）：区間［ｃ１，ｔ］の分数階微分値の計算
第１の分数階微分値演算部３２により、区間［ｃ１，ｔ］の分数階微分値が算出される。まず、計算区間［ｃ１，ｔ］（ｓ＝０番目の中区間）のメモリ数Ｊ０がＪ０＝Ｌ０となるまで（ステップＳ７においてＮｏと判断されている間）は、タイムステップｉとメモリ数Ｊ０とを１つずつ歩進させながら（ステップＳ４）、時間幅ｈ刻みのデータｆｉ＝ｆ（τｉ）を入力部３０を介して例えば外部装置１１（図１）から読み込み（ステップＳ５）、台形則などの従来法で［ｃ１，ｔ］の分数階微分値Ａを算出し、算出した分数階微分値Ａを記憶部３７に記憶させる（ステップ６）。この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ１に相当する。なお、ステップＳ５におけるデータｆ（τｉ）の読み込みは、区間［０，ｔ］の全データをあらかじめ入力して記憶部３７に記憶させたものを読み出す方式や、信号制御等のようにリアルタイムで時間更新時のデータを外部装置１１等から順次取り込む方式等、目的に応じて適宜適用することができる。
【００５９】
（ステップＳ８）：区間［０，ｃ１］の小区分の重み付き積分の計算
メモリ数Ｊ０がＪ０＝Ｌ０となった時点で（ステップＳ７においてＹｅｓと判断されたとき）、区間［０，ｃ１］の分数階微分値Ｂの計算に移行する。まず、ステップＳ８において、重み付き積分演算部３３により、区間［０，ｃ１］の分数階微分値Ｂの計算に用いる、区間［０，ｃ１］を構成する幅Δｓの小区間（ｓ＝１，２，．．．，Ｓｍａｘ−１）の重み付き積分値が算出される。重み付き積分値の計算は、まず、時刻ｔが更新されてＪ０＝Ｌ０となる度にｓ＝１である時刻ｔにおけるｃ１を上端とするΔ１幅の小区間についてのみ実行され、ｃ１は重み付き積分の計算が実行される度にｔ側にΔ１ずつ移動し、移動後の小区間［ｃ１−Δ１，ｃ１］に固有の重み付き積分値として記憶部３７に記憶される。さらに、ｓ≧２以上に移行して区間［０，ｃ１］内を幅Δｓの小区間（Δｓ＝ｍ１Δｓ−１）からなる中区間［ｃｓ＋１，ｃｓ］で構成し、前述のｃ１をｃｓ、前述のΔ１をΔｓとして同様の処理を実行する。こうして固有の重み付き積分値を有する各小区間［ｃｓ−（ｊ＋１）Δｓ，ｃｓ−ｊΔｓ］同士は、重複することなく時系列で連続して並んで区間［０，ｃ１］を構成し、各小区間の重み付け積分値は、時刻ｔに依存しないため、一度計算して記憶させておけば、時刻ｔの更新毎に再計算する必要はないため、非常に高速な計算処理を行うことができる。このステップＳ８における計算の詳細については、サブルーチンフローとして後述する。
【００６０】
（ステップＳ９）：区間［０，ｃ１］の分数階微分値Ｂの計算
記憶部３７に記憶されている区間［０，ｃ１］を構成する小区間の重み付き積分値を用い、式（１９）又は式（２０）（ｓ＝１固定の場合）によって分数階微分値Ｂが計算される。式（１９）及び式（２０）の内容については後述する。応用例として、あらかじめ入力して記憶しておいた式（１９）又は式（２０）中の右辺中の（ｔ−ｃｓ＋ｊΔｓ）−ｑ−ｋ−１のベキ計算値を読み出しながら式（１９）又は式（２０）の計算を実行してもよい。
【数１９】
2013024810.tif 000021
【数２０】
2013024810.tif 000022
【００６１】
（ステップＳ１０）：［０，ｔ］の分数階微分値の計算
全分数階微分値算出部３５により、分数階微分値Ａと分数階微分値Ｂとを合算して［０，ｔ］の分数階微分値（時刻ｔにおける分数階微分値）が算出される。
【００６２】
（ステップＳ１１）：計算終点判断
次いで計算終点判断が行われ、終点条件に達しない場合は、ステップＳ４へ戻り、タイムステップｉとメモリ数Ｊ０とを１つずつ歩進（時間更新に相当）させながら計算終点条件までステップＳ４〜Ｓ１１を繰り返して実行することにより、目的の時刻ｔの分数階微分値を算出し、終点条件に達した時点（ステップＳ１１においてＹｅｓと判断したとき）にこの処理を終了する。なお、リアルタイム制御のように、分数階微分値計算のフローが制御フローのサブルーチンとして適用される場合には、制御フローの計算終点条件によって分数階微分値計算フローが終了される。
【００６３】
次に、ステップＳ８の処理内容を図２ｂに示すサブルーチンフローを用いて詳細に説明する。
【００６４】
まず、区間［０，ｃ１］に移行した時点でｓ＝１が与えられ（ステップＳＳ１）、幅Δ１の小区間［ｃ１，ｃ１＋Δ１］のｍ個のｈ区間のデータｆ（τ）が記憶部３７から読み出され、その重み付き積分値Ｉ１，ｋｆ（ｃ１＋Δ１）が積分上端をｍｈとした式（２１）によって算出され、又はｍ＝２の場合には式（２１）は、式（２２）によって算出される（ステップＳＳ２）。これら式（２１）及び式（２２）の内容については後述する。次いで、ｃ１＝ｃ１＋Δ１として、分割点ｃ１を現時刻ｔ側にΔ１だけ移動させ、重み付き積分値Ｉ１，ｋｆ（ｃ１＋Δ１）が重み付き積分値Ｉ１，ｋｆ（ｃ１）として記憶部３７に記憶される（ステップＳＳ３）。この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ２に相当する。
【数２１】
2013024810.tif 000023
【数２２】
2013024810.tif 000024
【００６５】
次いで、Ｊ１をＪ１＋１と歩進させ、合わせて区間［ｃ１，ｔ］の幅ｈの区間数に相当するメモリ数Ｊ０を重み付き積分として纏めたｍ個分を減じてＪ０＝Ｊ０−ｍとする（ステップＳＳ４）。このステップＳＳ４の処理は、時間更新されて図２ａのステップＳ４に戻ったときに、区間［ｃ１，ｔ］の幅ｈの区間数に相当するメモリ数Ｊ０がＭ個からステップＳ６の分数階微分値Ａの計算が再開され、ｍステップ進むとステップＳ８のサブルーチンに移行するように実行するものである。次いで、幅Δ１の小区間数に相当するメモリ数Ｊ１がＬｓに達したかが判断され（ステップＳＳ５）、Ｎｏと判断された場合はサブルーチンを完了して図２ａのステップＳ９へ移行する。
【００６６】
さらに時間更新して図２ａのステップＳ４〜Ｓ１１の処理が繰り返し実行される流れの中で、その都度、図２ｂのサブルーチンにおけるステップＳＳ１〜ＳＳ５の処理が繰り返され、各小区間に固有の重み付き積分値が順次記憶部３７に記憶される。このサブルーチンにおける、この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ３に相当する。なお、図２ａのステップＳ９において、式（１９）によって分数階微分値Ｂが計算される場合においては、現時刻ｔを基準に、ｊ、ｓによって総和計算を行うため、幅Δ１の小区間の重み付き積分値はｃ１側から最新の重み付け積分値をｊ＝０番目データとして、過去に向かって順次１番目〜（Ｊ１−１）番目まで識別ラベルを時間更新毎に貼り替える処理（代入計算）をして記憶部３７に記憶される。このラベルの貼り替え処理には任意の方法が適用できる。
【００６７】
図２ｂに示すサブルーチンのステップＳＳ５において、Ｊ１＝Ｌｓであると判断された時点で（Ｙｅｓと判断したとき）、（Ｊ１−１）番目から（Ｊ１−ｍ１）番目の幅Δ１の小区間を幅Δ２の小区間に纏める処理に移行するために、ステップＳＳ６へ進む。
【００６８】
ステップＳＳ６においては、分点拡大処理を実行するかどうか判断し、実行する場合（Ｙｅｓと判断したとき）は、次のステップＳＳ７へ進む。ステップＳＳ７においては、ｓ＝ｓ＋１（＝２）として分点ｃ２が設定され、次いで、中区間［ｃ２，ｃ１］の（Ｊ１−１）番目から（Ｊ1−ｍ１）番目の幅Δ１の小区間を纏めた区間［ｃ２，ｃ２＋Δ２］（ここでΔ２＝ｍ１Δ１）の重み付き積分値が、既に計算、記憶された上述の（Ｊ１−１）番目から（Ｊ1−ｍ１）番目の小区間の重み付け積分値を用いて、ｓ＝２とした式（２３）によって算出され、ｍ１＝２の場合には式（２４）を用いて算出される（ステップＳＳ８）。これら式（２３）及び式（２４）の内容については後述する。なお、以上はｓ＝２として説明したが、本ステップでは、第（ｓ−１）中区間を構成する幅Δｓ−１の小区間のうち、（Ｊｓ−１−１）番目から第（Ｊｓ−１−ｍ１）番目の小区間の積分データを用いて区間［ｃｓ、ｃｓ+Δｓ］の重み付き積分値が式（２３）を用いて算出されるのである。
【数２３】
2013024810.tif 000025
【数２４】
2013024810.tif 000026
【００６９】
次いで、分点ｃ２をΔ２分だけ現時刻ｔ側に移動（ｃ２＝ｃ２＋Δ２として幅Δ２の小区間の上端をｃ２とする処理）し、重み付き積分値Ｉ２，ｋｆ（ｃ２＋Δ２）が重み付き積分値Ｉ２，ｋｆ（ｃ２）として記憶部３７に記憶される（ステップＳＳ９）。この処理は、後述する重み付き積分の計算アルゴリズムの説明における計算ステップ４に相当する。
【００７０】
その後、区間［ｃ３，ｃ２］中の幅Δ２の小区間の数に相当するメモリ数Ｊ２をＪ２＝Ｊ２＋１と更新すると共に、幅Δ２に纏められた幅Δ１の小区間数ｍ１個分をステップＳＳ４のＪ１から減じる処理（Ｊ１＝Ｊ１−ｍ１）が実行される（ステップＳＳ１０）。この処理により、次回のサブルーチン実行時に中区間［ｃ２，ｃ１］中の幅Δ１の小区間数Ｊ１は（Ｌｓ−ｍ１）個から開始される。次いで、幅Δ２の小区間数に相当するメモリ数Ｊ２がＬｓに達したかを判断し（ステップＳＳ１１）、Ｙｅｓの判断であれば、ステップＳＳ７に戻ってｓ＝ｓ＋１としてｓ≧３の場合に移行し、ｓ＝２の場合と同様に、ステップＳＳ７〜ＳＳ１１の処理が繰り返して実行される。
【００７１】
即ち、ステップＳＳ７において分点ｃｓが設定され、次いで第（ｓ−１）中区間を構成する幅Δｓ−１の小区間のうち、（Ｊｓ−１−１）番目から第（Ｊｓ−１−ｍ１）番目の小区間の積分データを用いて区間［ｃｓ、ｃｓ+Δｓ］の重み付き積分値が式（２３）を用いて算出される（ステップＳＳ８）。次いで、分点ｃｓをΔｓ分だけ現時刻ｔ側に移動（ｃｓ＝ｃｓ＋Δｓとして幅Δｓの小区間の上端をｃｓとする処理）し、重み付き積分値Ｉｓ，ｋｆ（ｃｓ＋Δｓ）が重み付き積分値Ｉｓ，ｋｆ（ｃｓ）として記憶部３７に記憶され（ステップＳＳ９）、その後、区間［ｃｓ＋１，ｃｓ］中の幅Δｓの小区間の数に相当するメモリ数ＪｓをＪｓ＝Ｊｓ＋１と更新すると共に、幅Δｓに纏められた幅Δｓ−１の小区間数ｍ１個分をＪｓ−１から減じる処理（Ｊｓ−１＝Ｊｓ−１−ｍ１）が実行され、次いで、幅Δｓの小区間数に相当するメモリ数ＪｓがＬｓに達したかの判断（ステップＳＳ１１）処理が繰り返される。
【００７２】
一方、ステップＳＳ１１における判断がＮｏの場合は、このサブルーチンを完了して図２ａのステップＳ９の処理へ移行し、以後、ステップＳ８の実行の度にサブルーチンＳＳ１〜ＳＳ１１が実行される。
【００７３】
このように、図２ｂのサブルーチン処理は、ｓ≧１においてＪｓがＬｓに達した中区間［ｃｓ，ｃｓ−１］について、ｃｓに直近の連続するｍ１個（ｓ＝１の場合はｍ個）の幅Δｓ−１の小区間を纏めて幅Δｓの小区間の重み付き積分値とする演算を実行する。また、ｓ−１の場合と同様に、幅Δｓの小区間の重み付き積分値（ｓ≧２）は、ｃｓ側から最新の重み付け積分値をｊ＝０番目データとして、もとからあるデータは過去に向かって順次１番目〜（Ｊｓ−１）番目まで識別ラベルを時間更新毎に張り替える処理（代入計算）をして記憶部３７に記憶される。このラベルの張り替え処理には任意の方法が適用できる。
【００７４】
また、図２ａにおけるステップＳ８のサブルーチンフローとして、図２ｃに示す他のサブルーチンフローを適用することができる。図２ｃに示す他のサブルーチンは、図２ｂのステップＳＳ２〜ＳＳ３及びステップＳＳ８〜ＳＳ９に代えて、ステップＮＳ２〜ＮＳ３及びステップＮＳ８〜ＮＳ９の処理を行うものである。即ち、先にｃ１及びｃｓの移動を含め旧重み付き積分値データの並べ替えを実行した（ステップＮＳ２及びＮＳ８）後に、新Δ１小区間の重み付き積分Ｉ１，ｋｆ（ｃ１）及び新Δｓ小区間の重み付き積分値Ｉｓ，ｋｆ（ｃｓ）（ここでｓ≧２）の計算を実行する（ステップＮＳ３及びＮＳ９）ものであり、その他のステップＮＳ１、ＮＳ４〜ＮＳ７、及びＮＳ１０〜ＮＳ１１は図２ｂのステップＳＳ１、ＳＳ４〜ＳＳ７、及びＳＳ１０〜ＳＳ１１とほぼ同様である。ただし、図２ｃのフローの方が図２ｂのフローに比べて数値計算の工程が少なくて済むため、より高速で処理することができる。
【００７５】
以上、本発明の第１の実施形態として、変数ｔを時間とし、積分下限ａをａ＝０とした高速演算装置について説明したが、本発明は、変数の変化に方向をもつものであれば時間以外の変数からなる分数階微分要素の分数階微分値の計算にも適用でき、積分下限も０に限定されない。また変数の方向も、変数を時間とした上記説明の例とは逆の現在から過去への流れ方向(負の方向という)の場合にも適用できる。時間以外の変数の例としては、分数階の空間微分モデル化した拡散現象の解析などに適用できる。また、後述する制御技術においては制御対象として、時間以外にも、温度、圧力、長さ（距離）、面積、体積、重量などがあり、これらの変化を計算し、適宜電気信号として制御器に入力し、対象物（製品、原料・反応物）の流量、タンク内のレベル、ポンプやコンプレッサーの回転数などを制御することができる。温度、圧力、長さ（距離）、面積、体積、重量のいずれかを変数とする場合は、上述の計算方法において、時間（時刻）の代わりに、これら温度などを用いることで全く同様に計算することができる。本発明の場合、変数が時間であれば、現時刻を基準として過去に遡って分割点を設けて区間を決め計算することにしているが、温度であれば、基準温度に対して低温側（あるいは高温側）に分割点を設けて区間を決め計算することになる。
【００７６】
（ＩＩ）高速演算装置における分数階微分値の計算方法
本実施形態の高速演算装置は、分数階微分要素を含むモデルにおける分数階微分要素を式（２５）のＲＬ定義、式（２６）のＧＬ定義、式（２７）のＣ定義のような積分形式の定義を用いて区間［ａ，ｔ］の分数階微分値を計算する。ただし、これらの式（２５）〜式（２７）において、正の非整数ｑに対して整数ｎｑは、ｑ＜ｎｑ＜ｑ＋１を満たすものである。また、積分方向が逆向き（変数の進みの方向が負）の場合には、式（２５）のＲＬ定義は、式（２８）、同様に式（２７）のＣ定義の場合には式（２９）となり、中区間及び小区間の並びを負の向きに並べ替えるなど、変数列（上述の例では時間列）を反転すれば、上述の変数の進みが正の場合と全く同じ手法で高速演算することができる。
【数２５】
2013024810.tif 000027
【数２６】
2013024810.tif 000028
【数２７】
2013024810.tif 000029
【数２８】
2013024810.tif 000030
【数２９】
2013024810.tif 000031
【００７７】
本実施形態においては、特に、
（１）図４に示すように、区間［ａ，ｔ］内のｃ点（分割点）において、区間［ａ，ｃ］と［ｃ，ｔ］とに分割する、
（２）区間［ａ，ｃ］の積分計算において、入力データｆ（τ）に重み付けした変数ｔに依存しない（変数ｔと独立な）τによる重み付き積分（以下、単に重み付き積分と称す）Ｉｓ，ｋｇ（ｂ）として記憶部に記憶する、
（３）この記憶した重み付き積分Ｉｓ，ｋｇ（ｂ）を用いて積分計算する、
ことを基本思想とするものである。これにより、区間［ａ，ｃ］の積分計算で記憶部に記憶するのはｆ（τ）ではなく、変数ｔに依存しない重み付き積分Ｉｓ，ｋｇ（ｂ）であるから、記憶部容量を小さくできると共に、過去の時刻における積分を再計算する必要がないので、計算時間を短縮でき、このときの計算時間は計算ステップ数ｎに対してｎｌｏｇｎで増加するので、従来のｎ２に比べて劇的に短くできるようにしたものである。なお、図４において、区間［ｃ，ｔ］のｈ刻みの点τは、ｔ側からτ０＝ｔ、τ１＝ｔ−ｈ、τ２＝ｔ−２ｈ、・・・、τＪ０−１＝ｔ−（Ｊ０−１）ｈのように全て異なるが、簡易的に区別せずにτとして関数値を全てｆ（τ）と記載している。また、上述した重み付き積分Ｉｓ，ｋｇ（ｂ）は、式（２５）〜式（２７）の各定義における区間［ａ，ｃ］の積分値Ｉｃｆ（ｔ）を一般化した場合の式（３０）における関数ｇ（τ）の重み付け積分であり、τ＝ｂ−ｕ、ｂ＝ｃｓ−ｊΔｓとして式（３１）で定義される。ここで、ＲＬ定義の場合、α＝−ｑ、ｇ（τ）＝ｆ（τ）であり、Ｇｒｕｎｗａｌｄ−Ｌｅｔｎｉｋｏｖの定義（以下ＧＬ定義）及びＣａｐｕｔｏの定義（以下Ｃ定義）の場合は、α＝ｎｑ−ｑ、ｇ（τ）＝ｆ（ｎｑ）（τ）である。この一般化については後述する。
【数３０】
2013024810.tif 000032
【数３１】
2013024810.tif 000033
【００７８】
さらに、本実施形態においては、区間［ａ，ｃ］をさらに小区間に分割する分割点ｃｓ（ここで、ｓ＝１，２，３・・・Ｓ）を設定したときに、（ｔ−ｃｓ）が十分に大きい条件（分数階微分値を求める変数ｔの点から十分に離れた領域）では、（ｔ−ｃｓ）の大きさ（変数が時間の場合には、現時刻ｔからどのくらい過去の時点か）に応じた範囲毎に既に計算されて記憶部に記憶されている重み付き積分Ｉｓ，ｋｇ（ｂ）を一纏めとした積分値に置き換えて記憶部に記憶し、それを用いて計算しても計算精度を維持することができ、分割点間隔を拡大できるので、重み付き積分Ｉｓ，ｋｇ（ｂ）の数を一層低減することによって計算時間の短縮化と精度の維持が実現できる。
【００７９】
（ＩＩ−１）高速演算化の原理と理論
本発明の分数階微分値の高速演算化の原理と理論について、積分区間が［ａ，ｔ］（ここでａ＜ｔ）において、変数ｔが時間であって、ｔの進む方向がａからｔの方向（正の方向）であるＲＬ定義を例にとり説明する。
【００８０】
式（２５）のＲＬ定義式において、積分する時間区間［ａ，ｔ］を分割点ｃで分割すると、式（２５）は式（３２）に変形することができる。
【数３２】
2013024810.tif 000034
【００８１】
式（３２）において、第１項の積分の上限ｃと第２項の積分の下限ｃとは共に時間とともに変化するが、第２項の積分の下限ｃに対する微分は、ちょうど第１項の積分の上限ｃの微分を相殺するので、式（３２）では、２つの積分のそれぞれ上限及び下限として固定されているとして計算する。そして、区間［ａ，ｃ］の積分に対応する右辺の第１項の積分に注目すると、この第１項の積分は、前述の通り積分区間の上下限が固定されているので、時刻ｔについての微分を積分の中に入れることが可能であり、第１項の積分をＩｃｆ（ｔ）と定義して式（３３）のように変形できる。
【数３３】
2013024810.tif 000035
【００８２】
通常、台形則やシンプソン則による数値積分の古典的な手法では、区間［ａ，ｃ］を短い間隔で分割した区間の積分を用いて数値計算されるので、今、時刻ｂ（ｂ∈（ａ，ｃ］）を上限とした小区間幅Δとして、この小区間の積分ＩΔｆ（ｂ）について考えると、式（３３）は式（３４）の形の積分の和に置き換えられる。ここで、ｂ−Δ≧０である。ただし、時刻ｂは、個々の小区間を識別するラベルのような役割を有しており、個々の小区間を識別できれば、任意の点（時刻）に設定できるものであって、上記のｂの定義はその一例であり、これに限定されない。
【数３４】
2013024810.tif 000036
【００８３】
この小区間の積分の計算について、本発明では、従来の方法では記憶しておいたｆ（τ）を用いて台形則などの近似方法に代えて以下の方法を適用する。即ち、時刻ｂにおける（ｔ−τ）−１−ｑ項について、式（３５）のようにＴａｙｌａｒ展開する。ここで、ｋは展開次数であり、式（３５）は｜τ−ｂ｜＜ｔ−ｂのときに収束する。
【数３５】
2013024810.tif 000037
【００８４】
この式（３５）を式（３４）に代入し、τ＝ｂ−ｕとする。ここでτ＝ｂ−ｕは、図５に示すように、区間［ｂ−Δ，ｂ］内の時刻τに関し、区間上限ｂを基準に定義したものであり、ｂを基準として変数ｕによってτが変化するので、変数τを変数ｕに置き換え、ｕの積分区間を［０，Δ］とすれば、式（３６）が導かれる。こうした置き換えは、区間幅Δの特定区間内における時刻τの積分計算をプログラムする場合においても有用である。
【数３６】
2013024810.tif 000038
【００８５】
こうして、ある特定区間（区間幅Δ）の上限が決まれば、区間幅を積分範囲とする変数ｕを導入した方がプログラム構成でも取り扱いやすくなる。
【００８６】
さらに、ｕ以外の項を積分外に移動して式（３６）を変形し、式（３７）を得る。ここで、ＩΔ，ｋｆ（ｂ）は式（３８）で定義される重み付け積分である。
【数３７】
2013024810.tif 000039
【数３８】
2013024810.tif 000040
【００８７】
式（３７）は、式（３３）のＩｃｆ（ｔ）が、式（３８）で定義されるｕの重み付き積分によって計算できることを示している。そして注目すべきは、式（３８）の重み付き積分が時刻ｔに依存しないので、ｂがｔよりも過去において、式（３８）の重み付き積分ＩΔ，ｋｆ（ｂ）を計算して記憶しておけば、時刻ｔが更新された後にも、その重み付き積分ＩΔ，ｋｆ（ｂ）を再計算することなく式（３３）のＩｃｆ（ｔ）の計算に利用できることであり、これが本発明の重要な第１の原理である。
【００８８】
小区間の幅Δは、式（３７）における級数の収束条件として、式（３９）の条件によって制限される。
【数３９】
2013024810.tif 000041
【００８９】
式（３７）は無限和であるが、数値計算においては近似値として有限の展開次数（ｋ＝Ｋ）まで、式（４０）を用いて計算する。
【数４０】
2013024810.tif 000042
【００９０】
このときの真値と近似値との誤差は、式（４１）に示す通りである。ここで、‖ｆ（τ）‖は、｜ｆ（τ）｜の上限であり、Ｏ（  ）は、Ｌａｎｄａｕのオーダー記号である。
【数４１】
2013024810.tif 000043
【００９１】
式（３８）の小区間の幅Δは、式（３９）の条件を満たせば、任意に設定できる。そして、振動波形の解析問題においては、波形の周期よりも広い間隔で積分すると精度が著しく低下する従来の方式と異なり、この小区間の幅Δは、式（３９）の条件においては振動波長よりも広く設定することさえでき、この点が従来法と全く異なる本発明の重要な第２の原理である。
【００９２】
本発明は、上述した第１の原理と第２の原理に基づいた計算方法の発明であり、これにより高い精度と処理時間の短縮を実現しようとするものである。
【００９３】
（ＩＩ−２）分割点間隔の拡大による更なる高速演算化の理論
次に、時間発展によってｔが増加する状況を考える。時間経過とともにある固定された時点ｂにおける時間間隔Δは、ｔ−ｂに比例して（近似の精度を落とさずに）増加することができる。これが本発明の重要な第３の原理である。
【００９４】
時刻τの単位長さあたりの幅Δの小区間数は、εを比例定数としてΔ＝ε（ｔ−ｂ）とおくと、１／ε（ｔ−ｂ）となる。そして、分割された全ての合計は、（１／ε）ｌｏｇ（ｔ−ａ）で表される時刻τの単位長さ当たりの数を統合して計算される。このように、式（３３）の区間［ａ，ｃ］の積分Ｉｃｆ（ｔ）を計算するための重み付き積分ＩΔ，ｋｆ（ｂ）の総数は、ｌｏｇ（ｔ−ａ）で対数的に増加する。
【００９５】
従って、本発明は、前述した第１の原理及び第２の原理に、さらに第３の原理を適用して、飛躍的に処理時間の短縮化を実現することができる。なお、本原理の導出過程の式（３５）においてＴａｙｌａｒ展開を適用したが、他の展開方法に適用しても同様に、本発明の各原理に帰着する。
【００９６】
（分割点間隔の拡大原理を適用した場合の区間［０，ｃ１］の分数階微分値計算理論）
この第３の原理は、「現時刻ｔから離れるほど（過去に遠ざかるほど）小区間の幅Δを広く設定しても計算精度が維持される」というものであるので、現時刻ｔ、式（３２）の積分下限ａ＝０とした区間［０，ｔ］において、現時刻ｔからの距離によって幅の異なる区間での積分を計算することを考える。
【００９７】
分割点ｃで分割した区間［０，ｃ］区間をさらに中区間［ｃｓ＋１，ｃｓ］（ここでｓ＝１、２、３、・・Ｓ−１）に分割する。ここで、ｃ１＝ｃであり、ｃＳをｃＳ＝０（ｓ＝Ｓ−１のとき）に固定する。なお、ｓ＝０のときには、ｃ０＝ｔとして、区間［ｃ１，ｔ］となり、Δ０＝ｈに相当する。そして、ｔ−Ｊ０ｈがｃ１に相当する。
【００９８】
さらに、中区間［ｃｓ＋１，ｃｓ］内をΔｓ＝２ｓｈの幅に分割した小区間における積分を考える。この幅ΔｓはΔｓ＝２Δｓ−１に相当するが、後述するように２倍に限定されず、２以上の整数倍であればよい。ここで、各小区間のΔｓの個数は１個から複数個まで取り得るので、分割点ｃｓ側から添え字ｊとして０、１、２・・・、Ｊｓ−１と番号付けて考える。
【００９９】
添え字ｓが増えることは現時刻ｔから遠ざかることに対応するので、現時刻ｔから遠い区間［ｃｓ＋１，ｃｓ］（中区間という）内のΔｓの幅も広くなる。
【０１００】
収束条件は、式（３９）において、ΔをΔｓと置き換え、ｂをｃｓと置き換えると、式（４２）となり（ここで、εは精度よい結果を得るための基準値）、その条件下において積分Ｉｃｆ（ｔ）は、式（３８）において、ΔをΔｓと置き、ｂをｃｓ−ｊΔｓと置き換えた式（４５）を用いて式（３７）を区間［ａ，ｃ］において和をとった式（４３）によって計算できる。なお、ここで、ｃｓは中区間の右端点、Δｓは小区間の幅であり、式（４３）の右辺の式（４４）の項は、忘却関数である。また、式（４３）では前記時刻ｂは、ｂ＝ｃｓ−ｊΔｓとおいている。前述のｂ＝ｃｓ−ｊΔｓとおいた場合の時刻ｂは、変数ｔの更新に伴う過去の分割点ｃｓの履歴点（時刻）であって、式（４３）は、ｂ＝ｃｓ−ｊΔｓを上端とする小区間の重み付け積分値Ｉｓ，ｋｆ（ｃｓ−ｊΔｓ）と式（４４）の忘却関数との積を計算して、合算することによって分数階微分値が算出されることを示している。また、式（４２）において、ｓ＝１のとき、式（１）及び式（１０）〜式（１３）となる。また、式（４２）の分母は、積分区間が［ａ，ｔ］の場合にはｔ−ｃｓであり、変数の進みが逆(負の方向)である積分区間［ｔ，ａ］の場合にはｃｓ−ｔである。
なお、式（４３）の計算を実行するプログラムは、各小区間の重み付け積分値の識別は、必ずしもｂ＝ｃｓ−ｊΔｓで行う必要はなく、例えば（ｓ，ｊ）の二次元マトリックスデータとして識別処理する手法の他、公知のコーディング手法を適用してもよい。
【数４２】
2013024810.tif 000044
【数４３】
2013024810.tif 000045
【数４４】
2013024810.tif 000046
【数４５】
2013024810.tif 000047
【０１０１】
式（４３）の中括弧｛｝内が中区間［ｃｓ＋１，ｃｓ］の積分値を与えるので、［０，ｃ１］が複数の中区間［ｃｓ＋１，ｃｓ］で分割されているときには、それらの中区間の積分値を加算すればよいのである。また、式（４３）はｋについて無限和であるため、数値計算においては有限の展開次数（ｋ＝Ｋ）までコンピュータに計算させればよく、近似値として式（４６）によって計算できる。
【数４６】
2013024810.tif 000048
【０１０２】
このときの真値と近似値の誤差は、｜ｆ（ｔ）｜の上限値をＬ、式（４２）の誤差εを用いて式（４７）となる。ここでＯ（  ）は、Ｌａｎｄｏｕのオーダー記号である。
【数４７】
2013024810.tif 000049
【０１０３】
なお、ｓ＝１に限定し、小区間Δ１の個数の上限を設けない場合には、式（４６）は式（４８）となり、区間［０，ｃ１］が複数のΔ１幅からなる小区間で構成される場合には、各小区間の重み付き積分値を合算することで区間［０，ｃ１］の分数階微分値が計算できることになる。
【数４８】
2013024810.tif 000050
【０１０４】
ここで、Ｉ１，ｋｆ（ｃ１−ｊΔ１）は、ｃ１−ｊΔ１＝ｂとして式（４９）で計算できる。
【数４９】
2013024810.tif 000051
【０１０５】
これは、前述した本発明の第１の原理及び第２の原理に基づいた最も基本的な計算方法に相当する。
【０１０６】
以上まとめると、本発明は、前述の第１の原理、第２の原理に基づいて、前述の変数ｔが幅Δ１進む毎に特定条件を満たしながら更新される分割点ｃ１を設定することにより、区間［ａ，ｔ］を区間［ａ，ｃ１］と区間［ｃ１，ｔ］とに分割し、区間［ｃ１，ｔ］は従来法で分数階微分値を算出し、
区間［ａ，ｃ１］においては、変数ｔがΔ１更新される度に、幅Δ１の小区間内のｆ（τ）を用いて、小区間における、変数ｔの更新に依存しない重み付き積分値を計算して記憶し、各小区間について各小区間の重み付け積分値と、各小区間に対応する忘却関数の積を計算したのち、全ての小区間について合算して分数階微分値を算出し、
区間［ａ，ｃ１］と区間［ｃ１，ｔ］の各分数階微分値を合算して、全区間［ａ，ｔ］の分数階微分値を算出することに特徴点を有している。
さらに、前述の小区間の重み付け積分値として、第３の原理に基づき区間［０，ｃ１］をさらに中区間［ｃｓ＋１，ｃｓ］（ここでｓ＝１、２、３、・・Ｓ−１）に分割した場合には、［ｃｓ＋１，ｃｓ］中の小区間Δｓの重み付け積分値は、中区間ｃｓに隣接して連続する中区間［ｃｓ，ｃｓ−１］の小区間Δｓ−１の重み付け積分値を結合して算出したものを適用できることに特徴点を有している。なお、前述の第１から第３の原理やそれらに基づく上記特徴点は、前述したとおり、分数階積分の計算にも適用できるので、分数階微積分値の計算方法の原理及び特徴点といえる。
【０１０７】
（ＩＩ−３）ＲＬ定義による分数階微分式以外への拡張
以上は、ＲＬ定義に沿って説明した本発明の第１の原理〜第３の原理であるが、これらの原理は、分数階微分を積分形式で表記した他の定義を適用した場合においても適用できる。例えば、式（２５）〜式（２７）における、区間［ａ，ｃ］の積分値は、前述の通り、式（３０）として一般表記できる。
【０１０８】
このとき、式（４９）の積分下端ａ＝０とすると、式（４０）は次の式（５０）の通り一般化され、本発明の第１の原理及び第２の原理が、積分項を含む他の分数階微分の定義にも適用できる。
【数５０】
2013024810.tif 000052
【０１０９】
ここで、式（５０）中のＩΔ，ｋｇ（ｂ）は、式（５１）である。
【数５１】
2013024810.tif 000053
【０１１０】
同様に、区間［０，ｃ］の分数階微分値を計算する式（４６）は、次の式（５２）の通り一般化される。
【数５２】
2013024810.tif 000054
【０１１１】
ここで、式（５２）中のＩｓ，ｋｇ（ｃｓ−ｊΔｓ）は、ｂ＝ｃｓ−ｊΔｓとおいた式（３１）である。また、式（５２）中の忘却関数を、ｂ＝ｃｓ−ｊΔｓとおくと式（５３）となる。なお、変数の進みが負の場合には、式（５３）の分子（ｔ−ｂ）α−ｋ−１は、（ｂ−ｔ）α−ｋ−１となる。
【数５３】
2013024810.tif 000055
【０１１２】
さらに、ｓ＝１に限定し、小区間Δ１の個数の上限を設けない場合である式（４８）及び式（４９）は式（５４）及び式（５５）で一般化される。
【数５４】
2013024810.tif 000056
【数５５】
2013024810.tif 000057
【０１１３】
式（５２）〜式（５５）において、ＲＬ定義の場合、α＝−ｑ、ｇ（τ）＝ｆ（τ）とし、ＧＬ定義及びＣ定義の場合は、α＝ｎｑ−ｑ、ｇ（τ）＝ｆ（ｎｑ）（τ）とすれば、各定義における区間［０，ｃ］（又は［０、ｃｓ］）の積分値が得られ、この積分値から分数階微分値が得られ、また、α＞０、ｇ（τ）＝ｆ（τ）とすれば分数階積分の内区間［０，ｃ］における値が計算できる。
【０１１４】
なお、変数の進みが負の場合には、変数の進みが正の方向の場合と逆に、変数が進むと変数は小さくなる（例えば、変数ｔを起点とするとｔ、ｔ−ｈ、ｔ−２ｈ、・・・）。そして、分割点ｃ１はｃ１＝ｔ＋ｍｈ、分割点ｃｓ（ｓ≧２）はｃｓ＝ｃｓ−１＋（Ｌｓ−ｍ１）Δｓ−１であり、分数階微分値を算出は、式（５２）Ｉｓ、ｋｇ（ｃｓ−ｊΔｓ）をＩｓ、ｋｇ（ｃｓ＋ｊΔｓ）に、忘却関数の分子の（ｔ−ｃｓ＋ｊΔｓ）α−ｋ−１を（ｔ＋ｃｓ−ｊΔｓ）α−ｋ−１のように変数列（上述の例では時間列）を反転して計算すればよい。
【０１１５】
（初期値について）
一般に、分数階微分値計算においては、τ≦ａのｆ（τ）の振る舞いを無視しているために、特にａ＝０の場合の事象を取り扱う場合には、τが０以前の履歴がないので、分数階微分方程式の解の解釈に曖昧さが生じる。
【０１１６】
本発明においては、τ≦ａで初期値ｆ（τ）＝０、かつτ＝ａにおいてｆ（ａ）のあらゆる微分の微分係数が０、としている。即ち、ｆ（ａ）＝０、ｆ（１）（ａ）＝０，・・・，ｆ（ｎｑ−１）（ａ）＝０（τ＝ａ、微分階数ｐ＞０）である。この初期値設定の条件であっても本発明の分数階微分値計算の精度は確保されている。
【０１１７】
なお、式（２６）のＧＬ定義には、第２項に初期値項を有するが、式（２６）の第１項の計算における初期値と上述の第２項の初期値項は相殺されるため、式（４９）及び式（５１）には式（２６）の初期値項に起因する項は含まれない。そして、ｆ（ａ）＝０、ｆ（１）（ａ）＝０，・・・，Ｄｎｑ−１ｆ（ａ）＝０のとき、ＧＬ定義及びＣ定義は同じ計算式で計算できる。
【０１１８】
（ＩＩ−４）重み付き積分の計算アルゴリズム
次に、式（４６）の重み付き積分Ｉｓ，ｋｆ（ｂ）の計算方法について、４つの計算ステップに分けてさらに詳細に説明する。
【０１１９】
ここで、現時刻ｔの近傍の時刻τ＝ｂとし、さらに、式（５６）の条件（ここで、εは精度よい結果を得るための上限値）を満たす最小の０より大きいｍの整数倍であるＭを定義する。
【数５６】
2013024810.tif 000058
【０１２０】
式（５６）中のＭは、区間［ｃ１，ｔ］におけるΔ０に等しい時間幅ｈ区間の基準個数である。また、式（５６）中のｍは、幅Δ１の小区間に置き換えられる区間幅ｈの小区間の個数（又は積分ステップの個数）であり、２以上の整数が好ましく、さらにｍ＝２又はＫの値に依存して、ｍは２以上であると本発明の第１の原理〜第３の原理に基づいて計算が最も効率よくできるので特に好ましい。ｍ＝２以外の例として、Ｋ＝２の場合はｍ＝４が計算効率がよい。以下各状態に沿って図６ａ〜図６ｅを用いて説明する。なお、図６ａ〜図６ｅにおいて、区間［ｃ，ｔ］又は区間［ｃ１，ｔ］のｈ刻みの点τは、ｔ側からτ０＝ｔ、τ１＝ｔ−ｈ、τ２＝ｔ−２ｈ、・・・、τＪ０−１＝ｔ−（Ｊ０−１）ｈのように全て異なるが、簡易的に区別せずにτとして関数値を全てｆ（τ）と記載している。
【０１２１】
（計算ステップ１）
計算ステップ１は、現時刻ｔが、ｍｈ／（ｔ−ｍｈ）≧ε（又は、Ｊ０＜Ｍ＋ｍ、ｍは２以上の整数）を満たし、ｃｓ＝ｃ１＝０となる区間の計算ステップである（図６ａ参照）。
【０１２２】
この区間は、本発明の分割点間隔拡大の手法（第１の原理〜第３の原理）を適用すると誤差が大きくなるので、従来法を用いて等間隔の時間幅ｈで数値計算を行う。従って、この計算ステップにおいては、ｆ（τ）の値は全て記憶部に記憶される。
【０１２３】
なお、上述のｍｈはｍｈ＝Δ１に相当する。以下ｍ＝２の条件にて説明するが、ｍは３以上の整数でも構わない。
【０１２４】
（計算ステップ２）
計算ステップ１の状態から時間が経過して、ある現時刻ｔで、初めて２ｈ／（ｔ−２ｈ）＜ε又は、Ｊ０がＪ０＝Ｍ＋２（Ｊ０＝Ｍ＋ｍにおいてｍ＝２としている）を満たした状態となったとき、この時点で、τ＝０、ｈ、２ｈのｆ（τ）を使って式（５７）で積分計算する。これは、式（３８）及び式（４５）において、ｂ＝２ｈ、Δ＝２ｈ（式（４５）においてはΔ１＝２ｈ）に該当する（図６ｂ参照）。
【数５７】
2013024810.tif 000059
【０１２５】
そして、得られた重み付き積分値は、式（５８）のようにラベルを貼って記憶部に記憶させる（プログラミングにおける代入して保存する作業に相当）。ここで記号（：＝）は、右辺を左辺に置換する操作と定義している。
【数５８】
2013024810.tif 000060
【０１２６】
このとき、Ｊ１＝１である。そして、区間［ｃ１，ｔ］の時間幅ｈ（＝Δ０）による分割点数Ｊ０は２個減少してＭ個となる。
【０１２７】
現時刻ｔが（Ｍ＋２）ｈを超えたとき、区間［ｃ１，ｔ］の間隔Δ０のｆ（τ）＝ｆ（ｔ−ｊｈ）（ここでｊ＝０，１・・・Ｊ０）の値は高速演算装置の記憶部に記憶される。これらのデータは、式（３２）の右辺の第２項の計算に使用する。また、これらのデータは、時間の更新によって、順次、式（４６）中の重み付き積分Ｉ１，ｋｆ（ｃ１）の計算にも使用する。
【０１２８】
重み付き積分Ｉｓ，ｋｆ（ｃｓ−ｊΔｓ）の計算の道筋は２つあり、第一の筋道として、最初にｓ＝１における幅Δ１の小区間の重み付き積分Ｉ１，ｋｆ（ｃ１−ｊΔ１）が、区間［ｃ１，ｔ］のｆ（τ）の値を用いて計算される計算ステップ３と、次に、第二の筋道として、ｓ≧２の場合に、第ｓ中区間中の間隔Δｓにおける重み付き積分Ｉｓ，ｋｆ（ｃｓ−ｊΔｓ）が、第（ｓ−１）中区間中のΔｓ−１の積分値Ｉｓ−１，ｋｆ（ｃｓ−１−ｊΔｓ−１）を用いて計算される計算ステップ４とである。
【０１２９】
（計算ステップ３）
計算ステップ３は、時刻ｔが進んで、２ｈ／（ｔ−ｃ１−２ｈ）＜ε（ただし、ｃ１≧０）、又はＪ０がＪ０＝Ｍ＋２（Ｊ０＝Ｍ＋ｍにおいてｍ＝２としている）の条件を満たした状態となったときであって、時刻ｔが進む以前に算出して記憶部に記憶されているｊ＝０，１，２，．．．，Ｊ１−１個の重み付き積分Ｉ１，ｋｆ（ｃ１−ｊΔ１）に加えて、τ＝ｃ１，ｃ１＋ｈ，ｃ１＋２ｈ（ｊ０＝Ｍ＋２，Ｍ＋１，Ｍ又はｊ０＝Ｊ０，Ｊ０−１，Ｊ０−２に対応）におけるｆ（τ）を用いて、ｓ＝１であるΔ１の重み付き積分Ｉ１，ｋｆ（ｃ１＋２ｈ）を式（５９）によって計算し、高速演算装置の記憶部に記憶させる計算ステップである。この式（５９）は、ｂ＝ｃ１＋２ｈ、Δ１＝２ｈとしたときの式（４５）の計算に相当する（図６ｃ参照）。
【数５９】
2013024810.tif 000061
【０１３０】
次に、既に第ｓ＝１中区間に記憶されているＩ１，ｋｆ（ｃ１−ｊΔ１）のラベルを式（６０）の第２式のように貼り替え、式（５９）で得られた値をＩ１，ｋｆ（ｃ１）の値として式（６０）の第３式のように代入する。ただし、式（６０）の第２式及び第３式の右辺のｃ１は、第１式によってｃ１の値を変更（貼り替え）する前のものである。なお、図６ｃにおいては、式（６０）の第２式の右辺の重み付け積分Ｉ１，ｋｆ（ｃ１−ｊΔ１）は、ｂ＝ｃ１−ｊΔ１としてＩ１，ｋｆ（ｂ）と表している。
【数６０】
2013024810.tif 000062
【０１３１】
ここで注目すべきは、重み付き積分及び時刻ｔが進む以前に算出して記憶部に記憶されている重み付き積分は、ラベルのみ貼り替えており再計算しないで、代入記憶のみ実行することである。また、計算ステップ３において時刻ｔがさらに進んで（更新されて）時間経過したときにも、重み付き積分の計算は、更新された［ｃ１，ｃ１＋２ｈ］の区間のみ式（５８）で計算し、時間更新前に計算された重み付き積分は、順次式（６０）に倣ってラベルを貼り替えて記憶部に記憶させていくことである。
【０１３２】
式（６０）におけるＪ１及びＪ０（ｓ＝０，１におけるＪｓに相当）のラベルの貼り替えは、ｓ＝０の幅Δ０の小区間からｓ＝１の幅Δ１の小区間への変換の数の変化として反映されている。
【０１３３】
（計算ステップ４）
この計算ステップ４は、時間を更新して、現時刻ｔが、Δｓ／（ｔ−ｃｓ−Δｓ）＜ε（又はＪｓ−１がＪｓ−１＝（ｍ１−１）Ｍ／ｍ＋ｍ１）の条件を満たした状態となったときであって、ｓ番目の中区間［ｃｓ＋１，ｃｓ］（ここでｓ＝２，３，・・・，Ｓ−１）について、重み付き積分ＩＳ，ｋｆ（ｃｓ＋Δｓ）を、既に計算されて記憶部に記憶された（ｓ−１）番目の中区間の幅Δｓ−１小区間のうち、ｃｓから直近の連続するｍ１個の幅Δｓ−１小区間の重み付き積分値を用いて計算する計算ステップである（図６ｄ及び図６ｅ参照）。ここで、ｍ１は、２以上の整数であって、後述する通り、ｍ１＝２が最も計算効率がよいので、以下ｍ＝２、ｍ１＝２の例で説明するが、ｍ及びｍ１は３以上の整数でも構わない。なお、図６ｄはｓ＝２の場合であって、式（６０）の第２式の右辺の重み付け積分Ｉ１，ｋｆ（ｃ１−ｊΔ１）は、ｂ＝ｃ１−ｊΔ１としてＩ１，ｋｆ（ｂ）と表している。
【０１３４】
重み付き積分Ｉｓ，ｋｆ（ｃｓ＋Δｓ）は、式（６１）で計算する。なお、同式中の２Δｓ−１の２がｍ１に相当する。
【数６１】
2013024810.tif 000063
【０１３５】
式（６１）の右辺は、式（６４）の通り、次の式（６２）及び式（６３）に示すｓ−１番目の中区間を構成するｊ＝Ｊｓ−１−１番目とＪｓ−１−２番目の小区間Δｓ−１の重み付き積分値を合算して求める。
【数６２】
2013024810.tif 000064
【数６３】
2013024810.tif 000065
【数６４】
2013024810.tif 000066
【０１３６】
ＩＳ，ｋｆ（ｃｓ＋Δｓ）が算出された後、式（６５）の通りラベルを貼り替えて代入記憶する。次に、既に算出して記憶されている第ｓ中区間を構成する小区間のΔｓ重み付け積分値Ｉｓ，ｋｆ（ｃｓ−ｊΔｓ）のラベルを式（６５）の第２式のように貼り替え、式（６３）で得られた値をＩｓ，ｋｆ（ｃｓ）の値として式（６５）の第３式のように代入する。ただし、式（６５）の第２式及び第３式の右辺のｃｓは、第１式によってｃｓの値を変更（貼り替え）する前のものである。なお、図６ｅにおいては、式（６５）でｓ−１としたときの第２式の右辺の重み付け積分Ｉｓ−１，ｋｆ（ｃｓ−１−ｊΔｓ−１）は、ｂ＝ｃｓ−１−ｊΔｓ−１としてＩｓ−１，ｋｆ（ｂ）と表している。
【数６５】
2013024810.tif 000067
ここで、Δｓ＝２Δｓ−１である。
【０１３７】
一方、ｍ１が３以上の場合、式（６１）及び式（６４）は、それぞれ式（６６）及び式（６７）で表される。なお、式（６７）において、ｍ１＝２とするとν＝２となり式（６３）となる。
【数６６】
2013024810.tif 000068
【数６７】
2013024810.tif 000069
【０１３８】
なお、式（６７）は分数階微分定義の一般形として式（６８）で表される。
【数６８】
2013024810.tif 000070
【０１３９】
式（６５）のラベル貼り替え作業によって、ラベルＪｓは１個増加し、Ｊｓ−１はＪｓ−１＝Ｍ／２（Ｊｓ−１＝（ｍ１−１）Ｍ／ｍにおいてｍ＝２、ｍ１＝２の場合に相当）となるように２個減少し、また、ｍ１が３以上の場合は、Ｊｓ−１＝（ｍ１−１）Ｍ／２となるようにｍ１個減少する。
【０１４０】
もしｓ＝Ｓの場合には、τ＝０を指定するラベルの添え字は、ＳからＳ＋１に置き換わり、ｃＳ＋１＝０、となる。
【０１４１】
このラベル貼り替え作業は、計算ステップ３と同様に、ラベルのみ貼り替えており再計算しないで、代入記憶のみ実行すればよく、時刻ｔがさらに進んで（更新されて）時間経過したときにも、重み付き積分の計算は、更新された［ｃｓ，ｃｓ＋２Δｓ−１］の区間のみ式（６３）で計算し、時間更新前に計算された重み付き積分は、順次、式（６５）に倣ってラベルを貼り替えて記憶部に記憶させていくのである。
【０１４２】
以上のように、時間の進行に合わせて、区間［ｃ１，ｔ］（ｓ＝０の区間に相当）の時間幅ｈ毎のｆ（τ）＝ｆ（ｔ−ｊｈ）のデータを使って、計算ステップ１から計算ステップ４によって区間［０，ｃ１］中の小区間の重み付き積分値Ｉｓ，ｋｆ（ｃｓ−ｊΔｓ）（ここで、ｊ＝０，１，２，・・・Ｊｓ−１、ｓ＝１，２，３，・・・，Ｓ−１）に置き換えながら記憶部に記憶させ、記憶された重み付き積分値を用いて式（４６）によって、区間［０，ｃ１］の分数階微分値が算出される。
【０１４３】
なお、計算ステップ３及び計算ステップ４において、計算した重み付き積分値は、ラベルを貼り換える手法で時間ｔの進行に対応した時系列データとして識別できるようにしているが、本手法に限定されるものではなく、任意のプログラム手法を適用できることはいうまでもない。
【０１４４】
（全区間［０，ｔ］の分数階微分値の計算）
このように計算ステップ１〜４で計算された重み付き積分値を用いて、区間［０，ｃ１］の分数階微分値は、式（４６）（一般化した表記では式（５２））で計算し、区間［ｃ１，ｔ］の分数階微分値は計算ステップ１のｆ（τ）で従来の方法で計算して、最後に区間［ｃ１，ｔ］と区間［０，ｃ１］の分数階微分値を足し合わせて全区間［０，ｔ］の分数階微分値を得る。
【０１４５】
なお、このアルゴリズムから明らかなように、現時刻ｔにおける分数階微分値の計算において、データｆ（τ）は区間［ｃ１，ｔ］分のみ必要であって、τ≦ｃ１−ｈのｆ（τ）は重み付き積分に置き換わるので記憶部から削除（又は重み付き積分で上書き）して、記憶部容量をさらに節約することができる。
【０１４６】
以上のアルゴリズムにおいて、式（４６）中の右辺中の（ｔ−ｃｓ＋ｊΔｓ）−ｑ−ｋ−１のベキ計算をｓとｊとの二元表に対応するデータ群としてあらかじめ与えておけば（計算データとして計算機に読み込んでおけば）、その分，記憶部容量は必要となるが、計算速度をさらに飛躍的に速くできる。ベキ計算をあらかじめｓとｊとの二元表に対応するデータ群としておくという方法は、単一のプログラムで同一時刻の分数階微分値を多数求める必要のある有限要素解析（ＦＥＭ）や、統計データを求めるため繰り返し計算の必要な場合に、特に有効である。また、ただ一度の計算のみの場合も速度が重要視される計算に対して有効である。さらには、式（５３）の忘却関数自体を先に計算して記憶しておくことも計算速度を向上するために有効である。
【０１４７】
また、別のアルゴリズムとして、図７に示すように、現時刻ｔの更新に合わせて分割点ｃ１＝ｔ−Ｍｈとなるように設定して計算するアルゴリズムとしてもよい。この場合ｓ＝１の最初の区間幅Δ’は、ｈからｍｈの間を変化し、ｍｈのときにΔ’＝Δ１となる。上述のΔ’がｈから（ｍ−１）ｈの間は、区間［ｃ１，ｔ］の分数階微分値計算と同様にｈ刻み毎のｆ（τ）を用いて時間に依存する重み付き積分計算として計算する（図７（ｂ）及び（ｄ））。これは、上述の最初の区間幅Δ１の上端を積分下限として現時刻ｔまでをｈ刻み毎のｆ（τ）を用いて、時間に依存する重み付き積分計算することに等しい。上述のΔ’がｍｈの場合にはΔ１となって、現時刻ｔの更新に依存しない重み付け積分として、上述したアルゴリズムと実質的に同様に計算される（図７（ｃ）及び（ｅ））。なお、図７ａ〜図７ｅは、図６ａ〜図６ｅと同様に、区間［ｃ，ｔ］又は区間［ｃ１，ｔ］のｈ刻みの点τは、ｔ側からτ０＝ｔ、τ１＝ｔ−ｈ、τ２＝ｔ−２ｈ、・・・、τＪ０−１＝ｔ−（Ｊ０−１）ｈのように全て異なるが、簡易的に区別せずにτとして関数値を全てｆ（τ）と記載している。
【０１４８】
また、別のアルゴリズムとして、ｂをｂｓ，ｊ＝ｃｓ−ｊΔｓとおいた式（５２）の重み付き積分値Ｉｓ，ｋｇ（ｃｓ−ｊΔｓ）に代えて、式（５２）中の１／（ｋ！Γ（α−ｋ））を式（５２）に乗じた式（６９）の重み付き積分値Ｉ’ｓ，ｋｇ（ｃｓ−ｊΔｓ）として記憶し、式（７０）で近似計算をすることによって、さらに高速で演算することができる。同様にｓ＝１の場合に対応する式（５４）及び式（５５）はそれぞれ式（７１）及び式（７２）となる。この方法の場合には、記憶した重み付き積分値を纏めて新たな拡大した幅の小区間の重み付き積分値を算出するステップにおけるｍ１＝２の場合の式（６２）は、式（４９）の一般化した形式として式（７３）となり、ｍ１が３以上の場合は、式（７４）から算出できる。なお、式（７４）において、ｍ１＝２とするとν＝２となり式（７３）となる。
【数６９】
2013024810.tif 000071
【数７０】
2013024810.tif 000072
【数７１】
2013024810.tif 000073
【数７２】
2013024810.tif 000074
【数７３】
2013024810.tif 000075
【数７４】
2013024810.tif 000076
【０１４９】
さらに、より計算精度を向上できる方法として、区間［０，ｃ１］の分数階微分値計算においては、ＲＬ定義を適用し、かつ区間［ｃ１，ｔ］の分数階微分値計算をＣ定義とすることがより好ましい。これは、ＲＬ定義においては、その積分核は積分区間を分割することによって現時刻ｔ近傍以外では（ｔ−τ）−ｑ−１とできるため、過去の影響（特に誤差の影響）を小さくでき、一方、Ｃ定義は、現時刻ｔでのｎｑ次までの微分係数が既知であればｔでの分数階微分係数定義式から直接求められるため、両者の特徴を活かして全区間［０，ｔ］の分数階微分値の精度を高めることができるのである。
【０１５０】
この方法における微分定義式は式（７５）となり、ｃが一定の場合にもｔ−ｃが一定の場合にも成り立つ。そして、式（７５）の第１項の計算は、これまでに説明したＲＬ定義での分数階微分値計算のアルゴリズムを適用すればよい。
【０１５１】
なお、図７のアルゴリズムの場合には、現時刻ｔに依存しない重み付き積分の計算領域をＲＬ定義、現時刻ｔに依存する積分計算の領域にＣ定義を適用すればよい。
【数７５】
2013024810.tif 000077
【０１５２】
ちなみに、式（７５）は以下の通り導出される。まず、式（２５）のＲＬ定義式を分点ｃで分割した式（３２）をａ＝０として、式（７６）の通り変形すると、ｃが定数の場合、式（７６）の右辺第１項は式（７７）となる。一方、式（７６）の第２項は式（７８）の通りＣ定義の微分式を関係付けられ、式（７７）及び式（７８）をそれぞれ式（７６）の第１項及び第２項に代入すると式（７５）が得られる。ｃがｔ−ｃで一定の場合にも導出の途中の変形が異なるが、結果として式（７５）となる。
【０１５３】
【数７６】
2013024810.tif 000078
【数７７】
2013024810.tif 000079
【数７８】
2013024810.tif 000080
【０１５４】
（ＩＩ−５）精度について
上記の計算アルゴリズムにおける精度を確保するためには、幅Δｓの小区間は、前述の式（４２）によって制限されるが、さらにこの制限条件について詳細説明する。
【０１５５】
前述した計算ステップ２及び計算ステップ３から理解されるように、第１番目（ｊ＝０番目）の幅Δ１の小区間は、式（４２）から式（７９）と書ける。
【数７９】
2013024810.tif 000081
【０１５６】
第１番目の小区間の幅Δ１はΔ１／（ｔ−ｃ１）＜εを満足する。ｃｓ（ここで、ｓ＝１，２，・・・、Ｓｍａｘ−１）の値は、計算時刻ｔ＝ｎｈ（ここでｎはタイムステップ数）によって変化する。しかしながら、ｓ番目の中区間を構成する幅Δｓの小区間の最少個数は、ｍ＝ｍ１＝２の場合Ｍ／ｍであり、一方、小区間の幅ΔｓはΔｓ＝２ｓｈである。このようにｔ−ｃｓの最小値は、ｍ＝ｍ１＝２として、式（８０）の通り見積もられる。
【数８０】
2013024810.tif 000082
【０１５７】
Δｓ＝２ｓｈと式（８０）との比から明らかなように、式（４２）の左辺Δｓ／（ｔ−ｃｓ）の比は、常にεよりも小さい２／Ｍ以下であり、任意の整数ｓにおいて常に式（４２）を満足するので、分割点の間隔を拡張しても精度が確保できるのである。
【０１５８】
（ＩＩ−６）最小の計算量となるｍ及びｍ１条件
前述の通り、本発明の分数階微分値の計算方法によれば、従来法に比べてデータポイント数と分数階微分値計算に必要な総データ数を減らせるため、計算装置に記憶されるメモリ数も少なくて済むため処理時間を短縮できる。ここまではアルゴリズムを理解し易いように簡略的にｍ＝ｍ１＝２の場合で説明したが、次に記憶部容量と計算量が最小となるｍとｍ１について説明する。
【０１５９】
前述の計算ステップ２及び３において、ｍ≧２の整数ｍ、上述のｍの倍数である整数Ｍにおいて、現時刻ｔ（又はＪ０）がＪ０＝Ｍ＋ｍ，ｍ≧２の条件を満たしたとき、ｃ１＋ｊｈ，ｊ＝０，１，２，・・・，ｍにおける各ｆ（τ）を使ってＩ１，ｋｆ（ｃ１＋ｍｈ）が計算される。
【０１６０】
また、前述の計算ステップ４においては、ｓ−１番目（ここで、ｓ＝２，・・・、Ｓｍａｘ−１）の区間のＪｓ−１がＪｓ−１＝ｍ１＋（ｍ１−１）Ｍ／ｍ，ｍ１≧２のとき、連続するｍ１個分の小区間の重み付き積分Ｉｓ−１，ｋｆ（ｃｓ＋ｊΔｓ−１）（ただし、ｊ＝０，１，２，・・・，ｍ１、ｋ＝０，１，２，・・・，Ｋ）を用いて、ｓ番目の中区間の新たな小区間の重み付き積分値Ｉｓ，ｋｆ（ｃｓ＋Δｓ）として計算される。
【０１６１】
これらの手順では、式（４６）の積分の精度を固定するための式（８１）の条件が課される。この条件は、ｔ−ｃｓが最少となる式（８１）の条件の場合であっても、満足されなければならない。ここで、式（８２）を満たすｓ番目の中区間の小区間の幅Δｓは、式（８３）である。
【数８１】
2013024810.tif 000083
【数８２】
2013024810.tif 000084
【数８３】
2013024810.tif 000085
【０１６２】
時刻ｔ＝ｃｓにおける計算に必要なデータポイント数は、以下の通りである。
（Ａ）ｒ番目の区間の最小幅は、（ｍ１ｒ−ｍ１ｒ−１）Ｍｈ（ｒ＝１，２，・・・，Ｓ−１）であり、小区間Δｒの幅はΔｒ＝ｍｍ１ｒ−１ｈである、
（Ｂ）このように、中区間ｒ番目の幅Δｒの小区間の最小の個数は、（ｍ１−１）Ｍ／ｍであって、ｒには依存しない、
（Ｃ）必要なΔｒの小区間数は、ｒ番目の区間（ｒ＞０）において、（ｍ１−１）Ｍ／ｍ＋ｍ１であり、一方、ｓ＝０番目の区間（区間［ｃ１，ｔ］）では、データの組合せを考慮するとＭ＋ｍ＋１である。
（Ｄ）総データ量をＮｍｅｍとして、展開次数がＫ（Ｋ次）の近似値は、式（８４）で与えられる。
【数８４】
2013024810.tif 000086
【０１６３】
時刻ｔ＝ｃｓに至る分数階微分値の計算の総計算ステップ数ｎは、式（８２）よりｎ＝ｍ１ｓ−１Ｍであることからｓ−１は式（８５）となる。
【数８５】
2013024810.tif 000087
【０１６４】
式（８５）を式（８４）の（ｓ−１）に代入すると、式（８６）が得られる。
【数８６】
2013024810.tif 000088
【０１６５】
式（８６）中、ε０＝ｍ／Ｍは固定された値である。ここで、ｍ１／ｌｏｇ２ｍ１と（ｍ１−１）／ｌｏｇ２ｍ１とは、ｍ１≧２で単調増加する。ｍを固定した場合の総データ量が最小となるｍ１の数は、ｍ１＝２において総データ量が最小となる。ｍ１＝２を式（８６）に代入すると、式（８７）となる。
【数８７】
2013024810.tif 000089
【０１６６】
総データ量が最小となるｍの数は、ｍを実変数として、式（８７）を微分して得た式（８８）の値が０となる条件から、式（８９）に示す値となる。
【数８８】
2013024810.tif 000090
【数８９】
2013024810.tif 000091
【０１６７】
このように、総データ量Ｎｍｅｍは、ｍｍｉｎに最も近い整数ｍのときに最小となる。記憶部に記憶されるデータ量は、次の通り見積もられる。
【０１６８】
今、ｍ＝２ν＋１、またｍ＝２の場合にＭ＝Ｍ２とおく。ここで、νは、分割点ｃ１において、ｍ個の小区間でｆ（τ）を積分した場合、２個の小区間でｆ（τ）を積分するのに比べて減少する中区間の数を意味する。
【０１６９】
ここで、Ｍ／ｍ＝Ｍ２／２と仮定することは精度の点から自然である。このことは、異なるｍ間での積分精度が等しいことを保証している。このことから、Ｍ＝Ｍ２・２νとなる。
【０１７０】
ｍ＝２の場合に、総データ数（Ｎｍｅｍ）２は、ｍ１＝２の式（８４）から式（９０）となる。
【数９０】
2013024810.tif 000092
【０１７１】
一方、任意のｍにおいて分数階微分値計算に必要な総データ数（Ｎｍｅｍ）ｍは、式（９１）となる。
【数９１】
2013024810.tif 000093
【０１７２】
式（９０）と式（９１）との差は、Ｍ／ｍ＝Ｍ２／２を用いて、式（９２）で表される。
【数９２】
2013024810.tif 000094
【０１７３】
式（９２）の差の値は、Ｋ（展開次数）、ｍ及びＭ２／２に依存し、タイムステップには依存しない。
【０１７４】
（総データ数の低減）
精度εＴは、ｍ１の値によらず、式（９３）の条件で拘束される。
【数９３】
2013024810.tif 000095
【０１７５】
展開次数として、式（９３）を満たす最小の整数ＫεＴとすると、式（９４）の右辺を式（９１）に代入すると、式（９５）の通り、与えられた精度を満たす最小のデータ数が求まる。
【数９４】
2013024810.tif 000096
【数９５】
2013024810.tif 000097
【０１７６】
式（８６）は、ｍ１＝２において式（９６）の通りに書き換えられ、総データ数はｎ／Ｍ≫１となるｎに対して対数的に増加する。
【数９６】
2013024810.tif 000098
【０１７７】
これによれば、例えば、Ｋ＝２、ｍ＝２、Ｍ＝４０の条件の場合、分数階微分値計算に必要な総データ数（従来の方法を１として）は、１２５０番目のタイムステップ（ｓ＝５）時点においては約１／３．４、１０２４０番目のタイムステップ（ｓ＝８）時点においては約１／１８、さらに、１６００００番目のタイムステップ（ｓ＝１２）時点においては約１／２００となり、タイムステップが多いほど従来法に比べて総データ量が少なくなるため、記憶部容量が節約されると共に、処理時間が劇的に短縮されるのである。
【０１７８】
以上、分数階微分値の計算原理、アルゴリズムを説明してきたが、分数階積分値の計算への適用について説明すると、分数階積分ａＩαｆ（ｔ）は、式（９７）で定義され、忘却関数とｆ(τ)の積の積分として定義される。式（９７）の積分区間を［ａ，ｃ］と［ｃ，ｔ］とに分割したときの区間［ａ，ｃ］の積分は式（３０）と等しい。従って、上記で述べた分数階微分値の高速計算の原理、理論が適用できるので、分数階積分値の高速計算にも有効である。また変数の進みが負の方向の場合（変数を時間とした場合には未来から現在への流れ）にも、中区間及び小区間の並びを負の向きに並べ替えるなど、変数列（上述の例では時間列）を反転すれば、上述した通常の変数の進み方向（変数の進みが正）の場合の分数階微積分を全く同じ手法を適用できる。変数の進みが逆の方向の場合、変数をｔとすると、分数階積分の定義式は式（９８）である。
【数９７】
2013024810.tif 000099
【数９８】
2013024810.tif 000100
【０１７９】
（ＩＩＩ）第２の実施形態（機器制御システム）
本発明の第２の実施形態として、第１の実施形態における高速演算装置を用いて分数階微分値の高速計算を行う機器制御システムについて説明する。この機器制御システムは、その分数階微分値計算が第１の実施形態の場合と同様に高速で行われるため、高い精度を確保しつつかつ高速な機器制御システムが実現される。
【０１８０】
本実施形態の機器制御システムのハードウェア構成は、第１の実施形態の場合と同様に、図１に示す、コンピュータ１０と、このコンピュータ１０に接続される外部装置１１とから構成されている。
【０１８１】
このようなハードウェア構成の機器制御システムにおいて、ＣＰＵ１０ｂは、作動時は、まず、ＲＡＭ１０ｄ内にプログラム記憶領域、データ記憶領域及びワークエリアを確保し、ＨＤＤ１０ｅからプログラム及びデータを取り込んで、プログラム記憶領域及びデータ記憶領域に格納する。次いで、このプログラム記憶領域に格納されたプログラムに基づいて、図９、及び図２ａ〜図２ｃに示す処理を実行する。ＣＰＵ１０ｂがプログラムを実行することによって、図１０に概略的に示すごとき構成の主要部を備えた機器制御システムが構築される。
【０１８２】
即ち、本実施形態の機器制御システムは、設定データである入力データと制御対象のシステム出力（実測値）とを読み込んで、その差分を計算し、その差分と制御器の分数階微分項を有する伝達関数とから制御器出力値を計算して制御対象に出力し、この制御器出力値と制御対象の伝達関数とからシステム出力を出力して制御対象である被制御機器を作動させると共に、システム出力と次の設定値とから上述の制御フローを繰り返しながら制御するように構成されている。その際、上述の伝達関数は、分数階微分要素の分数階微分式が積分式で定義され、この積分式の分数階微分値が、第１の実施形態における高速演算装置と同様に演算される。このように、制御器の分数階微分項を有する伝達関数に基づいた出力によって機器を制御するものであり、その出力計算に本発明の高速演算装置における分数階微分値計算手法が適用される。以下、フィードバック制御の一つであるＰＩＤ（比例・積分・微分）制御に適用される本実施形態の機器制御システムについて説明する。
【０１８３】
（ＰＩＤ制御）
一般に、ＰＩＤ制御は、設定値と実測値との差である偏差に対して、この偏差に比例して入力値を変化させる動作（比例動作又はＰ動作）と、この偏差の積分に比例して入力値を変化させる動作（積分動作又はＩ動作）と、この偏差の微分に比例して入力値を変化させる動作（微分動作又はＤ動作）とを組み合わせて制御する技術である。そして、制御対象が分数階の動的なシステムの場合には、積分動作と微分動作とをそれぞれＩλ、Ｄμとしたλ階の積分器とμ階の微分器とを有する制御器とすることにより、精度のよい制御が可能となることが知られている。この代表的なものには、Ａ．Ｏｕｓｔａｌｏｕｐ等のＣＲＯＮＥ  Ｃｏｎｔｒｏｌｌｅｒ等があり、盛んに研究されている。
【０１８４】
このようなＰＩＤ制御は、エアコンでは温調制御に、また化学プラントでは、液面制御、圧力制御、又は反応物供給量制御などに適用されているが、積分動作では、周期的変動があり、遅れが生じやすいことが指摘されている。
【０１８５】
ＰＩＤ制御において精度良く制御するには制御時間が重要であり、ＰＩＤの性能向上の一因として、制御時間（分数階微分処理の計算時間）の改善が有効である。
【０１８６】
（ＰＩＤ制御システムの構成）
図８は、本実施形態の機器制御システムとして、ＰＩＤ制御システムの一例をブロック図で示している。
【０１８７】
ＰＩＤ制御システムでは、まず、設定データである入力データＷ（ｓ）と制御対象のシステム出力Ｙ（ｓ）（実測値に相当）とを読み込んでその差分Ｅ（ｓ）を計算し、その差分Ｅ（ｓ）と制御器の伝達関数ＧＣ（ｓ）とから制御器出力値Ｕ（ｓ）を計算して制御対象に出力し、制御対象は制御器出力値Ｕ（ｓ）と制御対象の伝達関数Ｇ（ｓ）とからシステム出力Ｙ（ｓ）を出力して、以降、システム出力Ｙ（ｓ）と次の設定値Ｗ（ｓ）とから上述の制御フローを繰り返しながら制御対象が制御される。なお、ここで、ｓはラプラス演算子であり、上記システムの説明は、時間ｔの関数をラプラス変換したｓ領域で行っている。
【０１８８】
制御対象の伝達関数Ｇ（ｓ）は、実際には各制御対象ごとに異なるが、分数階微分系制御においては、Ｇ（ｓ）は式（１０１）のようなモデルに基づいて、式（１０２）のような分数階微分モデルであると仮定するが、実際の制御においては、式（１０２）の計算は不可能なので、結果として観測されるシステム出力Ｙ（ｓ）を上述の制御フローに適用している。よって、制御器において制御器出力値Ｕ（ｓ）を計算する部分に分数階微分値計算を実施する。
【０１８９】
制御器の伝達関数ＧＣ（ｓ）は、Ｐ動作の比例ゲインＫＰ、Ｉ動作の積分ゲインＫＩ、及びＤ動作の微分ゲインＫＤを用いて表記すると式（９９）となり、時刻ｔの微分方程式として、式（１００）となる。なお以降、実施形態の説明においては、微分又は積分の端点を表す下付き添字は省略する。
【数９９】
2013024810.tif 000101
【数１００】
2013024810.tif 000102
【０１９０】
この式（１００）の分数階微分項の計算には、第１の実施形態における高速演算装置の分数階微分値計算手法を使用する。これにより、高速かつ精度よく分数階微分制御が可能となる。
【数１０１】
2013024810.tif 000103
【数１０２】
2013024810.tif 000104
【０１９１】
本実施形態のＰＩＤ制御システムの一部動作を説明するフロー及びソフトウェアによって構築されるＰＩＤ制御システムの主要部の構成が図９及び図１０に示されている。
【０１９２】
図１０に示すように、ＰＩＤ制御システムは、制御パラメータが入力される制御パラメータ入力部１００と、制御対象において目標とする出力の設定値Ｗ（ｓ）及びその実測値Ｙ（ｓ）が入力されるＷ（ｓ）及びＹ（ｓ）入力部１０１と、これら設定値Ｗ（ｓ）及び実測値Ｙ（ｓ）の差分Ｅ（ｓ）を計算する差分Ｅ(ｓ)演算部１０２と、この差分Ｅ（ｓ）と制御器の伝達関数ＧＣ（ｓ）とから制御器出力値Ｕ（ｓ）を計算する制御器出力値Ｕ（ｓ）演算部１０３と、この制御器出力値Ｕ（ｓ）を出力するＵ（ｓ）出力部１０４と、記憶部１０５とを備えたコンピュータ１０６、並びに制御対象１０７から構成される。差分Ｅ(ｓ)演算部１０２、制御器出力値Ｕ（ｓ）演算部１０３及びＵ（ｓ）出力部１０４は、記憶部１０５に記憶された入力値又は記憶値を用いて計算処理及び出力処理を実行する。ここで、出力されたＵ（ｓ）は、制御対象１０７に入力され、制御対象１０７はその出力Ｕ（ｓ）に応じて作動して応答実測値Ｙ（ｓ）をフィードバックする。
【０１９３】
このうち、制御器出力値Ｕ（ｓ）を計算する制御器出力値Ｕ（ｓ）演算部１０３に、第１の実施形態で述べた高速演算装置が適用されており、本発明の分数階微分値計算アルゴリズムに従って記憶部１０５に記憶された差分の履歴であるデータや重み付き積分値を読み出して、本発明の分数階微分値計算が実行される。
【０１９４】
このＰＩＤ制御システムの動作を図９のフローに従って説明すると、まず、制御パラメータの入力を行って制御条件の設定を行う（ステップＳ２１）。これは、第１の実施形態の高速演算装置における一次パラメータの入力（ステップＳ１）、二次パラメータの計算（ステップＳ２）及び初期条件の設定（ステップＳ３）の処理とＰＩＤ制御システムとしての必要な制御条件の設定処理とを行うものである。次いで、設定値Ｗ（ｓ）と実測値Ｙ（ｓ）とが読み込まれる（ステップＳ２２）。次いで、これら設定値Ｗ（ｓ）と実測値Ｙ（ｓ）との差分Ｅ（ｓ）が算出され、記憶装置に記憶される（ステップＳ２３）。
【０１９５】
その後、記憶装置に記憶された差分Ｅ（ｓ）の履歴データと、制御器の伝達関数ＧＣ（ｓ）との積からなる分数階項を有する式が、第１の実施形態の高速演算装置の場合と同様に計算され、制御器出力値Ｕ（ｓ）が算出される（ステップＳ２４）。算出された制御器出力値Ｕ（ｓ）が制御対象１０７へ出力され（ステップＳ２５）、その後、制御終了の判断が行われる（ステップＳ２６）。終了条件に達しない場合は、ステップＳ２２へ戻り、ステップＳ２２〜Ｓ２６の処理が繰り返して実行される。終了条件に達した場合は（ステップＳ２６においてＹｅｓと判断したとき）にこの処理を終了する。
【０１９６】
本発明の分数階微分値計算方法を適用したＰＩＤ制御は、電気炉やエアコンなどの温調機能を有する機器の温度制御、ハードディスクドライブ（ＨＤＤ）のヘッダの位置を合わせる制御等の各種アクチュエータの駆動制御を含む既存のＰＩＤ制御技術に広く適用できる。本実施形態の機器制御システムは、こうしたＰＩＤ制御の他にも、画像処理としてはフラクタル画像処理の他、データ圧縮方法の一つであるフラクタル圧縮の処理などにも適用することができる。
【０１９７】
（ＩＶ）第３の実施形態（シミュレーションシステム）
本発明の第３の実施形態として、第１の実施形態における高速演算装置を用いて分数階微分値の高速計算を行うシミュレーションシステムについて説明する。このシミュレーションシステムは、その分数階微分値計算が第１の実施形態の場合と同様に高速で行われるため、高い精度を確保しつつかつ高速なシミュレーションが実現される。
【０１９８】
本実施形態のシミュレーションシステムは、分数階微分要素を含む運動方程式でモデル化できる物理現象、例えば振動現象や、ブラウン運動及び拡散現象等、任意の物理現象の解析シミュレーションに応用でき、高い精度を確保しつつ、かつ高速な処理を実現できる。以下、本実施形態のシミュレーションシステムとして、分数階微分でモデル化された粘弾性体の外部応力による形状変形挙動シミュレーションについて説明するが、本発明はこれに限られるものではない。
【０１９９】
本実施形態のシミュレーションシステムのハードウェア構成は、第１の実施形態の場合と同様に、図１に示す、コンピュータ１０と、このコンピュータ１０に接続される外部装置１１とから構成されている。
【０２００】
このようなハードウェア構成のシミュレーションシステムにおいて、ＣＰＵ１０ｂは、作動時は、まず、ＲＡＭ１０ｄ内にプログラム記憶領域、データ記憶領域及びワークエリアを確保し、ＨＤＤ１０ｅからプログラム及びデータを取り込んで、プログラム記憶領域及びデータ記憶領域に格納する。次いで、このプログラム記憶領域に格納されたプログラムに基づいて、図１２、及び図２ａ〜図２ｃに示す処理を実行する。ＣＰＵ１０ｂがプログラムを実行することによって、図１３に概略的に示すごとき構成の主要部を備えたシミュレーションシステムが構築される。
【０２０１】
粘弾性材料の変形解析を行うシミュレーションとしては、公知の種々のシミュレーション法が適用できる。例えば、有限要素解析法（ＦＥＭ）を用いた分数階微分モデルの粘弾性解析が盛んに研究されているため、ＦＥＭ解析手法を例に説明する。その分数階微分要素の計算に本発明の分数階微分値計算方法が適用される。
【０２０２】
（解析モデル）
図１１は粘弾性体の振動シミュレーションを行う場合を説明するためのモデル図である。一例として、粘弾性体に荷重印加される１自由度系の振動子モデルを挙げると、運動方程式は、式（１０３）となる。ただし、式（１０３）において、ｍは質量、ｘ（ｔ）は時間の関数である変位、ｋはバネ定数、ｃは粘弾性係数、Ｐ（ｔ）は時間の関数である外力、Ｄは微分演算子（ｄ／ｄｔ）である。なお以降、実施形態の説明においては、微分又は積分の端点を表す下付き添字は省略する。
【数１０３】
2013024810.tif 000105
【０２０３】
一般に、３次元の等温・等方性の粘弾性体の動的挙動を支配する運動方程式は，有限要素により定式化すると、要素単位では式（１０４）で表される。
【数１０４】
2013024810.tif 000106
【０２０４】
この式（１０４）において、Ｍは質量マトリックスであり、ｄ（ｔ）は節点変位、Ｆ（ｔ）は要素内力ベクトル、Ｐ（ｔ）は等価節点外力ベクトルである。これらＭ、Ｆ（ｔ）、及びＰ（ｔ）は、形状関数をＮ、変位とひずみに関係するマトリックスをＬ、密度をρとすると、体積力ベクトルｆ及び表面力ベクトルＴ（ｔ）を用いて、式（１０５）で定義される。なお、式（１０５）において、Ωは物体の体積Ｖに関する二重積分、Ｓは物体の表面の面積に関する三重積分、上付き添字Ｔは転置を意味する。
【数１０５】
2013024810.tif 000107
【０２０５】
図１１のモデルを３自由度に拡張した場合、例えば、日本機械学会［Ｎｏ．００−６］Ｄｙｎａｍｉｃｓ   ａｎｄ  Ｄｅｓｉｇｎ  Ｃｏｎｆｅｒｅｎｃｅ２０００  ＣＤ−ＲＯＭ論文集［２０００．９．５−８東京］２０７「分数階微分モデル粘弾性体の基礎理論」に開示されているように式（１０５）のＦは、式（１０６）で表され、運動方程式は式（１０７）となる。これは、偏差変形が分数階Ｋｅｌｖｉｎ−Ｖｏｉｇｔ則に，体積変形がＨｏｏｋｅ則に従うモデルである。ただし、式（１０６）及び式（１０７）において、Ｋ’は剛性マトリックス、Ｃ１は粘弾性減衰マトリックス、τは一般化緩和時間、ｑは分数階の微分階数である。また、式（１０７）において、Ｂは変位歪みマトリックスであり、Ｄ０及びＤ１は緩和関数マトリックスであって、Ｄ０は線形の応力歪みマトリックスに対応し、Ｄ１は粘弾性の効果を表すマトリックスに対応する。
【数１０６】
2013024810.tif 000108
【数１０７】
2013024810.tif 000109
【０２０６】
さらに、式（１０７）の左辺に時刻ｔにおける未知の変位ベクトルｄ（ｔ）、右辺に既知のベクトルがくるように変形すると、式（１０８）が得られる。
【数１０８】
2013024810.tif 000110
ここで、Ｋは式（１０９）で表される一般化剛性行列であり、Ｐ’（ｔ）は式（１１０）で表される時刻ｔにおける既知のベクトル（以下、一般化等価節点力ベクトルという）である。
【数１０９】
2013024810.tif 000111
【数１１０】
2013024810.tif 000112
【０２０７】
ここで、Ｓ（ｔｉ−１）は、式（１０８）の左辺中の時刻ｔｉ−１における既知量を纏めた項であり、Ｎｅｗｍａｒｋ−β法の数値積分を適用すると、ｔ０を初期値として式（１１１）で表され、節点変位ベクトルｄ（ｔ）の分数階微分の記憶積分項を含んでいる。ただし、式（１１１）中のα、βは、Ｎｅｗｍａｒｋ−β法のパラメータである。
【数１１１】
2013024810.tif 000113
【０２０８】
ＦＥＭ手法による解析では、現時刻のＰ’（ｔ）と剛性行列Ｋとから現時刻の節点変位ベクトルｄ（ｔ）、節点速度ベクトル、節点加速度ベクトルを計算し、その変位挙動や歪み、応力などの挙動を解析する。
【０２０９】
（シミュレーションシステムの構成）
図１３に示すように、シミュレーションシステムは、粘弾性体の形状メッシュデータ、材料データ及び外力データをシステムの記憶部１３６や入力ファイルから入力するデータ入力部１３０と、入力されたデータを用いて質量マトリックスＭ、粘弾性減衰マトリックスＣ１、剛性行列Ｋ’を作成して記憶部１３６に記憶させるマトリックス演算部１３１と、節点変位ベクトルｄ（ｔ）の履歴を用いて分数階微分値を計算し、形状メッシュデータに対応した一般化等価節点外力ベクトルＰ’（ｔ）を算出して記憶部１３６に記憶させる演算部１３２と、一般化剛性行列Ｋを算出する一般化剛性行列Ｋ演算部１３３と、記憶部１３６に記憶されている一般化剛性行列Ｋと一般化等価節点外力ベクトルＰ’（ｔ）とを用いて剛性方程式を解いて節点変位ベクトルｄ（ｔ）を算出して記憶部１３６に記憶させる変位ベクトル演算部１３４と、入力結果や演算結果をディスプレイ等に出力する出力部１３５と、記憶部１３６とを備えたコンピュータから構成される。データ入力部１３０、マトリックス演算部１３１、一般化等価節点外力ベクトルＰ’（ｔ）演算部１３２、一般化剛性行列Ｋ演算部１３３、変位ベクトル演算部１３４及び出力部１３５と、記憶部１３６とは、入力データや各種算出データを記憶、及び／又は読み出し作業がなされる構成となっている。さらに、節点変位ベクトルｄ（ｔ）、節点変位速度ベクトル、節点変位加速度ベクトルを用いて算出される物性値の挙動を計算する場合には、物性予測値演算部を含んでもよい。なお、上述した一般化等価節点外力ベクトルＰ’（ｔ）演算部１３２において、節点変位ベクトルｄ（ｔ）の分数階微分値Ｄｑｄ（ｔ）を求める際に、第１の実施形態で述べた高速演算装置による分数階微分値計算が実行される。
【０２１０】
（解析処理フロー）
以下、本実施形態における有限要素法（ＦＥＭ）による粘弾性材料の外力による形状変形挙動のシミュレーションを図１２のフローチャートを用いて説明する。線形の分数階微分モデルにおいては、節点変位ベクトルｄ（ｔ）の分数階微分値Ｄｑｄ（ｔ）の計算に、第１の実施形態で述べた高速演算装置による分数階微分値計算手法が適用される。
【０２１１】
まず、粘弾性体の形状データ、材料データ、等価節点外力ベクトルＰ（ｔ）データ、初期条件、一般化剛性行列Ｋの拘束条件（幾何学的拘束条件）、一般化等価節点ベクトルＰ’（ｔ）の拘束条件（力学的拘束条件）、一般化緩和時間τ等の解析条件を読み込み、記憶部１３６に記憶させる処理等を含む解析初期条件の設定を行う（ステップＳ３１）。
【０２１２】
次いで、時間繰り返しインデックスｉをｉ＝０に初期設定する（ステップＳ３２）。
【０２１３】
次いで、記憶部１３６に記憶されている粘弾性体の形状データ、材料データ、等価節点外力ベクトルＰ（ｔ）データを読み出して、質量マトリックスＭ、粘弾性減衰マトリックスＣ１、剛性行列Ｋ’を式（１０５）中のＭの定義式及び式（１０７）中のＣ１及び式（１０９）中のＫ’の定義式によって算出して記憶部１３６に記憶させる（ステップＳ３３）。
【０２１４】
次いで、時間繰り返しインデックスｉをｉ＝ｉ＋１として、時間を更新した後（ステップＳ３４）、節点変位ベクトルｄ（ｔ）分数階微分値Ｄｑｄ（ｔ）を算出して記憶部１３６に記憶させる（ステップＳ３５）。
【０２１５】
その後、記憶された節点変位ベクトルｄ（ｔ）の分数階微分値Ｄｑｄ（ｔ）を記憶部１３６から読み出してＳ（ｔｉ−１）を算出し、記憶部１３６に記憶させる（ステップＳ３６）。
【０２１６】
次いで、記憶部１３６から等価節点外力ベクトルＰ（ｔ）とＳ（ｔｉ−１）とを読み出して、式（１１０）から一般化等価節点外力ベクトルＰ’（ｔ）を算出して記憶部に記憶させる（ステップＳ３７）。
【０２１７】
次いで、記憶部１３６から質量マトリックスＭ、粘弾性減衰マトリックスＣ１、剛性行列Ｋ及び一般化緩和時間τを読み出して、式（１０９）により一般化剛性行列Ｋ’を算出し、記憶部１３６に記憶させる（ステップＳ３８）。
【０２１８】
次いで、一般化剛性行列Ｋ’に幾何学的拘束の処理と、一般化等価節点ベクトルＰ’（ｔ）に力学的拘束の処理を行った後、記憶部１３６に記憶させる（ステップＳ３９）。
【０２１９】
次いで、記憶部１３６に記憶されている幾何学的拘束処理された一般化剛性行列Ｋ’と、力学的拘束処理された一般化等価節点外力ベクトルＰ’（ｔ）とを読み出し、式（１０８）によって、時刻ｔの節点変位ベクトルｄ（ｔｉ）、節点速度ベクトル、節点加速度ベクトルを算出し、記憶部１３６に記憶させる（ステップＳ４０）。
【０２２０】
その後、計算終了かどうかの判断を行い（ステップＳ４１）、時間繰り返しインデックスｉがｉ≧ｎではない場合（Ｎｏの場合）は、ステップＳ３４へ戻って時間を更新し、ステップＳ３４〜Ｓ４１の同様のループで、各時刻毎の節点変位ベクトルｄ（ｔｉ）、節点速度ベクトル、節点加速度ベクトルを算出し、各時刻の履歴データとして記憶部１３６に記憶させる。この変位ベクトルｄ（ｔ）の各時刻の履歴データは、変位ベクトルの分数階微分値Ｄｑｄ（ｔｉ−１）において、式（３１）の分数階微分定義式におけるｆ（τ）に相当し、この履歴データについて本発明の分数階微分値の計算が適用される。これにより分数階微分値Ｄｑｄ（ｔ）の分数階微分値が精度よく高速で得られ、その結果、シミュレーションの精度を保ちつつ従来に比較して短時間で解析を行うことが可能となる。
【０２２１】
ステップＳ４１において、時間繰り返しインデックスｉがｉ≧ｎと判断した場合（Ｙｅｓの場合）は、算出されて記憶部１３６に記憶されている各時刻のデータ、即ち変位ベクトルｄ（ｔ）、変位速度ベクトル、変位加速度ベクトルが、記憶部１３６から読み出され、公知のデータ加工手段によって、グラフィック表示などの形式で出力される（ステップＳ４２）。また、算出して記憶部に記憶されている節点変位ベクトルｄ（ｔ）、節点速度ベクトル、節点加速度ベクトルを用いて、節点歪みベクトルや要素内応力ベクトルを算出することもできる。
【０２２２】
なお、図１２においては、分数階微分値Ｄｑｄ（ｔｉ−１）の算出とＰ’（ｔｉ−１）の算出とを分けて実行するように記載されているが、プログラムにおいては、Ｐ’（ｔｉ−１）を算出する計算ステップの中でサブルーチンとして分数階微分値Ｄｑｄ（ｔｉ−１）の算出を実行させてもよい。
【実施例】
【０２２３】
本発明の高速演算装置の処理時間効果を検証するために、式（１１２）で与えられる分数階微分モデルの伝達関数と、式（１１３)で与えられる分数階微分モデルの制御システムとからなる閉ループ系ＰＤμ制御（α＜μ＜β）における出力の計算時間について、本発明を適用した実施例１と、従来法を適用した比較例１及び２とを比較した。
【数１１２】
2013024810.tif 000114
【数１１３】
2013024810.tif 000115
【０２２４】
（実施例１）
前述のＰＤμ制御は、出力をｙ（ｔ）、入力をｗ（ｔ）として、式（１１４）の微分方程式でモデル化できる。式（１１４）の微分方程式は、式（１１５）を用いて式（１１４）を（１−β）階微分して式（１１６）へと変形し、式（１１６）を離散化した式（１１７）（添え字「０」はｎ−１ステップ目の値を示す）から現時刻ｔの出力ｙ（ｔ）の値として出力ｙｎを数値計算する。
【数１１４】
2013024810.tif 000116
【数１１５】
2013024810.tif 000117
【数１１６】
2013024810.tif 000118
【数１１７】
2013024810.tif 000119
【０２２５】
具体的には、前述の離散化により、式（１１７）右辺の各分数階微分項は式（１１８）及び式（１１９）の分数階積分として計算できるので、式（１１７）中のＡのうち式（１１８）及び式（１１９）中のｊ＝０からｎ−１までの和の部分をＡ０ｙｎ（実際には出力ｙｎにあたる項は含まれない）、ｊ＝ｎ項目の部分をＡｎｙｎとすると、式（１１７）は式（１２０）となり、式（１２０）の両辺を［１＋ｈＡｎ／２］で除して出力ｙｎを数値計算した。式（１１６）の係数ａ０及びａ１は、ａ０＝ａ１＝１とし、ゲイン値Ｋ、ＴｄはそれぞれＫ＝２０．５、Ｔｄ＝３．７３４３とした。また、μ＝１．１５である。
【数１１８】
2013024810.tif 000120
【数１１９】
2013024810.tif 000121
【数１２０】
2013024810.tif 000122
【０２２６】
分数階微分値の計算は、式（９４）の条件を満足するパラメータとして、ｍ＝ｍ１＝２、Ｍ＝４０、Ｋ＝２、ε０＝１／２０とし、ＲＬ定義に基づく式（４６）によって、区間［０，ｃ１］の分数階微分値を算出し、区間［ｃ１，ｔ］の分数階微分値は、式（３２）の第２項のＲＬ微分の値を台形則によって求め、区間［０，ｃ１］の分数階微分値と区間［ｃ１，ｔ］の分数階微分値とを合算して全区間［０，ｔ］の分数階微分値を求めるプログラムとした。なお、区間［ｃ１，ｔ］と［０，ｃ１］の積分中のベキ式の値は、あらかじめ計算したデータを記憶部に記憶しておいて計算を実行した。【０２２７】
（比較例１）
比較例１として、実施例１においてＲＬ定義による分数階微分値の計算を「台形則」で計算を行った。「台形則」とは、台形公式はニュートン・コーツの公式と呼ばれる数値積分の公式のひとつであって、被積分関数を一次関数で近似し、台形の面積の公式に帰着させて積分の近似値を求める方法であり、式（１２１）の要領で計算するものである。このＲＬ微分中のベキ式の値は、あらかじめ計算したデータを記憶部に記憶しておき、これを読み出して計算を実行した。
【数１２１】
2013024810.tif 000123
【０２２８】
（比較例２）
比較例２として、実施例１において従来のＧＬ定義による分数階微分値の計算を行った。「ＧＬ定義」による計算法は、式（１２２）を用いて計算する。ＧＬ定義による計算において、式（１２２）の係数である式（１２３）は、あらかじめ計算して記憶部に記憶させておき、これを読み出して式（１２２）の計算を実行している。
【数１２２】
2013024810.tif 000124
【数１２３】
2013024810.tif 000125
【０２２９】
実施例１並びに比較例１及び２について、計算時間（ＣＰＵ時間）を計算ステップ数に対してプロットした結果を図１４に示す。図１４より、本発明の高速計算装置では、計算の数千ステップ以上における台形則、ＧＬ定義の場合に比べてＣＰＵ時間が低く抑えられていることが分かる。このことから、本発明の高速演算装置を用いた制御は、従来法よりも極めて短時間で出力計算できるので、追従性の高い制御が実現できることが分かる。計算の精度は、一次の精度でステップ数が等しければ、本発明の実施例１、比較例１である台形則、及び比較例２であるＧＬ法の精度はほぼ同等である。【０２３０】
本発明の高速演算装置の処理時間効果を検証するために、シミュレーションシステムについて、式（１２４）の分数階微分を含む粘弾性体の１自由度系の振動子モデルの分数階微分値計算の処理時間について、本発明を適用した実施例２と、従来法を適用した比較例３及び４とを比較した。
【数１２４】
2013024810.tif 000126
ただし、式（１２４）において、ｍｓは質量、ｘ（ｔ）は変位、ｔは時間、ｃは粘弾性係数、ｑは分数階数（０より大きい実数）、ｋｓはバネ定数、ｆ（ｔ）は質量に与えられる外力を表す。【０２３１】
（実施例２）
まず、時刻ｔ＞０において、式（１２４）にｆ（ｔ）＝ｓｉｎ（ｔ）を入力した場合の応答を計算する。その他の式（１０９）のパラメータは、ｍ＝ｃ＝ｋ＝１，ｑ＝０．５とおいた。本発明の高速計算装置を使用して、式（１２４）の分数階微分値の計算を行った。この実施例２の分数階微分値の計算は、式（９４）の条件を満足するパラメータとして、ｍ＝ｍ１＝２、Ｍ＝４０、Ｋ＝２、ε０＝１／２０とし、ＲＬ定義に基づく式（４６）によって、区間［０，ｃ１］の分数階微分値を算出し、区間［ｃ１，ｔ］の分数階微分値は、式（３２）の第２項のＲＬ微分の値を台形則によって求め、区間［０，ｃ１］の分数階微分値と区間［ｃ１，ｔ］の分数階微分値とを合算して全区間［０，ｔ］の分数階微分値を求めるプログラムとした。なお、区間［ｃ１，ｔ］と［０，ｃ１］の積分中のベキ式の値は、あらかじめ計算したデータを記憶部に記憶しておいて計算を実行した。【０２３２】
（比較例３）
実施例２において、分数階微分値の計算を比較例１と同じ台形則で計算した。【０２３３】
（比較例４）
実施例２において、分数階微分値の計算を比較例２と同じＧＬ法で計算した。【０２３４】
実施例２並びに比較例３及び４について、計算時間（ＣＰＵ時間）を計算ステップ数に対してプロットした結果を図１５に示す。図１５より、本発明の高速計算方法では、計算の数千ステップ以上における台形則、ＧＬ定義の場合に比べてＣＰＵ時間が低く抑えられていることが分かる。各ステップ数までに必要なメモリ数は、台形則ではステップ数に等しい。ＧＬ法ではｘ（ｔ）のための記憶部と係数のための記憶部が必要で、ステップ数の２倍である。従って、台形法とＧＬ法はＮ２アルゴリズムの結果となっている。一方、本発明によれば、必要なメモリ数は、Ｍ＝４０、ｍ＝ｍ１＝２、Ｋ＝２の場合、ステップ数＝４０×２ｎに対して４３＋３×２２ｎであり、ＮｌｏｇＮアルゴリズムである。計算の精度は、一次の精度でステップ数が等しければ、本発明の実施例２、比較例３である台形則、及び比較例４であるＧＬ法の精度はほぼ同等である。【０２３５】
このように本発明により高速計算アルゴリズムを使用することによって、これまで実現が難しいとされてきた分数階微分方程式の数値計算が飛躍的に短時間で計算でき、分数階微分要素を含む現象の解析などを実用的な時間で行うことができる。【０２３６】
以上述べた実施形態は全て本発明を例示的に示すものであって限定的に示すものではなく、本発明は他の種々の変形態様及び変更態様で実施することができる。従って本発明の範囲は特許請求の範囲及びその均等範囲によってのみ規定されるものである。
【産業上の利用可能性】
【０２３７】
本発明の分数階微分値の計算を行う高速演算装置は、分数階微分要素を含む様々な現象の経時的な挙動や特性解析（シミュレーション）や信号制御、画像処理、経済推移予測などに幅広く適用でき、この装置を用いれば、精度よく短時間で処理できるので、商品開発の効率アップや、生産性、製品の品質向上、さらに経済変動の予測情報の提供などにも利用することができる。
【符号の説明】
【０２３８】
１０  コンピュータ
１０ａ  バス
１０ｂ  ＣＰＵ
１０ｃ  ＲＯＭ
１０ｄ  ＲＡＭ
１０ｅ  ＨＤＤ
１０ｆ  音声処理部
１０ｇ  画像処理部
１０ｈ  ディスク駆動装置
１０ｉ  入出力インタフェース
１０ｊ  スピーカ
１０ｋ  ディスプレイ
１０ｌ  キーボード
１０ｍ  マウス
１０ｎ  プリンタ
１１  外部装置
３０  入力部
３１  二次パラメータ演算部
３２  第１の分数階微分値演算部
３３  重み付き積分演算部
３４  第２の分数階微分値演算部
３５  全分数階微分値算出部
３６、１３５  出力部
３７、１０５、１３６  記憶部
１００  制御パラメータ入力部
１０１  Ｗ（ｓ）及びＹ（ｓ）入力部
１０２  差分Ｅ(ｓ)演算部
１０３  制御システム出力値Ｕ（ｓ）演算部
１０４  Ｕ（ｓ）出力部
１０６  コンピュータ
１０７  制御対象
１３０  データ入力部
１３１  マトリックス演算部
１３２  一般化等価節点外力ベクトルＰ’（ｔ）演算部
１３３  一般剛性行列Ｋ演算部
１３４  変位ベクトル演算部
【図１】
2013024810.tif 000140
【図２ａ】
2013024810.tif 000141
【図２ｂ】
2013024810.tif 000142
【図２ｃ】
2013024810.tif 000143
【図３】
2013024810.tif 000144
【図４】
2013024810.tif 000145
【図５】
2013024810.tif 000146
【図６ａ】
2013024810.tif 000147
【図６ｂ】
2013024810.tif 000148
【図６ｃ】
2013024810.tif 000149
【図６ｄ】
2013024810.tif 000150
【図６ｅ】
2013024810.tif 000151
【図７】
2013024810.tif 000152
【図８】
2013024810.tif 000153
【図９】
2013024810.tif 000154
【図１０】
2013024810.tif 000155
【図１１】
2013024810.tif 000156
【図１２】
2013024810.tif 000157
【図１３】
2013024810.tif 000158
【図１４】
2013024810.tif 000159
【図１５】
2013024810.tif 000160
【国際調査報告】
2013024810.tif 000161
2013024810.tif 000162