UNISYS TECHNOLOGY REVIEW

雑感

第 129 号，SEP. 2016

ENIAC で実行された数値計算について
長

1.

はじめ

島毅に

この度，ENIAC 完成 70 周年を記念して米国にて発刊された ENIAC in Action の邦訳「ENIAC

─

現代計算技術のフロンティア［1］

─」（以下，本書）が刊行されたが，筆者は，ご縁を得て，その翻訳作業の一部をお手伝いさせて頂く過程で，ENIAC 開発当時の様々な物語に初めて接することができた．
世界で最初に実用化された電子計算機 ENIAC は，米国陸軍の資金援助の下，1943 年にペンシルバニア大学のムーアスクールで開発開始され，同校の J. プレスバー・エッカー
ト Jr. とジョン・W・モークリーが主導して完成した．稼働後，陸軍の弾道研究所に移設され，1955 年に退役するまで，全米各地から訪れる研究者達の持ち込む様々な研究課題を解析するための共用設備として長く活躍した．何もかもが手探り状態の中から，電子計算機のアーキテクチャを構成すること，それを実装するのに必要な電子部品を選定すること，汎用性を実現するためのモジュラー構造の採用，計算を指示するプロセスを，ループや分岐，サブルーチンなどを含むプログラミングという技術として確立すること，解析すべき問題の定式化，定式化された問題の数値解法の考察，…計算機にまつわる，ありとあらゆる概念が，ほぼ同時に並行してこの時代に生まれ，成熟して，今日もほとんどその本質を変えることなく利用されていることは驚異以外の何物でもない．本書においては，その全ての過程が，それを担った多くの天才達の人生とともに，余すところ無く描かれていて，大変興味深い．多くの皆様に是非一読をお勧めしたいと思う．
本稿では，本書で言及されている ENIAC で実行された数値計算の代表的なもののうち，本文だけでは分かりにくいと感じたいくつかの計算テクノロジーについて，筆者なりに補足説明とその後の発展についての所感を述べてみたい．
2. モンテカルロ・シミュレーションについて
本書の第 8 章・9 章に，ENIAC 稼働初期に大々的に実施された，核分裂現象のモンテカルロ法を用いたシミュレーションをめぐる様々なエピソードが紹介されているが，数値計算法としての説明が不十分なので，補足説明をさせて頂くことにする．
まず，シミュレーションという言葉を簡単に説明しておこう．計算機上に実世界の事象（システム）の数学的モデルを作り，そのモデルを使って事象の挙動を模倣することをシミュレーション（数値実験）と呼ぶ．システムの特性が確定的ではなく，確率分布に従ってバラつく要素を含む場合があるが，このようなシステムを確率的システムと呼ぶ．
モンテカルロ・シミュレーションとは，確率的システムの挙動を，効率的に把握することができる非常に有用な数値解析法である．この手法では，システムの数学的モデルの確（115）69

70（116）

率的な要素に対して，乱数を利用して無作為に選び出した有限個の標本値を当てはめ，その挙動から実システムの挙動を推測する．名称のモンテカルロは，カジノで有名なモナコ公国のモンテカルロに由来するが，カジノでサイコロを振る動作と，乱数を用いて無作為に標本値を選び出す操作の連想から命名されたものと思われる．
特性値にバラつきがある要素を含む数学的モデルを正攻法で解析する場合にまず思いつくのは，バラつきのある要素それぞれについて，そのバラつき範囲をいくつかの区間に分割し，それぞれの区間の代表値を取り出し，それぞれの代表値同士の組合せについて総当りで解析して，系の性能値の最良値，最悪値，平均値，特性曲線の包絡線などを求める方法であろう．ただし，この方法では，バラつく要素の数が多数になると解析すべきケース数が級数的に大きくなり，現実的な計算時間では処理できなくなる欠点を持つ．一方，統計学において，母集団から無作為に抽出したサンプル集団の分布は，サンプル数を適切に選べば，母集団の分布を十分な精度で近似することができるという性質が知られている．
具体例として，世論調査などで何百万人，何千万人もの人達の意見をアンケートで調査する場合に，無作為に抽出した 100 人〜 1000 人程度にアンケートを実施して，十分な精度で世論を推定できるという事実が挙げられる．この性質を利用して，性能にバラつきのある構成要素を多数含むシステムの挙動解析の効率化を図ったのが，モンテカルロ法である．この方法においては，各々のバラつきのある要素について，乱数を利用して無作為なサンプル値を選び出し，要素の値を確定したサンプルシステムの挙動を解析する計算を，サンプルを選びなおしながら複数回繰り返し，総当りアプローチより桁違いに少ない回数の解析で，十分な信頼性を確保しつつ，システムの性能の推定値を得ることが可能になる．
3.

