大学の卒業研究でFortranのプログラムを書きました。（いつのこっちゃ.1986年）

現代制御理論の研究室です。

連立微分方程式でPade近似の論文を拝見しました（参考文献参照）。
手元のFortranのライブラリに入っていないとのことでした。
卒業研究の題材として書き始めました。

<この項は書きかけです。順次追記します。>

手作業による最適化

解法は、大きく前半と後半に分かれています。
最初は、それぞれの部分で必要のない計算はなるべくしないようにループなどを定義していました。
計算時間はだんだん速くなっていきました。

次に、前半だけで使っている変数と、後半だけで使っている変数を、同じ変数を使ってみました。
すると、プログラムの大きさが、ぐっと小さくなりました。

しかし、コンパイラの最適化を選ぶと、違う答えが出るようになりました。

同じ変数の処理の順番を最適化で変更したのだろうと想定できました。

二つの版の更新

どちらの方式がいいか理論的な勉強が不足していました。
ひとまず、変数を分けて、コンパイラの最適化を利用する版と、
同じ変数を使いまわしてコンパイラの最適化を選ばない版を更新しました。

しばらくすると、両方手を入れるのは、面倒だけでなく、間違いがあり、計算結果が違うと原因調査に時間がかかることがわかりました。

まだ当時は、自動生成のプログラムは書いたことがなく、一方から他方を自動生成するところまで至っていませんでした。
（自動生成のプログラムを書き始めたのは翌年です）
自己再構築(self reconstruction)による算譜(program)自動生成器(generator)自動生成器(generator)の作成
https://qiita.com/kaizen_nagoya/items/0facd1ed443947e53cc9

手の最適化の版だけ更新することにしました。

誤差の評価

Fortranの計算で、誤差の評価をしようとしました。
倍精度で計算していたものを、四倍精度に変更して計算してみることにしました。

四倍精度で、全件の誤差評価の計算をさせると１日経っても、２日経っても、計算が終わりませんでした。
そこで、途中までの処理結果を出力するようにしました。

出力結果を見てみると、計算量から逆算すると、最終的に誤差評価が終わるのが３年後だと推測できました。

このプログラムの全数検査による誤差評価は中止しました。

行列指数関数の近似計算精度
https://qiita.com/ascmat/items/2c0055139730204d567a

python-control を用いたfeedback制御の基礎
https://qiita.com/hikaruyaku/items/021769a76e72ed34a881

今なら

MatlabのSymbolic Math Toolbox™ の関数 pade
https://jp.mathworks.com/help/symbolic/pade-approximant.html

今なら実施する誤差評価

四倍精度を計算しようとしたのが間違い。
普通の精度と倍制度の差の分布から、四倍精度と倍制度の誤差を推測するという方法を取っても良かったかもしれない。そして、部分を計算してみて全体を予測してもよかったかもしれない。　

そして、今なら全件計算してもたぶん１日で終わる。計算速度の工場で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
https://www.amazon.com/dp/0521450071

Padé approximant
http://www.scholarpedia.org/article/Padé_approximant

Padé Approximant
https://mathworld.wolfram.com/PadeApproximant.html

Padé approximation of models with time delay
https://www.mathworks.com/help/control/ref/pade.html

Padé Approximations。
https://www.cfm.brown.edu/people/dobrush/am33/Mathematica/ch3/pade.html

むだ時間の Pade近似と連分数展開による実現法
計測自動制御学会論文集 Vol. 21, No. 11, 昭和60年11月(1985)
https://www.jstage.jst.go.jp/article/sicetr1965/21/11/21_11_1171/_pdf

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コンパイラとモデルに基づく設計
https://qiita.com/kaizen_nagoya/items/5e50b3c9c0004fe92bd9

模型駆動開発(Model Driven Design)への道
https://qiita.com/kaizen_nagoya/items/bb4d73bfb3cbba88727f

一覧

物理記事　上位100
https://qiita.com/kaizen_nagoya/items/66e90fe31fbe3facc6ff

量子(0) 計算機, 量子力学
https://qiita.com/kaizen_nagoya/items/1cd954cb0eed92879fd4

数学関連記事１００
https://qiita.com/kaizen_nagoya/items/d8dadb49a6397e854c6d

統計(0)一覧
https://qiita.com/kaizen_nagoya/items/80d3b221807e53e88aba

品質一覧
https://qiita.com/kaizen_nagoya/items/2b99b8e9db6d94b2e971

言語・文学記事　１００
https://qiita.com/kaizen_nagoya/items/42d58d5ef7fb53c407d6

医工連携関連記事一覧
https://qiita.com/kaizen_nagoya/items/6ab51c12ba51bc260a82

