大学の卒業研究でFortranのプログラムを書きました。（いつのこっちゃ.1986年）
現代制御理論の研究室です。
連立微分方程式でPade近似の論文を拝見しました（参考文献参照）。
手元のFortranのライブラリに入っていないとのことでした。
卒業研究の題材として書き始めました。
<この項は書きかけです。順次追記します。>
手作業による最適化
解法は、大きく前半と後半に分かれています。
最初は、それぞれの部分で必要のない計算はなるべくしないようにループなどを定義していました。
計算時間はだんだん速くなっていきました。
次に、前半だけで使っている変数と、後半だけで使っている変数を、同じ変数を使ってみました。
すると、プログラムの大きさが、ぐっと小さくなりました。
しかし、コンパイラの最適化を選ぶと、違う答えが出るようになりました。
同じ変数の処理の順番を最適化で変更したのだろうと想定できました。
二つの版の更新
どちらの方式がいいか理論的な勉強が不足していました。
ひとまず、変数を分けて、コンパイラの最適化を利用する版と、
同じ変数を使いまわしてコンパイラの最適化を選ばない版を更新しました。
しばらくすると、両方手を入れるのは、面倒だけでなく、間違いがあり、計算結果が違うと原因調査に時間がかかることがわかりました。
まだ当時は、自動生成のプログラムは書いたことがなく、一方から他方を自動生成するところまで至っていませんでした。
（自動生成のプログラムを書き始めたのは翌年です）
自己再構築(self reconstruction)による算譜(program)自動生成器(generator)自動生成器(generator)の作成
https://qiita.com/kaizen_nagoya/items/0facd1ed443947e53cc9
手の最適化の版だけ更新することにしました。
誤差の評価
Fortranの計算で、誤差の評価をしようとしました。
倍精度で計算していたものを、四倍精度に変更して計算してみることにしました。
四倍精度で、全件の誤差評価の計算をさせると１日経っても、２日経っても、計算が終わりませんでした。
そこで、途中までの処理結果を出力するようにしました。
出力結果を見てみると、計算量から逆算すると、最終的に誤差評価が終わるのが３年後だと推測できました。
このプログラムの全数検査による誤差評価は中止しました。
行列指数関数の近似計算精度
python-control を用いたfeedback制御の基礎
今なら
MatlabのSymbolic Math Toolbox™ の関数 pade
今なら実施する誤差評価
四倍精度を計算しようとしたのが間違い。
普通の精度と倍制度の差の分布から、四倍精度と倍制度の誤差を推測するという方法を取っても良かったかもしれない。そして、部分を計算してみて全体を予測してもよかったかもしれない。
そして、今なら全件計算してもたぶん１日で終わる。計算速度の工場で1000倍にはなっている。
ソフトウェアテストで全件試験は不可能だと書いている文書が時々ある。
1986年の16bit CPU 40MHzでは、16bit値の全数検査は事実上不可能だったかもしれない。
当時でも、8bit変数は全件試験をかけていた。
理由は、どこに0割が発生するかわからなかったし、
どこで誤差が大きくなるかわかっていなかったから。
全数計算して、誤差分布を調べ、誤差が大きくなる計算は、誤差が大きくならない計算に分岐できるようにしたりと考えた。
p.s.
手計算で、いくつかの計算結果の誤差の妥当性の確認をしています。
参考文献・URL
Chudnovsky G.V. Hermite Pade' approximations to exponential functions and elementary estimates of the measure of irrationality of π, Lect. Notes Math 925, 299-322, Zbl518.4014, 1982
https://www.researchgate.net/publication/225553208_Hermite-pade_approximations_to_exponential_functions_and_elementary_estimates_of_the_measure_of_irrationality_of_p
full text申請中。
Chudnovsky G.V. Pade' approximations to solutions of linear differential equatios and applications to Diophantine analysis. Lect. Notes Math 1052, 85-167, Zbl 536.10029, 1985
http://www.pnas.org/content/pnas/80/16/5158.full.pdf
y(x) - \dfrac{P(x)}{Q(x)}
is regular at
x = \infty
Padé Approximants (Encyclopedia of Mathematics and its Applications, Series Number 59
Padé approximant
Padé Approximant
Padé approximation of models with time delay
Padé Approximations
むだ時間の Pade近似と連分数展開による実現法
計測自動制御学会論文集 Vol. 21, No. 11, 昭和60年11月(1985)
Fortran
https://qiita.com/search?page=3&q=FORTRAN
【fortran】行番号で制御はやってはならない
https://qiita.com/Bluepost59/items/83fe23cfae31dba99ecf
fortranで配列の平均値計算を高速化
https://qiita.com/ykatsu111/items/f8d421b315e27c6a7588
Fortran 2008 の編集記述子まとめ (出力時)
https://qiita.com/chirimen/items/f867855f8707ed4e104c
Fortran で zip 的連結
https://qiita.com/cure_honey/items/6d17a1a6696435abd241
自己参照
連立微分方程式のPade近似と、Cコンパイラとモデルに基づく設計
模型駆動開発(Model Driven Design)への道
https://qiita.com/kaizen_nagoya/items/bb4d73bfb3cbba88727f
＜この記事は個人の過去の経験に基づく個人の感想です。現在所属する組織、業務とは関係がありません。＞
文献履歴
ver. 1.00 初稿 20180227
ver. 1.01 参考文献追記 20180310
ver. 1.02 検算方法の案を追記、TeX表記追記 20180325
ver. 1.03 参考文献追記 20220716
ver. 1.04 参考資料追記 20221112
ver. 1.05 表題変更 20221224
ver. 1.06 ありがとう追記 20230312
最後までおよみいただきありがとうございました。
いいね 💚、フォローをお願いします。
Thank you very much for reading to the last sentence.
Please press the like icon 💚 and follow me for your happy life.