モンテカルロ法による核分裂現象の解析
次に第 8 章に記述されている，モンテカルロ法を用いた核分裂現象の解析について，分かりやすく解説してみたい．
ENIAC が本格的に稼働し始めたのは，第二次世界大戦直後で，当時米国は既に原子爆弾を完成させていたものの，その性能は不完全で貴重なウラン資源の 99％は未使用のまま浪費されるレベルであった．世界は冷戦下の軍事力拡張競争の真っ只中で，米国は次世代核兵器，すなわち水素爆弾の完成に躍起になっていた．水素爆弾は超高温高圧の環境下で，プラズマ化した重水素あるいは三重水素の原子核が融合してヘリウム原子核を生成する際に放出される莫大なエネルギーを，破壊力に転化した兵器である．その当時，核融合の起爆に必要な超高温高圧環境を実現することができる唯一の手段は，核分裂反応であった．そこで考案された水素爆弾の構造は，中心に小規模の原子爆弾をおき，その周りを重水素と三重水素を含んだ殻構造で囲み，まず原子爆弾を爆発させ，その爆発によって得られる超高温高圧の環境を利用して，核融合反応を二次的に起爆するというものである．
当時は米国のウラン濃縮能力は小規模で，戦略的高濃度ウランの蓄積量も微少であった．従って，少量のウランを効率よく核分裂させるための最適な核爆弾の設計技術を，実験や試作による貴重なウラン資源の浪費をできるだけ抑えた研究方法で確立するというのが，国家的急務であった．

雑感

ENIAC で実行された数値計算について（117）71

原子爆弾はウランやプルトニウムのような不安定な放射性物質の小球のまわりを，中性子をはじき返すタンパーと呼ばれる物質と，中性子を吸収する物質の混合層で同心球状に覆った形状を基本としている．核分裂が連鎖的に拡大するように，放射性物質の球殻構造を最適化するのが設計の命題となる．
この問題をモンテカルロ法で効率よく解析したのが，ENIAC による 4 期にわたるモンテカルロ・ランと呼ばれる解析プロジェクトである．
この解析では，以下のような中性子の確率モデルが使用された．
①極座標系で定義され，その初期値を，乱数を用いて与えられた，位置座標をもつ．
②大きさと角度を分布関数で確率的に規定された，速度をもつ．
③中性子は，所定の時間後に図 1 に示されているような 4 種類の状態のいずれかに遷移するが，それぞれに遷移する確率が定義されている．ただし，この遷移確率は固定的なものではなく，確率的に変化する中性子の速度の関数として変動するが，傾向として速度が大きくなれば，中性子がウランの原子核に衝突する確率が大きくなるので，「派生」の状態の発生確率が大きくなり，その分それ以外の状態の発生確率は減少する．

図 1 中性子の遷移状態「離脱」は，中性子の速度が大きくて解析対象領域から飛び出した状態であり，これ以降この状態の中性子は解析対象から除外される．「吸収」は他の粒子に衝突して，吸収されてしまう状態であり，やはりこれ以降の解析対象からは除外される．
「散乱」は障害物に衝突して位置が変わるが，解析領域内に存在し続ける状態で，引き続き解析対象となる．
「派生」はまさに核分裂に対応する状態であり，ウランの原子核に衝突してそれを分裂させ，派生的に新しい四つの中性子を生み出す状態で，自身はそこで解析対象から外れるが，派生された四つの中性子が新たに解析対象に加えられる．
解析対象とするすべての中性子モデルについて，乱数を用いて，それぞれの次の瞬間の位置，及び遷移状態を決定する．遷移した四つの状態のうち，「離脱」と「吸収」となった中性子モデルを対象から除外し，「派生」された新しい中性子に対応する新しい中性子モデルを解析対象に加えて，また次の瞬間について同じ操作を繰り返す（図 2）．所定の時間に到達するまでこの過程を繰り返して，一定時間経過後に解析対象として存在する中性子モデルの総数を計測する．その数が初期の 100 に比べて十分に大きくなっていれば，核分裂が臨界状態を超え，爆発的連鎖反応を起こしていると判定できる．この過程を繰り72（118）

