SSブログ
数値解析 ブログトップ
前の10件 | -

SIRモデルの改良 [数値解析]

SIRモデルの改良

 

 

ケルマックとマッケンドリンクによって1927年に提出されたSIRモデル

  ISIR-001.png

に出産率b、自然死亡率μを導入し、次のように修正する。

  ISIR-002.png

 

(2)の辺々を足し合わせると、

  

したがって、全人口に関して、次の式を得る。

  

とすると、安定な平衡値として

  

が得られる。

 

問 次の微分方程式を解き、

  

となることを確かめよ。

【解】

  

両辺にを掛けると、

  

両辺をtで積分すると、

  ISI$-003.png

μ>0だから、

  ISIR-004.png

となり、全人口の平衡値に一致する。

(解答終)

 

そこで、はじめから全人口はであると仮定し、

  ISIR-005.png

について考えればよい。

 

感染症の流行初期の段階では、S(t)≒Nだから、(4)の第2式は次のように線形化可能(マルサスのモデル)。

  ISIR-006.png

したがって、基本再生産数

  

と定義することができ、R₀>1ならば感染は拡大し、R₀<1ならば感染は自然消滅することになる。

 

(4)式でになる(S,I)の点――このような点を平衡点という。力学系の話なので、詳しい話は、物理担当のddt³さんに聞くとよい(^^ゞ――を求めると、

  

という定常解が得られる。

前者の(N,0)は感染者がいない場合の定常解である。後者の定常解はR₀>1のときにのみ意味をもつ解で、これは集団の中に感染者が常に存在するエンデミックを表す定常解である。

 ――シンプルなSIRモデルでは、R₀の値にかかわらずとなるので、エンデミックな状態を再現できないが、修正した新たなモデルでは、出産やヒトの流入などで、集団の感受人口が常に供給されるので、エンデミックが起きる。――

 

エンデミックの状態をとすると、

  

 

ネムネコは、感染症の流行モデルの研究者じゃなく、ただのネコなので、こんな抽象的な話をされても具体的なイメージがわかなくて、わかった気分になれない。

――あくまで”わかった”気分で、本当にここまでの内容を理解できているかどうかは、また、別の話!!――

だ・か・ら、

数値的であろうがどのような方法であろうが、(3)を解き、それをグラフにし、視覚化しないと!!

 

というわけで、例によって、出生率b、自然死亡率μ、感染率β、隔離率(回復率)γに適当な値を入れ、オイラー法を用いて、連立微分方程式(4)を解いてみたにゃ。

――オイラー法ならば、誰でも表計算ソフトを使うことで簡単に計算できる!!。根性を出せば、手計算や、そろばん、タイガー計算機などを使って、この問題を数値的に解くことさえできる。だって、この計算で使うのは、小学校で習う加減乗除の四則演算のみで、ルートや三角関数といった難しい数学は一切不要なんだから――

 

 

出生率b=1、自然死亡率μ=0.01とすると、安定な人口の平衡値Nは、

  

そして、β=0.005γ=0.1とすると、基本再生産数R₀

  

t=0における感受性人口S(0)を99、感染者人口I(0)=1とし、オイラー法を用いて(4)を解くと、次のような結果が得られる。

 

 

横軸にS、縦軸にIを取った位相図は次のとおり。

 

 

 

t=400におけるS(t)I(t)の値は

  

で、(7)式を用いて計算したエンデミックな定常状態

  

を再現していることがわかるにゃ。

 

「オイラー法は精度が悪い」と、多くのヒトがオイラー法を酷評するけれど、この問題の場合、意外に正確に計算できているにゃ。

だ・か・ら、

定性的に、そして、ある程度、定量的に(3)の解の挙動を知りたないのならば、簡単なオイラー法で十分だケロ。

数値計算の講義を担当する(大学の)先生方などは、決まって、「オイラー法は精度が悪いから使っちゃダメ。4次のルンゲ・クッタを使え」と、オイラー法を悪者扱いするけれど、オイラー法はそうそう捨てたもんじゃないってこと。

 

ではあるが、

4次のルンゲ・クッタ法を用いて、計算領域をt=500まで広げ、より高精度に計算した結果は次のとおり。

 

 

この計算結果によると、少し振動しながら、エンデミック状態に近づいてゆくみたいだね。

このことは、位相図を見ると、よくわかるにゃ。

 

 

反時計回りの渦(?)を描きならが、平衡点に収束するみたいだケロね。

そして、t=500くらいまで計算すると、先に求めたエンデミック状態の解(22.00, 7.09)に到達するのであった。

 

この図だと、うずまき具合がわかりづらいので、R₀=2.38の場合について、オイラー法で数値的に解いた結果を示すにゃ。

 

 


これだと、渦巻きを描きながら、パンデミックの定常解に近づくことがよくわかると思うにゃ。

 

記事の内容とは関係ないけれど、この動画も埋め込んでおこう。

 

 

例によって、計算に使用したスプレッドシートを公開しておいたので、興味ある奴はダウンロードし、値をいろいろ変えて計算してみるといいにゃ。

 

https://docs.google.com/spreadsheets/d/e/2PACX-1vSgi4iiUf2_0dDm6Fo6En7zTylPnSfjEKrUYHJ6-sJ9m1yBdotJ_PCGOqERQjFLr0iNNrO03eUWNhyc/pub?output=xlsx

 

https://docs.google.com/spreadsheets/d/e/2PACX-1vSgi4iiUf2_0dDm6Fo6En7zTylPnSfjEKrUYHJ6-sJ9m1yBdotJ_PCGOqERQjFLr0iNNrO03eUWNhyc/pub?output=ods


nice!(0)  コメント(0) 

感染症の流行モデル SIRモデル [数値解析]

感染症の流行モデル SIRモデル

 

 

つい最近知ったのだけれど、感染症の流行過程をあらわす古典的なモデルにSIRモデルというものがあるそうだ。

 

Sを感染感受者(非感染者)、Iを感染者、Rを回復者(隔離者)とするとき、SIRモデルは次の連立常微分方程式で表される。

  SIR-001.png

ここで、βは感染率、γは回復率(隔離率)であり、tはある時刻を起点にした時間。

 

(1)式の第1式、第2式、第3式の辺々を足し合わせると、

  

になるので、

  

そこで、

  

とおこう。Kは、人の出入りがない、孤立した地域の総人口だと考えればいい。

 

特に、感染症の流行初期においては、感染者が少なく、

  

と近似することができ、感染者数は次のマルサスモデルに従う。

  

よって、流行初期において、感染者が増加するためには

  

でなければならない。

このパラメータR₀を基本再生産数とよび、R₀>1のとき、流行初期の段階で感染は拡大し、感染者数は指数関数的に増加し、R₀<1ならば、感染の流行は発生せず、感染者数は下の図のように自然に減衰してゆく。

 

 

連立常微分方程式(1)は、一般に、解けない。より正確に表現すると、(1)の解であるS(t)I(t)R(t)を初等的なtの関数の組み合わせで表すことができない。

しかし、γ=0で、R(0)=0のとき、

  

となるので、これから

  

となり、これを(1)の第2式に代入すると、

  SIR-05.png

ここでとおくと、

  SIR-06.png

という解が得られる。

したがって、

  

 

回復率(隔離率)γ=0のときは、住民の全てが感染しまう。

一度感染すると誰一人として回復することがなく、しかも、感染者は非感染者に永遠にうつし続けるのだから、この結果は当然。

 

なお、(6)式の分子分母をN₀で割り、βN=r、さらに、I(t)=Nとおくと、(3)式は、次のロジスティック関数になる。

  SIR-07.png

このとき、Kは環境収容力とか何とか呼ばれるにゃ。

 

ネムネコは、少し前に、新型コロナウイルスの感染者数は、概ね、ロジスティック関数(7)に従って増加し続けると記事に書いたと思うけれど、これは回復率γ=0としたSIRモデルの解ということになる。

この記事はうっかり消してしまったようだけれど…。

 

上の図ではK=10⁵と書いてあるけれど、K=5×10⁵が正しいにゃ。

中国当局の発表によると、中国での感染者数は約8万人なので、ネムネコの予想より随分と少ないけれど、実際の感染者は発表の10倍ともいわれているので、と苦しい言い訳をする。

 

S(t)I(t)R(t)の関係ならば、

(1)の第2式を第1式で割ると、

  

この両辺をSで積分すると、

  

t=~のとき

  

となるので、

  

よって、

  SIR-08.png

また、(2)式より

  

となるので、

  SIR-09.png

 

これを解いたといっていいのか、微妙だけれど(^^ゞ。

 

また、感染者R(t)が最大になるのは、

 SIR-model-eq-40.png

の時だにゃ。

この時の感受性保持者の数は、γとβの値だけで決まるんだケロよ。

 

さて、γ=0といった特殊なケースでないと、非線形の連立常微分方程式(1)の解をtの関数として求めることはできないので、これを数値的に解こうじゃないか。

そこで、

  

という初期条件で、

感染率β=0.05に固定し、回復率γ=01、2の場合について――この数値が適当なものであるかどうかについては知らない!!ーー、

オイラー法を用いて、Δt=0.1として、数値的に解いてみたにゃ。

 

 

 

 

ここでは、ウィキペディアの記事にならって、γを回復率としたけれど、隔離率、R(t)を隔離された人口(病気からの回復による免疫保持者、隔離者、死亡者)とするのが正式みたいだケロね。

なお、ここでいう隔離者人口は感染症に発症した人、病院や隔離施設に隔離された人のことではないみたいだから、注意するにゃ。

また、感染者数R(t)は、時刻tの時点において感染症に罹っている人数のことで、時刻tまでに感染症の人数、つまり、累積感染者数ではないので、この点も注意するにゃ。

 

ところで、

先述したように、γ=0のときI(t)はロジスティック関数になるので、この数値計算の結果の妥当性を調べるために、β=0.025Δt=0.1ときのオイラー法に数値解とロジスティック関数との比較を以下に示す。

 

 

ロジスティック曲線は、tがあまり大きくないとき、指数関数的な挙動をするので、Δt=0.1という粗い計算だと、どうしても誤差が大きくなってしまうけど、定性的にはロジスティック曲線を再現しているにゃ。

Δt=0.05くらいにすると

  

と、結構、合うようになるにゃ。

 

ちなみに、オイラー法を用いた計算法は、

  

とし、連立微分方程式をつぎの漸化式に近似し、k=1,2,3,・・・に対して、逐次計算すればいいにゃ。

  

計算精度はともかく、オイラー法ならば、プログラムを作ることなく、表計算ソフトで、簡単に、非線形の連立常微分方程式を解くことができるにゃ。

 

そんな物好きがいるかどうかわからないけれど、計算に使用したスプレッドシートをネットで公開したので、興味のある奴はダウンロードし、β、ガンマの値などを変化させたとき、計算結果がどうなるか確かめてみるといいにゃ。

アドレスは

 

https://docs.google.com/spreadsheets/d/e/2PACX-1vQbAtk0oqiyJWLnIxqSrbIS8kMPoLX4D78IWk6hXsnhjSMryHs-9liX5Y4R68Sq-D6YyYfMGzkxJj2V/pub?output=xlsx

 

 

オイラー法を用いた数値解は、何故か知らないけれど、

  

を満たしいて、定性的に、そしてある程度定量的に(1)の解の挙動を知るにはこれで十分でしょう。

 

より高精度な計算をしたい奴は、Δtを小さくするのではなく、2次や4次のルンゲ・クッタ法を用いて解くといい。

精度には4次のルンゲ・クッタ法を用いたほうがいいけれど、この程度の微分方程式ならば、計算精度的に2次のルンゲ・クッタ法で十分だと思う。

 

 


nice!(2)  コメント(3) 

「数値解析勉強中の大学生」さんからいただいた質問への回答 [数値解析]

「数値解析勉強中の大学生」さんからいただいた質問への回答

 

 

2)k2の式で、

  

となっていますが、何故にはではなくが使われるのかが完全に理解できておらず、ご教授頂きたいです。宜しくお願いいたします。

 

 

  

に対する、ルンゲ・クッタ法は

  

ですよね。

 

時刻t陽に含まない、

  

の場合――力学系のといいます――は、これに対するルンゲ・クッタ法は

  

になります。

 

これは、次のようなベクトル表記を用い、そのまま、連立常微分方程式に拡張が可能です。

すなわち、

  

これに対する、ルンゲ・クッタ法は

  

これを成分で表すと、

  hrk-007.png

に対するルンゲ・クッタ法は、i=1,2,・・・、nに対して、

  hrk-009.png

となります。

 

そして、

  

の場合は、

  

になるでしょう、って話なんですが。

 

いま振り返ると、

  

ではなく、

  

とし、これから、偏微分方程式(1)は

  

という連立常微分方程式に変換――というか近似――できる。

そして、

  

とおくと、

  

になり、これにルンゲ・クッタ法を適用すると、

  

と書いたほうが良かったのかもしれませんね。

 

この話、近似法には、正直、胡散臭いところがいくつかあるもので、この部分の話はあまりしたくない(^^ゞ。

 

それで、質問の

何故にはではなくが使われるのかが完全に理解できておらず、ご教授頂きたいです。

ですが、

時刻t+Δt/2の変化分(の近似値)はだから、にこれを加えるみたいは話なんですよ。

同様に、 時刻t+Δt/2の変化分(の近似値)はだから、にこれを加え、

時刻t+Δt/2の変化分(の近似値)はだから、にこれを加えるといった話です。

 


nice!(0)  コメント(1) 

ネムネコ、リープ・フロッグ(蛙飛び)法で単振動を計算する [数値解析]

ネムネコ、リープ・フロッグ(蛙飛び)法で単振動を計算する

 

物理で最も基本的な運動の一つに単振動(一次元調和振動子)というものがある。

ニュートンの運動方程式は

ここで、mは(質点の)質量、ωは角振動数(角速度)であり、xは(質点の重心の)位置、tは時刻である。

m>0なのだから、この方程式は

としてもいいよね。

 

さて、

とおくと、(2)式は次のようになる。

したがって、2階微分方程式(2)と次の1階の連立微分方程式は同じもの。

 

微分の合成公式から

になるので、(6)式は

これは何かといえば、(力学的)エネルギー保存則を表しているにゃ。

 

(2)式をさらに一般化し、

とすると、(7)は次の連立微分方程式になる。

 

そして、今回紹介するリープ・フロッグ(蛙飛び)法は、

という初期値問題を

という漸化式に書き換え、(10)、(11)から求まるを微分方程式の解の近似解にするもの。

 

(10)に関しては、テーラー展開

から直ちに出るし、と考えれば、微分方程式(9)は

となり、これを積分すると、

この右辺の積分を台形公式で近似すれば

となるので、

になる。

したがって、(10)、(11)の(離散化)誤差はO((Δt)³)で、修正オイラー法、2次のルンゲ・クッタ法と同程度ということになる。

 

通常、これとは違った導出をするのですが、同じ結果が得られるので、構いやしないや(^^ゞ。

 

Δt=0.1として、t=20までの200ステップを、このリープ・フロッグ法、蛙飛び法を用いて、

 

を解いてみたにゃ。

 

こちらがその結果。

 

 

 

特に見て欲しいのが、横軸にx、縦軸に速度vをとった位相図。

 

厳密解は、

だから、

になるのですが、リープ・フロッグ法で解いた真円と見間違うほど綺麗な円を描く。


さらに、t=200までの2000ステップまで伸ばし計算を進めても、ルンゲ・クッタ法などとは違って、リープ・フロッグ法はエネルギー保存則を常に一定の誤差範囲の中で満たすので綺麗な円軌道を描く。

 

 

これに対して、リープ・フロッグ法と同程度の離散誤差をもつ2次のルンゲクッタ法を用いてこの問題を解いた場合、時間の経過とともに、誤差の伝播のためにエネルギーが増加し、軌道がx²+v²=1から徐々に徐々に離れていく。

 

 

スゴイにゃ、リープ・フロッグ法。


nice!(1)  コメント(0) 

どれが修正オイラー法なの? [数値解析]

どれが修正オイラー法なの?

 

 

微分方程式の初期値問題

  

があるとする。

 

ねこ騙し数学では、

とするとき、

  dorega-001.png

を修正オイラー法(改良オイラー法)、

  doraga-002.png

を2次のルンゲ・クッタ法と呼んでいるが、(2)をホイン(Heun)法、(3)を修正オイラー法と呼ぶヒトがいるなど、その呼称はヒト、書籍などによって異なっているようである。

 

そして、新たに、修正オイラー法(?)と呼ばれる方法を紹介する。

 

まず、オイラー法を用いてにおける(1)の厳密解を予想する。

  

この値を求め、次のように修正する。

  

さらに、この値を用い

  

以下、同様に、

  

と計算を繰り返し、収束値

  

とする方法。

 

ネムネコが大学院時代から使っている数値解析の本ーー結構、有名な数値解析の古典で、教科書(海外の教科書の翻訳)ーーでは、この方法をホイン法(オイラー・ガウス予測修正子法)と呼んでいる。

混乱の極みにある。

そして、ネムネコは、カビが生えるくらい古臭いこの教科書ーーサンプルプログラムは、いつの時代の「ふぉーとらん」かわからないくらい古いもの(だって、この本が書かれたのは1960年代!!)で、そのサンプルプログラムを現代のFortranでコンパイルすると、エラーが出る(笑)ーーの呼称に従っている。リスペクトだにゃ。

 

さて、例によって、次の初期値問題を例にこの方法を具体的に紹介することにする。

 

  

 

この微分方程式の解がであることは言うまでもないだろう。

 

さて、(7)をオイラー法で離散化すると、

  

になる。

修正オイラー法(2)の場合は、

  

2次のルンゲ・クッタ法も同じく

  

そして、厳密解の場合は

  

(8)、(9)、(10)の括弧の中を注目して欲しいのだけれど、オイラー法は(10)式の2次以上の項を切り捨てたもので、修正オイラー法3次以上の高次の項を切り捨てたものになっていることがわかるだろう。

つまり、このことからも、(8)式の局所的な(離散化)誤差の程度はO(h²)、(9)の局所的な(離散化)誤差の程度はO(h³)であることが確かめられる。

 

ここで新たな修正オイラー法の場合を説明する。

まず、オイラー法(8)でを推測するので、

  

次に、(5)式で修正するので、

  

さらに、

  

となり、

  

になる。

(10)と(11)を比較すると、(10)の3次の項h³/6、(11)の3次の項はh³/4で食い違っているので、極限値(12)にどれくらいの意味があるのかは疑問だが、h=0.1のとき、

  

になり、これから

  

であることがわかるにゃ。

ちなみに、修正オイラー法は

  

で、だから、今回、あらたな修正オイラー法(ホイン法?)の方がより厳密解に近いことがわかる。

 

問 次のことを示せ。

  doraga-006.png

 

次に、h=0.1として、新たな修正オイラー法(ホイン法?)で、y(0.1)の近似値を求めるにゃ。

  

この問題の場合、2、3回程度修正すれば、よいことがわかる。

よりも真実の値から返って遠ざかり、精度を悪化させるので、この問題の場合、3回以上の修正は無意味だね。

次のを求めるためには、で得た1.10525とし、

  dorega-008.png

とする。

 

この新しい修正オイラー法(ホイン法?)は、厳密解との局所的な誤差を(2)式で定式化される修正オイラー法の1/2程度に抑えることが出来るんだケロ。

 

 

すごいケロ。

 

でも、だから、の方がよりもいい値なんだにゃ。

何とも皮肉な話である。

 

この他にも、オイラー・リチャードソン法なんてのもある・・・。

 

Improved Euler's Methodの訳語も、修正オイラー法、改良オイラー法の2種類があり、さらに、混乱に拍車をかけているように思う。

ホント、困ってしまう。

 

 Improved Euler's Method : 改良オイラー法(?)
 Modified Euler's Method :  修正オイラー法(?)

 

こう翻訳した(・・?

 

ネット上に存在する、 大学の先生が書いた数値計算の記事でも2次のルンゲ・クッタ法を修正オイラー法と呼んでいるものがあるので、これらの記事を読むとき、十分、気をつけたほうがいいケロよ。


nice!(1)  コメント(0) 

挿入法 [数値解析]

挿入法

 

基本選択法よりも、さらに、効率のよい整列法である挿入法について、{4,2,5,1,3}を例に説明する。

 

まず、{4,2,5,1,3}の部分列(部分集合)である{4,2}に注目し、最後尾の要素2の挿入箇所を探し、そこに2を挿入すると

  {4,2}→{2,4}

になり、{4,2}の並び替えが終了する。

次に、{2,4,5}の最後尾の5の挿入箇所を探し、そこに1を挿入する。すると、次のようになる。

  {2,4,5}→{2,4,5}

同様に、

  {2,4,5,1}→{1,2,4,5}

最後に、最後尾の3の挿入箇所を探し、そこに3を挿入すると、

  {1,2,4,5,3}→{1,2,3,4,5}

と、{4,2,5,1,3}を昇順に整列することができる。

 

一般に、

 

の場合、

i=2,3,・・・,nに対して、

部分列におけるの挿入箇所を探索し、そこにを挿入することによって、整列することができる。

 

この手順を、アルゴリズム的に書くと次のようになる。

 

do i=2, 3, ・・・, n

 do j=i, i−1, ・・・, 2

  a(j−1)>a(j)ならば、a(j−1)a(j)を交換

 end do j

end do i

 

この方針に基づいてプログラムを書くと、次のようになる。

 

parameter(n=5)
integer a(n)
integer w

a(1)=4; a(2)=2; a(3)=5 ; a(4)=1; a(5)=3

do i=2,n
    do j=i,2,-1
        if(a(j-1).gt.a(j)) then  ! a(j-1)とa(j)の比較
            w=a(j); a(j)=a(j-1); a(j-1)=w  ! a(j-1)とa(j)の交換
        end if
    end do
end do

write(*,*) a

end
 

 

このプログラムで昇順にデータを整列できるのですが、たとえば、

 {2,4,5}

のように整列の必要のない部分列に関しては、整列の処理を省くことができるので、次のようにプログラムを改良することできる。

 

parameter(n=5)
integer a(n)
integer w

a(1)=4; a(2)=2; a(3)=5 ; a(4)=1; a(5)=3

do i=2,n
    do j=i,2,-1
        if(a(j-1).gt.a(j)) then  ! a(j-1)とa(j)の比較
            w=a(j); a(j)=a(j-1); a(j-1)=w  ! a(j-1)とa(j)の交換
        else
            exit  ! 並び替えの必要がないのでループから抜ける
        end if
    end do
end do

write(*,*) a

end

 

 

これで無駄な、無意味な挿入箇所の探索省くことが可能となり、より高速に並び替えを行うことができるのですが、a(j−1)a(j)の交換処理には

  w=a(j); a(j)=a(j-1); a(j-1)=w

と3つの作業が必要になるので、まだまだ、無駄が多い。

 

そこで、アルゴリズムを次のように変更する。

 

do i=2, 3, ・・・, n

 t=a(i)

 do j=i, i−1, ・・・, 2

  a(j−1)>tならば、a(j)a(j−1)に代入 (a(j−1)を1つ後ろにずらす)

  そうでなければ、ループを抜ける

 end do j

 ta(j)に代入

end do i

 

このように変更すると、

  w=a(j); a(j)=a(j-1); a(j-1)=w

 a(j)=a(j−1)

と1つの作業にすることができる。

これで、より、高速に並び替えを行うことができる。

 

parameter(n=100)
integer a(n)
integer t

do i=1, n
    a(i)=nrandom(n)  ! 乱数を使って1〜nまでのデータの作成
end do

do i=2,n
    t=a(i) ! a(i)をtにコピー
    do j=i,2,-1
        if(a(j-1).gt.t) then  ! a(j-1)とtの比較
            a(j)=a(j-1)  ! a(j-1)を後方に1つずらす
        else
            exit  ! 並び替えの必要がないのでループから抜ける
        end if
    end do
    a(j)=t  ! tを挿入
end do

write(*,100) (a(i), i=1,n)

100 format (10(i8))

end


! 1〜nまでの乱数を作成
function nrandom(n)
    call random_number(x)
    nrandom = n*x+1
end


プログラムの実行結果



nice!(0)  コメント(0) 

基本選択法 [数値解析]

基本選択法

 

前回紹介したバブルソート(Bubble Sort)は、隣接する2つの大小を比較し、順序が逆転していたならば、2つの要素を交換することのよって、データを並び替える方法であった。

そして、データの数がn個であった場合、比較回数が

  

回であり、データの交換回数が最悪

  

回になり、並び替えの効率が悪い方法であるという話をした。

 

隣接する2つの要素の大小関係を比較し、順序が逆転していたら交換することによって何をしているのかといえば、昇順ならば要素の最大の値のものを、降順ならば要素の値が最小のものを、その最後尾に配置してるのである。であるならば、昇順ならば要素の最大値を、降順ならば最小値を要素の中から選び出し、それを最後尾に配置すればよいということになるであろう。

 

たとえば、

  {4,2,5,1,3}

であれば、5が最大値なので、

  {4,2,5,1,3}→{4,2,3,1,5}

と並び替える。

5は最後尾に配置したので、次に、5を取り除いた

  {4,2,3,1}

の最大値4を最後尾に配置する。

  {4,2,3,1}→{1,2,3,4}

これで並び替えが終了。

 

{4,2,5,1,3}の場合、

 バブルソート:6回

 基本選択法:2回

と、データの交換回数を減らすことができる。

 

の場合、この方法のアルゴリズムは次のようなものになるであろう。

 

do i=n,n−1,・・・,2,1

 i_max=1

 do j=2,3,・・・,i

  a(j)>a(i_max)ならばji_maxにする

 end do j

 a(i_max)a(i)を交換

end do I

 

このFortranプログラムは次のようになる。

 

parameter(n=5)
integer a(n)
integer dummy

a(1)=4; a(2)=2; a(3)=5; a(4)=1; a(5)=3

do i=n, 2, -1
    i_max=1
    do j=2,i
        if(a(j).gt.a(i_max)) then
            i_max=j
        end if
    end do
    ! データの交換
     dummy=a(i) !データの退避
     a(i)=a(i_max)
     a(i_max)=dummy
end do

write(*,*) a

end

 

i_max=iのとき、交換をする必要がないので、次のように変更すると、さらに速くなる。

 

parameter(n=5)
integer a(n)
integer dummy


a(1)=4; a(2)=2; a(3)=5; a(4)=1; a(5)=3

do i=n, 2, -1
    i_max=1
    do j=2,i
        if(a(j).gt.a(i_max)) then
            i_max=j
        end if
    end do
    ! データの交換
    if(i_max.ne.i) then
        dummy=a(i) !データの退避
        a(i)=a(i_max)
        a(i_max)=dummy
    end if
end do

write(*,*) icount
write(*,*) a

end

データを昇順、降順に並び替えるという作業は、結構、奥が深いんだケロ。

そして、現在、最速とされるクイックソートよりも速い並べ替え法を発見すると、歴史に名を残すことができる。のみならず、特許をとれれば、巨万の富を稼げるかもしれない。

 

次回は、さらに速い挿入法を紹介する。

 


nice!(0)  コメント(0) 

バブルソート(泡立ち法) [数値解析]

バブルソート(泡立ち法)

 

いま仮に次のような数の列があるとする。

  

これを小さい数から大きい数へならぶ順(昇順)の数の列に変えたいとする。

この最もシンプルな方法は、隣り合う要素の大小関係を比較し、大小が逆転していたら、隣り合う2つの要素を交換するという方法。

具体的には、まず、{3,2,1}の(3,2)を比較する。3と2の大小関係が逆転しているので、3と2を交換し(2,3)にする。すると、

  

になる。

次に、{2,3,1}の(3,1)の大小関係を比較すると、順序が逆転しているので、3と1を交換し、(1,3)にする。

すると、

  

になる。

これで、最も大きい3を右端に並べられことができた。

 

ということで、{2,1,3}の部分列、部分集合

  

の要素を昇順に並べ替えればいいだろう。

(2,1)を比較すると、これは順序が逆転しているので順序を交換し、(1,2)にすると、

  

になる。

以上のことをまとめると、

  

と、隣り合う要素の大小関係を比較し、順序が逆転していたら両者を交換することによって、{3,,1}を昇順に{1,,3}にできる。

こういう数の列の並び替えをバブルソート(泡立ち法)という。

 

では、もっと複雑な

  

の場合について考える。

(4,2)は順序が逆転しているので、2つの要素を交換すると、

  

になる。

次に、(4,5)を比較すると、この順序は逆転していないので、交換せずにそのまま。

  

次に、(5,1)を比較すると、順序が逆転しているので、交換する。

  

次に、(5,3)を比較すると、順序が逆転しているので、2つの要素を交換する。

  

これで最も大きな5が右端に並べられた。

というわけで、今度は、

  

を同様に昇順に並べ替えればよい。

  {2,4,1,3}→{2,4,1,3}→{2,1,4,3}→{2,1,3,4}

これで、最大の4を右端に持ってくることができたので、4を取り除いた部分集合

  {2,1,3}

を同様に並び替えると、

  {2,1,3}→{1,2,3}→{1,2,3}

になる。

これですべて昇順に並び替えることができたのだけれど、{1,2,3}の最大値3を取り除いた部分集合

  {1,2}

を昇順に並び替える。

  {1,2}→{1,2}

これで並び替えの作業はおしまい。

これで、{2,4,5,1,3}を昇順に{1,2,3,4,5}と並び替えることができた。

 

さらに一般的につぎのような数の列、集合があるとする。

  

これを昇順に並び替える方法、アルゴリズムは次のようなものになるだろう。

 

do i=nn−1,n−2、・・・,2

 do j=1,2,・・・,i−1

  a(j)a(j+1)の大小を比較してa(j)a(j+1)を交換

 end do j

end do i

 

フォートラン的な記述ですが・・・。

 

これを元にFortranのプログラムを作ると、次のようになる。

 

! バブルソート
parameter(n=5)
integer a(n)
integer dummy

a(1)=4; a(2)=2; a(3) = 5; a(4)=1; a(5)=3

do i=n, 2, -1
    do j=1,i-1
        if (a(j).gt.a(j+1)) then !順序が逆転していたら交換
            dummy = a(j)
            a(j)=a(j+1)
            a(j+1)=dummy
        end if
    end do
end do

write(*,*) a

end

 

バブルソートを用いると、昇順にデータを並び替えられるのですが、この方法は、データの比較・交換の回数が最悪、

  

n²/2回。

n=100ならば最悪5000回、n=1000ならば50万回、n=1万ならば5千万回と、データの数nが大きくなると、データの比較・交換回数が大きくなりすぎるので、実際の並び替え、ソーティングでバブルソートが使われることはない。

 

まぁ、パソコンといえども、現在は処理速度、計算速度がトンデモなく速いから、データ数が1万くらいあってもバブルソートですぐに並び終えるので、データ数が1万くらいまでならば、最も基本的なバブルソートで十分に実用に足りると思うけれど、コンピュータの能力が低かったときは、バブルソートの遅さは致命的だったんだケロよ。

 

次回は、バブルソートよりも少しだけ早くソーティングが可能になる方法を紹介するにゃ。

 

 


nice!(0)  コメント(0) 

定積分の矩形公式と台形公式の誤差の問題3の解答例 [数値解析]

問題3(積分の第1平均値の定理)

fが閉区間[a,b]で連続、g[a,b]で非負連続ならば、

  

を満たすξが存在することを示せ。

【証明】

になるのはg=0のときで、

  

したがって、a<ξ<bである任意のξに対して(1)が成立する。

次に、g[a,b]で恒等的にg(x)=0でないとする。

fは閉区間[a,b]で連続なので、最小値mと最大値Mをもつ。

条件g≧0より、

  

だから、

  

が成立する。

fが定数関数でないとき、

  

で辺々を割ると、

  

したがって、中間値の定理より

  

を満足するξが存在する。

fが定数関数であるとき、(1)が成立するのは明らか。

(解答終)

 

追加問題

関数f(x)が有界閉区間[a,b]で連続であれば、

  

が存在することを示せ。

【解答】

f(x)が定数関数のとき(2)が成り立つことは明らか。

そこで、f(x)は定数関数でないとする。

f(x)[a,b]で連続なので、f(x)[a,b]で最小値mと最大値Mをもつ。

したがって、

  

f(x)は定数関数でないので、

  

辺々をb−a>0で割ると、

  

したがって、中間値の定理より、

  

を満たすξが存在する。

(解答終)

 

本によっては、(1)ではなく、(2)を積分の平均値の定理と書いてあるので注意。

 

 


nice!(0)  コメント(0) 

オイラー法の誤差解析 [数値解析]

オイラー法の誤差解析

 

次の常微分方程式の初期値問題について考える。

  

この微分方程式の解はy=y(x)であるから、(1)式の右辺は

  

xだけの関数になるので、(1)式は次のようになる。

  

したがって、

  

となる。

右辺の定積分を

  

と近似すると、局所的な離散誤差は

  

となる。

そして、この式から大域的な誤差の評価式

  

が得られる。

先の仮定より

  

だから、

  

 

次に、

  

として得られる大域的な誤差について考える。

  

ここで、Aは定数。

リプシッツ条件

  

を仮定すると、

  

ここで、

  

とおくと、

  

ε₀=0だから、

  

したがって、

  

となり、オイラー法の大域的な誤差は1次オーダーであることがわかる。

 

 

 


nice!(0)  コメント(0) 
前の10件 | - 数値解析 ブログトップ

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。