自動車　記事　１００
https://qiita.com/kaizen_nagoya/items/f7f0b9ab36569ad409c5

通信記事１００
https://qiita.com/kaizen_nagoya/items/1d67de5e1cd207b05ef7

日本語（０）一欄
https://qiita.com/kaizen_nagoya/items/7498dcfa3a9ba7fd1e68

英語(0) 一覧
https://qiita.com/kaizen_nagoya/items/680e3f5cbf9430486c7d

転職(0)一覧
https://qiita.com/kaizen_nagoya/items/f77520d378d33451d6fe

仮説（0）一覧（目標100現在40）
https://qiita.com/kaizen_nagoya/items/f000506fe1837b3590df

音楽　一覧(0)
https://qiita.com/kaizen_nagoya/items/b6e5f42bbfe3bbe40f5d

「@kazuo_reve 新人の方によく展開している有益な情報」確認一覧
https://qiita.com/kaizen_nagoya/items/b9380888d1e5a042646b

Qiita(0)Qiita関連記事一覧（自分）
https://qiita.com/kaizen_nagoya/items/58db5fbf036b28e9dfa6

鉄道（０）鉄道のシステム考察はてっちゃんがてつだってくれる
https://qiita.com/kaizen_nagoya/items/26bda595f341a27901a0

安全（0）安全工学シンポジウムに向けて: 21
https://qiita.com/kaizen_nagoya/items/c5d78f3def8195cb2409

一覧の一覧( The directory of directories of mine.) Qiita(100)
https://qiita.com/kaizen_nagoya/items/7eb0e006543886138f39

Ethernet 記事一覧　Ethernet(0)
https://qiita.com/kaizen_nagoya/items/88d35e99f74aefc98794

Wireshark 一覧 wireshark(0)、Ethernet(48)
https://qiita.com/kaizen_nagoya/items/fbed841f61875c4731d0

線網（Wi-Fi）空中線(antenna)(0) 記事一覧(118/300目標)
https://qiita.com/kaizen_nagoya/items/5e5464ac2b24bd4cd001

OSEK OS設計の基礎　OSEK(100)
https://qiita.com/kaizen_nagoya/items/7528a22a14242d2d58a3

Error一覧 error(0)
https://qiita.com/kaizen_nagoya/items/48b6cbc8d68eae2c42b8

++ Support(0)　
https://qiita.com/kaizen_nagoya/items/8720d26f762369a80514

Coding(0) Rules, C, Secure, MISRA and so on
https://qiita.com/kaizen_nagoya/items/400725644a8a0e90fbb0

プログラマによる、プログラマのための、統計(0)と確率のプログラミングとその後
https://qiita.com/kaizen_nagoya/items/6e9897eb641268766909

なぜdockerで機械学習するか書籍・ソース一覧作成中 (目標100)
https://qiita.com/kaizen_nagoya/items/ddd12477544bf5ba85e2

言語処理100本ノックをdockerで。python覚えるのに最適。:10+12
https://qiita.com/kaizen_nagoya/items/7e7eb7c543e0c18438c4

プログラムちょい替え（0）一覧:4件
https://qiita.com/kaizen_nagoya/items/296d87ef4bfd516bc394

Python(0)記事をまとめたい。
https://qiita.com/kaizen_nagoya/items/088c57d70ab6904ebb53

官公庁・学校・公的団体（NPOを含む）システムの課題、官（０）
https://qiita.com/kaizen_nagoya/items/04ee6eaf7ec13d3af4c3

「はじめての」シリーズ　ベクタージャパン　
https://qiita.com/kaizen_nagoya/items/2e41634f6e21a3cf74eb

AUTOSAR(0)Qiita記事一覧, OSEK(75)
https://qiita.com/kaizen_nagoya/items/89c07961b59a8754c869

プログラマが知っていると良い「公序良俗」
https://qiita.com/kaizen_nagoya/items/9fe7c0dfac2fbd77a945

LaTeX(0) 一覧　
https://qiita.com/kaizen_nagoya/items/e3f7dafacab58c499792

自動制御、制御工学一覧（０）
https://qiita.com/kaizen_nagoya/items/7767a4e19a6ae1479e6b

Rust(0) 一覧　
https://qiita.com/kaizen_nagoya/items/5e8bb080ba6ca0281927

小川清最終講義、最終講義（再）計画, Ethernet(100) 英語(100) 安全(100)
https://qiita.com/kaizen_nagoya/items/e2df642e3951e35e6a53

＜この記事は個人の過去の経験に基づく個人の感想です。現在所属する組織、業務とは関係がありません。＞
This article is an individual impression based on my individual experience. It has nothing to do with the organization or business to which I currently belong.

文献履歴

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.