返すことにより，本質的には膨大な数の中性子の状態変化の様子を，ごく少数の統計的サンプルの挙動を解析するだけで，高い精度で近似的に求めることができ，貴重なウランやプルトニウム資源を実験で浪費することなく最適設計値を求めることに多大な貢献をした．

図2

核分裂過程のモンテカルロ解析（概念図）

余談になるが，ENIAC には大容量の記憶領域もディスク等の外部記憶装置もなかったので，核分裂モデルの多数の中性子に関するデータを計算機内に保持することができなかった．そこで，モンテカルロ・シミュレーションを一つの時間刻み分計算するたびに，各中性子の状態を 1 枚ずつカードにパンチ出力し，次のステップではそのカードを分別して，「離脱」と「吸収」の分を取り除き，残りを読み込みながら，次のステップの状態を計算してパンチ出力することによって，計算を進めていったという記述が本書にある．つまり，一つの時間刻みごとに大量のカードが打出され，人手でそれを取り出し，ソーターで分別し，再入力分をカード読取装置にセットし，計算を再開するというサイクルを，延々繰り返したのである．時間，手間，資源消費の全ての面において，今日では考えられないような光景であるが，それでも当時は，手計算では全く歯が立たない計算を画期的な速度で解いた快挙，という事なのである．
4.

モンテカルロ法における，乱数生成
上述の様にモンテカルロ法においては，統計的なパラメタの値を与えられた分布内か

ら，乱数を用いて無作為に一つ選び出す作業が多用される．この乱数が，統計学的に純粋な乱数であればあるほど，モンテカルロ法の根拠になっている「サンプル集団の統計的性質と母集団の統計的性質の相似性」が具現できるため，使用する乱数列の品質が大変重要になってくる．
初期には，乱数は乱数表から拾った値をカードにパンチして，プログラムの外から与え雑感

ENIAC で実行された数値計算について（119）73られていたが，解析の効率向上のために，プログラム内部で自動的に生成することが，初めて試みられた．プログラムのロジックである限り，本来は作為的であり再現性を持ってしまうことはやむを得ないが，その中でもできるだけ無作為に近い，互いに関連性の無い数の列を，与えられた範囲で均等にバラつくように生成する擬似乱数生成法が数値計算技術として確立され，用いられるようになっている．
ENIAC におけるモンテカルロ解析では，本解析の中心人物である，ジョン・フォン・ノイマンによって考案された以下のような演算が用いられたと，本書に紹介されている．
すなわち，N 桁の大きな数を最初の乱数として出力し，それを二乗して 2N 桁の数を作り，中央の N 桁の数字を取り出して次の乱数出力とする．以降このプロセスを繰り返す，という方法である．この方法は「平方採中（Mid Square）法」と呼ばれて，今でも数値計算法の教科書などに解説されている．
また，Rn+1 = mod [N1 ´ Rn + N2, NN]（ただし N1，N2 は十分に大きな互いに異なる素数，NN は計算機で表現可能な最大の整数，関数 mod は第一引数を第二引数で割った余りを出力する） という線形合同法と呼ばれる方法も良く用いられてきた．どちらの方法も，前と同じ値が一度出現すると，同じ系列を繰り返して発生してしまうという周期性を持つため，乱数としての完全性にかける．近年では，その周期性を使用する整数表現ビット数の範囲で極大化した「線形帰還シフトレジスタ法」（図 3）や，更に 2 進数表現された生成乱数値において，いくつかの位置の桁同士の値を入れ替えて周期をより長大化させる，「メルセンヌ・ツイスタ法」などが考案され，生成する擬似乱数の品質が向上している．

図 3 線形帰還シフトレジスタ法

