(57)【要約】【課題】  磁場のローレンツ力による荷電粒子の速度変化を解析する。【解決手段】  磁場のローレンツ力による荷電粒子の速度変化は、初期速度に対して所定時間後の速度を得るために、差分形式の運動方程式を２段階に分けて演算するボリスの方法が知られている。これに対し、運動解析装置内には、ボリスの方法を３段階以上に分けて解析できるよう拡張した一般式に基づいて各段階に応じた計算式を格納したモジュールを用意し、オペレータの指示または要求精度等に基づく自動設定によって解析段数を設定した上で、当該段数に応じたモジュールを使用して解析を実施する。  こうすることにより、ボリスの方法に比較して、高精度、高速度で解析を行うことが可能となる。【選択図】  図４
2019003457.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータによって、磁場のローレンツ力により荷電粒子のΔｔ時間後の速度ｕ＋を求める運動方程式をｎ（ｎは３以上の整数）段階のステップに分けて数値計算する運動解析方法であって、
前記コンピュータが実行する工程として、
（ａ）初期の速度ベクトルｕ―、磁場ベクトルＢ、前記荷電粒子の電荷ｑ、および前記荷電粒子の質量ｍを入力し、前記コンピュータのメモリの所定領域に格納する工程と、
（ｂ）前記ｎ段階のうちｌ番目（ｌは１以上ｎ以下の整数）の速度ベクトルｕτｌを次の数式により計算し、順次ｕτｌをメモリの所定領域に格納する工程と、
【数１】
2019003457.tif 000017
ここで
【数２】
2019003457.tif 000018
（上式においてｃは光速を表す）
（ｃ）前記工程（ｂ）をｎ回繰り返して実行して得られる速度ベクトルｕτｎを速度ｕ＋として出力する工程とを備える運動解析方法。
【請求項２】
請求項１記載の運動解析方法であって、さらに
予め複数のｎに対して、前記工程（ｂ）で用いる数式をそれぞれ用意し、前記コンピュータのメモリ領域に格納しておく工程と、
運動解析に用いるｎを設定する設定工程とを有し、
前記工程（ｂ）は、前記設定されたｎに対応する数式を用いて前記ｎ段階の計算を行う運動解析方法。
【請求項３】
請求項１記載の運動解析方法であって、
運動解析に用いるｎを設定する設定工程を有し、
前記工程（ｂ）は、ｌ番目（ｌは１以上ｎ以下の整数）の計算ごとに、順次β（ｌ，ｎ）をメモリの所定領域に格納することにより、前記数式を設定して前記ｎ段階の計算を行う運動解析方法。
【請求項４】
請求項２または３記載の運動解析方法であって、
予め精度とｎとの関係、および計算時間とｎとの関係を記憶したデータベースを用意する工程と、
前記工程（ｂ）の計算に先立って、解析の要求精度および許容計算時間を入力する工程とを備え、
前記設定工程は、前記要求精度および前記許容計算時間に応じて前記データベースを参照して前記運動解析に用いるｎを設定する工程である運動解析方法。
【請求項５】
コンピュータによって、磁場のローレンツ力により荷電粒子のΔｔ時間後の速度ｕ＋を求める運動方程式をｎ（ｎは３以上の整数）段階のステップに分けて数値計算し運動解析を行わせるためのコンピュータプログラムであって、
（ａ）初期の速度ベクトルｕ―、磁場ベクトルＢ、前記荷電粒子の電荷ｑ、および前記荷電粒子の質量ｍを入力し、前記コンピュータのメモリの所定領域に格納する機能と、
（ｂ）前記ｎ段階のうちｌ番目（ｌは１以上ｎ以下の整数）の速度ベクトルｕτｌを次の数式により計算し、順次ｕτｌをメモリの所定領域に格納する機能と、
【数３】
2019003457.tif 000019
ここで
【数４】
2019003457.tif 000020
（上式においてｃは光速を表す）
（ｃ）前記工程（ｂ）をｎ回繰り返して実行して得られる速度ベクトルｕτｎを速度ｕ＋として出力する機能とをコンピュータに実現させるコンピュータプログラム。
【請求項６】
磁場のローレンツ力により荷電粒子のΔｔ時間後の速度ｕ＋を求める運動方程式をｎ（ｎは３以上の整数）段階のステップに分けて数値計算する運動解析装置であって、
（ａ）初期の速度ベクトルｕ―、磁場ベクトルＢ、前記荷電粒子の電荷ｑ、および前記荷電粒子の質量ｍを入力し、記憶する解析条件入力部と、
（ｂ）前記ｎ段階のうちｌ番目（ｌは１以上ｎ以下の整数）の速度ベクトルｕτｌを次の数式により計算し、順次ｕτｌを記憶する解析実行部と、
【数５】
2019003457.tif 000021
ここで
【数６】
2019003457.tif 000022
（上式においてｃは光速を表す）
（ｃ）前記解析実行部での計算をｎ回繰り返して実行して得られる速度ベクトルｕτｎを速度ｕ＋として出力する結果出力部とを備える運動解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、磁場のローレンツ力による荷電粒子の速度変化を解析する運動解析方法に関する。
【背景技術】
【０００２】
電磁界環境の中で荷電粒子またはプラズマ粒子の運動を解析するための数値解析の手法として、ボリスの方法が知られている。以下、ボリスの方法について、簡単に説明する。
【０００３】
まず電磁界環境下での粒子の運動方程式は次式で与えられる。
【数１】
2019003457.tif 000003
ここで
ｘは荷電粒子の位置ベクトル；
ｔは時間；
ｖは荷電粒子の速度；
ｍは荷電粒子の質量；
ｑは荷電粒子の電荷；
Ｅは電界ベクトル；
Ｂは磁界ベクトル
を意味する。
【０００４】
上記運動方程式を時間差分の形で記載すると、次式となる。
【数２】
2019003457.tif 000004
【０００５】
ここでｕ、Ｂの定義は次式の通りである。また、ｕ，Ｅ，Ｂの右上に示したt等の符号は、時刻ｔにおける速度（ｕ）、電界（Ｅ）、磁界（Ｂ）であることを意味する。ｃは光速を意味する。
【数３】
2019003457.tif 000005
【０００６】
さらに、新たな変数ｕ−、ｕ＋を次式の通り定義する。
【数４】
2019003457.tif 000006
このとき、時間差分の運動方程式は、次式のように書き換えられる。
【数５】
2019003457.tif 000007
【０００７】
ボリスの方法とは、上記運動方程式の解法であり、初期速度ｕ−からΔｔ時間後の速度ｕ＋を算出する際に、途中のｕτ１を経て次式で表される２段階の計算によって解を得るという方法である。
【数６】
2019003457.tif 000008
【０００８】
上式において、速度ｕと磁界Ｂの外積は、速度ｕを回転させる作用を表している。図１はボリスの方法における演算を模式的に示した説明図である。
第１段階では、初期速度ｕ−と磁界Ｂとの外積と、初期速度ｕ−との和によって速度ｕτ１が得られる。次に、第２段階では、初期速度ｕτ１と磁界Ｂとの外積と、初期速度ｕ−との和によって速度ｕ＋が得られる。初期速度ｕ−とｕ＋の大きさは等しい。このようにボリスの方法では、外積を２段階にわたって求めることになり、速度ｕ＋を得るために初期速度ｕ−を角度Ａ１１、Ａ１２の２段階で回転させることを意味している。角度Ａ１２は、Ａ１１の２倍となる。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】H.Qin, S.Zhang, J.Liu, Y.Sun, W.M.Tang, Why is Boris algorithm so good? Phys. Plasmas 20(2013)084503
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかし、ボリスの方法は、精度が悪く計算の所要時間も大きいという課題があった。図１に示した通り、速度ｕ＋からさらに所定時間経過した後の速度ｕ２を求めようとする場合、ボリスの方法では、速度ｕ＋をさらに、角度Ａ２１、Ａ２２のように２段階で回転させる計算を行う必要がある。即ち、長時間経過後の速度を計算する場合、ボリスの方法では、２段階の計算を何度も繰り返す必要が生じ、解析時間が長くなることがある。
一方、２段階の計算における回転角Ａ１１、Ａ１２等を大きくすることにより、演算の所要時間を抑える方法も考えられるが、かかる方法をとる場合には、精度が低くなるという課題が生じる。
本発明は、これらの課題に鑑み、コンピュータを利用して、荷電粒子の運動解析を行う際に、計算速度および精度の向上を図ることを目的とする。
【課題を解決するための手段】
【００１１】
本発明は、
コンピュータによって、磁場のローレンツ力により荷電粒子のΔｔ時間後の速度ｕ＋を求める運動方程式をｎ（ｎは３以上の整数）段階のステップに分けて数値計算する運動解析方法であって、
前記コンピュータが実行する工程として、
（ａ）初期の速度ベクトルｕ―、磁場ベクトルＢ、前記荷電粒子の電荷ｑ、および前記荷電粒子の質量ｍを入力し、前記コンピュータのメモリの所定領域に格納する工程と、
（ｂ）前記ｎ段階のうちｌ番目（ｌは１以上ｎ以下の整数）の速度ベクトルｕτｌを次の数式により計算し、順次ｕτｌをメモリの所定領域に格納する工程と、
【数７】
2019003457.tif 000009
ここで
【数８】
2019003457.tif 000010
（ｃ）前記工程（ｂ）をｎ回繰り返して実行して得られる速度ベクトルｕτｎを速度ｕ＋として出力する工程とを備える運動解析方法と構成することができる。
【００１２】
本願の発明者は、２段階で計算するボリスの方法を３段階、４段階というように拡張する研究を重ねた結果、ｎ段階（ｎは３以上の整数）で計算する際の一般式を見いだした。それが上記式である。図１に示したボリスの方法においては、角度Ａ１１、Ａ１２と２段階で速度を回転させる解析を行っていたのに対し、本発明の式を用いれば、３段階以上の任意の段階に分けて計算することができる。この結果、１段ごとの回転角を抑制できるため、解析の精度が向上する。また、ボリスの方法では、長時間経過後の状態を解析するためには、２段階の計算を繰り返す必要があり、全体としての計算量が増加していたが、本発明の方法ではｎ段階の計算で解を得ることができるため、計算量を抑制することもできる。
【００１３】
計算量について、一例として３段階でｕ＋を求める場合を考える。図２は本発明の方法における３段階での計算を模式的に示した説明図である。３段階の解は、上述の一般式にｎ＝３を当てはめることによって得られる。この場合の計算も、図２に示すように、ｕ−を角度Ａ３１，Ａ３２，Ａ３３だけ回転した結果としてｕ＋が求められる。角度Ａ３２は角度Ａ３１の２倍、角度Ａ３３は角度Ａ３２の２倍、即ち角度Ａ３１の４倍となっている。このことは、図２のｕ＋の解を得るため、ボリスの方法では、２段階を２回繰り返す必要があることを意味し、合計４段階分の計算を要することを意味している。
同様に、４段階の場合を考える。図３は本発明の方法における４段階での計算を模式的に示した図である。回転角度の関係だけを簡単に示した。図中の最外の円に示すように、４段階で計算する場合には、４回の回転を行うことになるが、それぞれ従前の回転角度の２倍ずつの角度で回転されることになる。従って、同様の解をボリスの方法で得ようとすれば、図示する通り８回（＝２３回）の計算が必要となることになる。
このように、一般的に、ｎ段階の計算と同様の結果をボリスの方法で得る場合には、２ｎ−１回の計算を要することになる。従って、本発明の方法によれば、ボリスの方法に比較して、大幅に計算量を抑制することができることになる。
【００１４】
本発明の運動解析方法においては、さらに
予め複数のｎに対して、前記工程（ｂ）で用いる数式をそれぞれ用意し、前記コンピュータのメモリ領域に格納しておく工程と、
運動解析に用いるｎを設定する設定工程とを有し、
前記工程（ｂ）は、前記設定されたｎに対応する数式を用いて前記ｎ段階の計算を行うものとしてもよい。
【００１５】
例えば、３段階で計算するための計算式、４段階で計算するための計算式などを予め用意しておく態様である。こうすることにより、３段階での解析を実行する際には、予め用意された３段階のための計算式を用い、４段階での解析を実行する際には、予め用意された４段階のための計算式を用いればよく、容易に各段数に応じた解析を実現することができる。
計算式を予め用意しておく方法としては、例えば、予め計算に必要な係数、または式全体をデータベースとして格納しておくようにしてもよいし、各段数に応じた解析を実現するためのモジュールを格納しておくものとしてもよい。
計算の段数ｎは、例えば、オペレータが指示するものとしてもよいし、後述するように要求精度等に応じて自動的に設定されるようにしてもよい。
本発明は、このように複数の段階に応じた計算式を用意する態様に限らず、当初から３段階のみ、４段階のみというように特定の段階に応じた計算を行う態様をとるものとしても良い。
【００１６】
本発明の運動解析方法においては、
運動解析に用いるｎを設定する設定工程を有し、
前記工程（ｂ）は、ｌ番目（ｌは１以上ｎ以下の整数）の計算ごとに、順次β（ｌ，ｎ）をメモリの所定領域に格納することにより、前記数式を設定して前記ｎ段階の計算を行うものとしてもよい。
指定された段数の計算を行う際に、逐次、必要な係数等を求めながら解析を実行する態様である。かかる態様によれば、予め段数に応じた計算式を用意する必要がないという利点がある。
【００１７】
本発明においては、
予め精度とｎとの関係、および計算時間とｎとの関係を記憶したデータベースを用意する工程と、
前記工程（ｂ）の計算に先立って、解析の要求精度および許容計算時間を入力する工程とを備え、
前記設定工程は、前記要求精度および前記許容計算時間に応じて前記データベースを参照して前記運動解析に用いるｎを設定する工程であるものとしてもよい。
【００１８】
本発明の解析においては、段数ｎを増やすことによって精度が向上する一方、計算の所要時間が増大することがわかっている。従って、要求精度に応じて段数ｎの下限値を求めることができ、許容計算時間に基づいて段数ｎの上限値を求めることができる。解析に用いる段数ｎはこうして得られた下限値、上限値の間で設定すればよい。候補の段数ｎが複数存在する場合には、オペレータに選択させるようにしてもよいし、候補の段数ｎの中で最小値または最大値を選ぶというように予め選択方法を設定しておいてもよい。
【００１９】
本発明は、上述した運動解析方法としての態様に限らず、かかる運動解析を実現するための解析装置、コンピュータによってかかる解析を実現するためのコンピュータプログラム、およびかかるプログラムを記録したコンピュータ読み取り可能な記録媒体として構成することもできる。ここでコンピュータ読み取り可能な記録媒体としては、フラッシュメモリ、ハードディスク、光ディスクなど種々の媒体を用いることができる。 
【図面の簡単な説明】
【００２０】
【図１】ボリスの方法における演算を模式的に示した説明図である。
【図２】本発明の方法における３段階での計算を模式的に示した説明図である。
【図３】本発明の方法における４段階での計算を模式的に示した説明図である。
【図４】運動解析装置の構成を示した説明図である。
【図５】段数と誤差、所要時間の関係を示すグラフである。
【図６】解析処理のフローチャートである。
【図７】解析結果を示す説明図である。
【発明を実施するための形態】
【００２１】
図４は、運動解析装置の構成を示した説明図である。スーパーコンピュータに、図示する各機能を実現するコンピュータプログラムをインストールすることによって構成した。運動解析装置は、複数のコンピュータ、サーバ等をネットワークで接続して構成してもよい。
【００２２】
以下、図中の各機能ブロックについて説明する。
解析条件入力部１１は、運動解析に要する種々の条件を入力する。例えば、荷電粒子の位置ベクトル、時間、荷電粒子の初期速度、荷電粒子の質量、荷電粒子の電荷、電界、磁界などの条件が挙げられる。また、運動解析を何段階の演算で行うかの指示を入力するようにしてもよい。
【００２３】
解析段数設定部１２は、運動解析を実行する段数ｎを設定する。オペレータが段数ｎを指示している場合には、それを利用すればよい。オペレータが指示していない場合には、精度・所要時間データベース１３を参照して決定することになる。
【００２４】
精度・所要時間データベース１３は、段数と精度および計算所要時間の関係を格納するデータベースである。図５を用いて、その内容について説明する。図５は、段数と誤差、所要時間の関係を示すグラフである。
図５（ａ）は、横軸に運動解析後の速度の回転角度、縦軸に１回の時間ステップ当たりの回転角度の限界値を表したグラフである。解析の段数２〜８までのそれぞれについての解析結果が示されている。例えば、段数２の結果によれば、初期速度から２πだけ回転した速度が得られるような運動解析を実現しようとしても、１回の時間ステップ当たりの回転角はπ以下が限界であることが分かる。このことは、つまり、２πだけ回転した速度が得られるような運動解析を２段階の計算で行う場合には、２回以上の時間ステップが必要となることを意味しており、その分、長時間を要することを意味する。図５（ａ）に示される通り、１回の時間ステップで解析できる限界の角度は段数が増えるに従って増大し、概ね６段以上で計算を行えば、１回の時間ステップでも、２πだけ回転した速度が得られるような運動解析を実現することができることになる。
【００２５】
一方、図５（ｂ）は、横軸に１回の時間ステップでの回転角の限界値を示し、縦軸に誤差を示したグラフである。解析の段数２〜８までのそれぞれについての解析結果が示されている。この結果によれば、１回の時間ステップでの回転角が同一の場合、解析の段数が増大するほど誤差が小さくなることが分かる。また、所定の誤差を実現しようとした場合、解析の段数が小さくなるほど１回の時間ステップ当たりの回転角を小さくする必要があることが分かる。
【００２６】
以上の通り、図５のグラフは、解析の段数ｎと、所要時間（図５（ａ））および誤差（図５（ｂ））の関係を表していることになる。精度・所要時間データベース１３は、このような解析結果を予めデータベースとして格納しているのである。データの格納形式は、図５（ａ）、図５（ｂ）に示した形式に限られず、段階数ｎと精度および所要時間との関係を直接または間接に表すものであればよい。
【００２７】
図４に示した解析段数設定部１２は、精度・所要時間データベース１３を参照して、オペレータによって入力された要求精度および許容所要時間を満たす段数ｎを設定することになる。具体的には、段数ｎと精度との関係に基づき、要求精度を満たす段数の下限値を求め、段数ｎと所要時間との関係に基づき、許容所要時間を満たす段数の上限値を求め、両者の間で解析に用いる段数ｎを選択すればよい。複数の候補段数が存在する場合には、オペレータに選択させてもよいし、最小値・最大値・中間値など予め定めた規則に従って選択するようにしてもよい。
【００２８】
解析実行部２０は、荷電粒子の運動解析を実行する。図示する通り、本実施例では、解析実行部２０には、解析の段階数に応じたモジュールを用意した。２段階解析モジュール２２は２段階で解析を実行するため以下の計算式を予め組み込んだモジュールである。
【数９】
2019003457.tif 000011
【００２９】
３段階解析モジュール２３は、３段階で解析を実行するため以下の計算式を予め組み込んだモジュールである。
【数１０】
2019003457.tif 000012
【００３０】
４段階解析モジュール２４は、４段階で解析を実行するため以下の計算式を予め組み込んだモジュールである。
【数１１】
2019003457.tif 000013
【００３１】
５段階解析モジュール２５は、５段階で解析を実行するため以下の計算式を予め組み込んだモジュールである。
【数１２】
2019003457.tif 000014
【００３２】
６段階解析モジュール２６は、６段階で解析を実行するため以下の計算式を予め組み込んだモジュールである。
【数１３】
2019003457.tif 000015
【００３３】
７段階以上の解析を行う際には、一般式から段階数に応じて計算式を求めて解析を実行する。計算式設定モジュール２８は、以下に示す一般式を組み込んでいる。
【数１４】
2019003457.tif 000016
【００３４】
上式から分かる通りｌ番目の係数βは、ｌ−１番目、ｌ−２番目等の係数を用いて求められる形式となっている。従って、計算式設定モジュール２８は、ｌ番目の計算で用いた係数βをそれぞれメモリ内に格納し、次の段階の係数の算出に利用する。これらの係数のうち、予め算出可能なものについては、係数テーブル２９に予め計算した結果を格納しておくようにしてもよい。こうすることで、計算式の設定に要する時間を短縮することができる。
【００３５】
また、いずれの段階の計算においても速度の変化を２〜６段階、またはそれ以上の段階に分けて求めていくものであるから、各段階で得られた速度の演算結果などは、逐次、メモリに格納し、それを次の段階以降の計算で利用する。
【００３６】
解析実行部２０は、解析の段数に応じて、２段階解析モジュール２２〜６段階解析モジュール２６、または計算式設定モジュール２８を選択し、解析を実行する。
結果出力部３０は、以上で得られた解析結果を出力する。
【００３７】
本実施例では、２〜６段階の解析については、それぞれ個別のモジュールを用意したが、これらを省略し、全ての段階に対して計算式設定モジュール２８のみで対応するようにしてもよい。また、逆に、計算式設定モジュール２８を省略し、７段階以上についても個別のモジュールを用意するようにしてもよい。
【００３８】
図６は、解析処理のフローチャートである。図４に示した運動解析装置が実行する処理である。
解析処理を開始すると、運動解析装置は、解析条件を入力する（ステップＳ１０）。入力される解析条件は、解析条件入力部１１（図４参照）において説明した通りである。
次に、運動解析装置は、解析段数を設定する（ステップＳ１１）。解析段数がオペレータによって指定されている場合は、その段数を用いればよい。また、指定されていない場合には、精度・所要時間データベース１３を参照して、要求精度および許容所要時間を満たす範囲内で段数を設定すればよい。
なお、オペレータによって解析段数が指定されている場合においても、精度・所要時間データベース１３を参照して、得られる精度および所要時間をオペレータに提示するようにしてもよい。オペレータは、この結果を見て、指定した段数を要求精度等に適したものに修正することが可能となる。
【００３９】
運動解析装置は、設定された段数が６以下の場合（ステップＳ１２）、段数に応じたモジュールを用いて解析を実行する（ステップＳ１３）。即ち、２段階が設定されたときには２段階解析モジュール２２を使用し、３段階が設定されたときには３段階解析モジュール２３を使用するという具合である。
一方、設定された段数が６より大きい場合には（ステップＳ１２）、計算式設定モジュール２８を用いて、段数に応じた計算式を設定しながら解析を実行する（ステップＳ１４）。
【００４０】
運動解析装置は、以上で得られた結果を出力し（ステップＳ１５）、解析処理を終了する。
【００４１】
図７は、解析結果を示す説明図である。横軸に速度ベクトルの回転角度を示し、縦軸に速度の縦方向成分の変化を示した。解析は、図中の丸印、即ち、π／２の回転角度ごとに行った。解析結果の精度がよければ、縦方向の速度はｃｏｓ曲線（図中の破線）に一致するはずである。
図７（ａ）〜図７（ｆ）では、それぞれ３段階〜８段階における解析結果を実線で示している。図示する通り、３段階では、実線と破線のずれが確認され、計算誤差が生じていることが分かるが、段数が増すにつれて両者のずれは小さくなり、６段階以上ではほとんど誤差の影響が見られなくなっていることが確認される。
このように本実施例の解析手法によれば、段数を増すことで非常に良い精度で解析を実現することが分かる。
【００４２】
以上で説明した実施例の運動解析装置、運動解析方法によれば、電磁界の影響中での荷電粒子の運動を精度よく高速に計算することができる。
なお、本発明は、本実施例に示した態様に関わらず、その趣旨を変更しない範囲で種々の態様により実現可能である。
【産業上の利用可能性】
【００４３】
本発明は、磁場のローレンツ力による荷電粒子の速度変化を解析するために利用することができる。
【符号の説明】
【００４４】
１１    ：解析条件入力部
１２    ：解析段数設定部
１３    ：精度・所要時間データベース
２０    ：解析実行部
２２    ：２段階解析モジュール
２３    ：３段階解析モジュール
２４    ：４段階解析モジュール
２５    ：５段階解析モジュール
２６    ：６段階解析モジュール
２８    ：計算式設定モジュール
２９    ：係数テーブル
３０    ：結果出力部
【図１】
2019003457.tif 000023
【図２】
2019003457.tif 000024
【図３】
2019003457.tif 000025
【図４】
2019003457.tif 000026
【図５】
2019003457.tif 000027
【図６】
2019003457.tif 000028
【図７】
2019003457.tif 000029