それでも，上記全ての方法で生成される擬似乱数は過去の履歴が既知であれば未来の生成値が予測可能であるため，暗号論的に安全ではないとされる．暗号論的な安全性を確保するために，従来法で生成された擬似乱数を暗号論的ハッシュ関数に通すという，「暗号論的擬似乱数生成器」という方法も存在する．
いまや，公開されたブラックボックス関数として乱数生成機能が手軽に利用できるため，その生成方法や品質について，陽に意識することは少なくなってきたと思われるが，かつて乱数生成サブルーチンから自分でプログラミングして来た古い時代のプログラマの生き残りとして，「たかが乱数，されど乱数」という思いから，久しぶりに語ってみた．
5. 弾道計算と数値積分法
本書の第 1 章に，大砲の射程から砲身の仰角を求めるために軍が使用していた射表につ74（120）いて，かなり詳しい説明が載っている．この射表を計算するためには，弾道の運動特性を記述する常微分方程式を解く必要があったが，その式のなかに砲弾速度・温度・風速・気圧などに関して非線形の関数で表される空気抵抗項が存在するため，解析的には求解できないので，専用の手回し式の機械式計算機を用いて，微小時間幅毎の位置を計算していた．
軍の弾道研究所に集められた 100 人もの専任要員が，昼夜 2 交代で機械式計算機を回し続けても，軍の需要を満足できず計算の積み残しは増大する一方だった．電子計算機開発を目指していたムーアスクールの研究者はこの事実に着目し，射表計算の高速化を実現するための電子式計算装置の提案を軍に提出し，逼迫する研究費を軍事予算から調達することに成功し，開発物語が始まる訳である．
弾丸の軌道は，以下のような常微分方程式で表現される．
dx=udt

x t =0 = 0

du= -Dx (u, v,T , P ,W )dtu t =0 = V0 cos( 0 )

dy=v
dt

y t =0 = 0dv= - g - Dy (u, v , T , P , W )dtv t =0 = V0 sin( 0 ),0：空気抵抗， ：温度， ：気圧，：風力，：初速度，θ0：打出し仰角， ：重力加速度

これを，初期条件から時間について逐次積分していくことで，軌道の時間変化が計算できる．計算機による数値積分計算は，数値計算の代表的な技術のひとつであり，今日まで精度と効率を改善するために様々な技法が開発されている．
数値積分法は，その解法の特性によって，いくつかに分類できる．代表的な分類軸として，以下のものがある．
①陽解法と陰解法②単段法と多段法③固定時間刻み幅法と可変時間刻み幅法
これらの軸は互いに独立で，組合せ可能である．例えば，有名なルンゲクッタ・ギル法
は，単段型・時間刻み幅固定・陽解法と分類できる．また，予測子修正子型の代表解法であるミルン法は多段型・時間刻み可変・陽解法である．本書 38 ページ図 1.5 の中に，ENIAC で実施された弾道計算の差分式が示されているが，使われている数値積分公式は単段・刻み幅固定・陽解法のオイラー法であったと思われる．それでは，これから簡単にそれぞれの分類軸の意味とその特性について解説してみる．

雑感

ENIAC で実行された数値計算について（121）75

5. 1 陽解法と陰解法
陽解法は， 番目の積分値を=-1+ Δ という形で計算していくが，このとき増分Δ をその時点で既知である - 1 番目以前の積分値と微分係数値だけを使用して陽的に算出できるため，計算は容易であるが，積分刻み幅 Δ を大きく取りすぎると誤差が発散する危険性をもつ．一方，陰解法では，Δ を計算する際にこれから求めるべき積分値（ ）を必要とするため，Δ を陽的には計算できず，番目のを未知数として陰に含む方程式を解きながら積分を進める必要がある．微分方程式が連立していて，積分値 がベクトル値となる場合には，1 ステップごとに連立非線形方程式をニュートン・ラプソン法などの繰返し法によって解かねばならない場合もありプログラムの実装は面倒である．半面，積分刻み幅をいくら大きくとっても誤差が発散しないので，電子回路など，きわめて短周期で微小振動する変数を含みながら安定的に静定しているために長時間分の計算を必要とするいわゆるスティフな系（硬い系）の解析では必須な数値積分法である．スティフな系を陽解法で積分すると，固定刻み幅法では解が発散し，可変刻み幅法では刻み幅が極限まで縮小されて計算が停留してしまうなどの現象に陥り，まったく実用性が無くなってしまうという事実は，筆者もかつて実務で体験し，陰解法積分を実装して切り抜けた記憶がある．
5. 2 単段法と多段法
積分計算を一つの時間刻み幅分実行する際に， 番目の計算値と - 1 番目の計算値だ
けを使用して計算を完結することができる積分法を単段法とよび，一方 - 1 番目以外に，- 2 番目， - 3 番目，…を使用しないと計算を完結できない方法を多段法と呼ぶ．単段法は単純だが，Δ を大きく取ると積分の精度を上げるのが難しくなり，ルンゲクッタ法などのように，1
Δ に於ける微分係数の評価を計算過程に加えるなどの工夫が必要であ2

る．反面，多段法は履歴の管理などが煩雑になる代わりに，過去の履歴を利用して，高次の近似多項式を用いた高精度で効率の良い積分計算を，分かりやすく実装可能である．
5. 3

固定時間刻み幅法と可変時間刻み幅法

初期の積分計算法は，固定の時間刻み幅について実施されていたが，この刻み幅をどういう値に設定するかで，積分計算の精度や計算時間に大きな影響が及ぶため，経験の無い利用者には使いづらい技術であった．この点を改良したのが，可変刻み幅積分法である．
この手法では，積分計算の過程で派生する情報をもとに，各計算ステップにおいて積分計算の打切り誤差を推定し，計算開始時に利用者から与えられた許容誤差情報と比較することにより，時間刻み幅を広げたり，縮めたりして，許容誤差範囲内の打切り誤差を実現しながら計算を進める．可変時間刻み幅を実現する積分手法として知られているものに，予測子修正子法と呼ばれる方法がある．この方法では，まず過去の数段階の積分値列に予測子公式を適用して次の積分値を予測し，予測値を用いてその時間の微分係数を評価しその値と過去の履歴に対して，修正子公式を適用し，予測値を修正する．この時，予測子と修正子の差を，積分誤差の推定値として利用する．他にも，単段法のルンゲクッタ・ギル法を改良して，時間幅 Δ で行った積分と，時間幅を1Δ として 2 回行った積分との計算
2

76（122）

結果を比較することにより，時間幅を自動調節する方法も利用されている様である．
6.

気象シミュレーションと偏微分方程式の数値解法
本書の第 10 章に，「ENIAC を用いて行われた最も有名な計算は気象の数値シミュレーションである」という書き出しの一節がある．地球を取り巻く大気の挙動を，流体力学解
析で予測しようという試みは，電子計算機がこの世に産声を上げたときから，その最も重要な適用領域であり，また現代においても，スーパーコンピュータと呼ばれる最新型の超大型並列計算機の重要な適用問題の一つであり続けている．気象の挙動を正確にモデル化し，それを計算することによって天候の変化を正確に予測できるようになったときに，人類が得られる利益は計り知れないものになるのであるから，当然である．然るに，今日現在においても，数値気象予測情報は，気象予報士の大事な情報源の一つではあるものの，予報士を失業させるには至っていない．気象解析の場合は，広大な地球表面全体を解析領域とするため，現代の最大級の計算機の記憶容量と計算速度をもってしても，計算できるのは，1km 間隔の離散点上における 1 時間に数回程度の間隔の気圧・風向・風速の変化だけである．したがって局地的な雷雨や竜巻などせいぜい数 km の範囲で激しい時間的変化を伴う現象を予測するには，まだまだ計算点も計算時間間隔も粗過ぎて精度が不足なのである．反面，大域的・長期的な気候変化，例えば偏西風のうねりの変化やそれに伴う極地寒気団の動向の予測を根拠とした気候の長期的な予報の精度向上には著しいものがある．これには，計算機の能力の向上による計算点の増加や計算時間間隔の短縮も大切であるが，もっと重要なのが，気象モデルの境界条件の精度の向上と思われる．気象モデルの境界条件は，地球表面の温度分布と大気圏最上部の温度分布が重要であるが，これらの値は，従来は海洋測候船やラジオゾンデから得られたまばらな情報を補間して入力するしかなかったが，近年は気象衛星の赤外線画像から画素単位の解像度で補間無しに入手できるようになった．
気象解析の本質である，流体解析のベースとなるのは，ナビヤー・ストークス方程式と呼ばれる偏微分方程式である．
 1
Dv
1= F - p +   + v
3Dt
ここでæ¶D¶¶¶ö
v+w + ÷÷÷ ：全微分演算子= ççç u +Dt èç ¶x¶y¶z¶t ø÷v = (u, v, w) ：速度ベクトル， F = ( Fx , Fy , Fz ) ：外力ベクトル，：圧力，ρ：密度，  =æ ¶u = div (v ) = ççççè ¶x+：動粘性係数，μ：粘性係数¶v ¶w ö÷÷÷ ：圧力による歪+¶y¶z ÷ø雑感

ENIAC で実行された数値計算について（123）77 ¶ ¶ ¶   = ex+ ey+ ez：ナブラ演算子（ ex , ey , ez は各座標軸方向の¶x¶y¶z単位ベクトル）æ ¶2¶2¶ 2 ö÷
 = 2 = ççç 2 + 2 + 2 ÷÷ ：ラプラス演算子¶y¶z ÷÷øèç ¶x

この方程式は，大変複雑な方程式であるため，解析する流体の特徴に応じて簡略化されて解析するのが通常である．気象解析の場合は，非圧縮性を仮定して，右辺第 4 項（圧縮項）を無視し，密度 ρ を定数とするのが通例である．
偏微分方程式は，解析領域を離散化して，離散計算点上の状態変数を求める連立方程式として数値解析するのが普通である．離散化の代表的なものとして，次の三つをあげることができる．
①差分法②有限要素法③粒子法
本書 275 ページ図 10.4 に ENIAC で行われた気象解析の手順書が示されているが，「渦度」と「フーリエ変換」というキーワードから，使用された離散化手法は差分法であったことが想像できる．同じ図 10.4 の流れ図には，パンチカード出力とかパンチカード操作と分類されたボックスがある．モンテカルロ・シミュレーションの時と同様，気象モデルの多数の格子点に関するデータを計算機内に保持することができなかったので，離散化した偏微分方程式を一つの時間刻み幅分計算するたびに，その値を 1 格子 1 枚のカードにパンチ出力し，次のステップではそのカードを読み込みながら，次のステップの状態の計算を進めていったことが分かる．以下に，それぞれの離散化手法について，簡単に解説する．
6. 1 差分法
空間を等間隔の格子点で分割し，微分方程式の中に現れる微分演算子を，隣接する格子の状態変数間の差分式で近似した，離散化方程式に帰着させる方法である．
p - pn¶p= n+1xn+1 - xn¶x

この手法は，分かりやすく実装も容易なため最も広く用いられており，気象解析にもこの方法が用いられる事が多い．欠点としては，矩形の格子点でモデル化するため，解析領域が複雑な曲面で囲まれているような場合は，境界領域をサイコロで積み上げるように近似するので，そこでの誤差が解析結果に悪影響を及ぼしてしまう事が挙げられる．この誤差を低減するには格子間隔を狭めて格子点を増やす必要があるが，3 次元問題では，格子間隔を半分にすると格子点の数は 8 倍となるため，計算負荷の増加も考慮する必要がある．
6. 2 有限要素法
複雑な境界面を持った解析領域を，その境界に沿った多面体などを積み上げて表現し，

78（124）

多面体の頂点を計算点とする方法である．多面体内部を連続的に補間する解析的関数を導入し，その関数の微分関数で偏微分方程式の微分演算子を近似離散化する（図 4）．3 次元
CAD システムの普及に伴い，複雑な解析形状の計算機モデルが自由に定義できるようになり，形状の自動分割機能も改良され，有限要素法による偏微分方程式の解析は，機械部品の強度解析，エンジン内部の冷却水の流路解析など，多くの応用領域で活用されている．

図4

有限要素法による機械部品の離散化モデル

6. 3 粒子法
粒子法は，近年大変注目されている新しい離散化手法である．上述した差分法や有限要素法が，空間上に固定された計算点と，計算点同士の隣接関係を前提に離散化を展開するのに対し，粒子法では，解析空間内に任意にばら撒いた計算点（粒子）について，隣接関係を利用せずに離散化を進める．また，この方法の大きな特徴として，対象とする偏微分方程式は，空間に固定したオイラー座標系を用いて定式化するのではなく，計算点と共に座標系が移動するラグランジェ座標系上で定式される．流体解析で違いを説明すると，オイラー系では空間に固定した視点で，そこを通過する流体の状態変化に着目するのに対し，ラグランジェ系では，同じ流体粒子に着目しその時間的な動きを追い続けることになる．粒子法では，計算ステップごとに計算対象とする粒子の周辺に存在する粒子を検索し，それらの持つ情報を利用して，微分演算子を離散化近似する．この方法は，解くべき問題ごとに微分演算子の離散化方法を個別定義しなければならないという限界はあるものの，①計算点の隣接関係に依存しないため，格子が乱れるような大変形に対しても計算の整合性が保障される．
②計算点ごとの演算の独立性が高いため，近年流行の GPGPU などを利用した，超並列計算との適合性がよい．などの特徴が魅力的なため，多くの研究者が適用領域の拡大や計算速度の向上に努力中であり，その成果が期待される．
─────────
参考文献 ［ 1 ］ Thomas Haigh，Mark Priestley，Crispin Rope 著，土居範久監修，羽田昭裕，川辺治之翻訳，「ENIAC ─現代計算技術のフロンティア─」，共立出版，2016 年 6 月．
（ユニシス技報 第 6 代編集委員長）