So-net無料ブログ作成
数値解析 ブログトップ
前の10件 | -

十進BASICとマチンの公式を用いて円周率πを100桁求める [数値解析]

十進BASICとマチンの公式を用いて円周率πを100桁求める

 

円周率πの近似値を100桁求めたところで何の役にも立たないけれど、十進BASICを使って、πの近似値を100桁求めてみた。

計算結果は次の通り。

 

Π=3.

14159 26535 89793 23846 26433 83279 50288 41971 69399 37510

58209 74944 59230 78164 06286 20899 86280 34825 34211 70679

82148 08651

 

すこし余分を見て、110桁まで表示してある。

 

この計算に使用したプログラムは次の通り。

 

REM Machiの公式を用いたπの近似計算

OPTION ARITHMETIC DECIMAL_HIGH

LET n=100

LET S=0

LET FMT$="#."&REPEAT$("#",110)

PRINT USING fmt$

FOR k=1 TO n

LET A=(-1)^(k-1)/(2*k-1)

LET B = 4/5^(2*k-1)

LET C = 1./239^(2*k-1)

LET S = S + 4*A*(B-C)

NEXT k

PRINT USING fmt$:s

END

 

 

これをプログラムと呼ぶのが恥ずかしい。

計算効率も悪いし・・・。

プログラムのわかりやすさ(可読性)を重視したということで、この点は大目に見てもらおう。

 

これくらいの計算回数ならば計算速度は問題にならないけれど、計算速度を重視するならば、FOR NEXTループのべき乗計算は絶対に避けるべき。

たとえば、

このべき乗計算は

という漸化式に書き換えることができる。

 

だから、たとえば、次のように書き換えたほうがいいと思うにゃ。

 

 

REM Machiの公式を用いたπの近似計算

OPTION ARITHMETIC DECIMAL_HIGH

 

LET N=100

LET S=0

LET S5=1/(5*5)

LET S239=1/(239*239)

LET f=1

LET B=1/5

LET C=1/239

 

FOR K=1 TO N

LET A=f/(2*k-1)

LET S = S + A*(16*B-4*C)

LET f=-f

LET B = S5*B

LET C = S239*C

NEXT k

 

LET FMT$="#."&REPEAT$("#",110)

PRINT USING fmt$

PRINT USING fmt$:S

 

END

 

 

FOR NEXTループの反復回数nを100万くらいすると、おそらく、計算速度の違いを体感できると思う。

 

このプログラムを使って、1000桁まで計算をしてみたのだが、996桁目の値が異なってしまった。十進BASICの1000桁の精度を誇る、計算拡張モードを使ったのだけれど、このあたりが十進BASICで計算できる円周率πの限界のようだ。

 


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

円周率の計算 その2 [数値解析]

円周率の計算 その2

 

§1 計算法の改良

 

前回、円周率πの近似計算に使用したライプニッツの級数(グレゴリーの公式)

  

は収束が非常に緩慢であった。

そこで、収束を速くする方法について考えることにする。

 

計算に使用する基礎式は、前回、同様に

  

 

いま、仮にとすると、

  

したがって、

  

そして、上の式と(2)式から次式を得ることができる((3)式はEulerが発見したもの)。


  

 

このようにして新たに得られた級数(4)を用いて円周率πを計算してみる。

 

 

  

と分解し、計算してある。

計算結果を見ればわかる通り、ライプニッツ級数とは異なり、わずか3、4回の計算でπ≒3.14が得られ、しかも、収束の速度が速い。

 

 

§2 マチンの公式

 

のとき、三角関数の倍角公式より

  

と、tan 4a1=tan π/4)に近い値になり、その差は1/119である。

そこで、4aπ/4の差を求めると、

  

となる。

したがって、

  

となる。

このようにして得られた

  

マチンの公式と呼び、円周率の(近似)値を求める際によく使用される。

 

(2)にx=1/51/239を代入すると、


  

これを(5)に代入すると、

  

が得られる。

この右辺の無限級数の部分和を計算することで、円周率πの近似値を求めることができる。

すなわち、

  

を用いて、円周率πの近似値を求めることができる。

 

(7)の右辺の級数の収束は速く、

  

と、πのよい近似値を与えることができる。

 

表計算ソフトを使い、(7)式で求めたπの計算結果を下に示す。

 

 

 

(6)の右辺の無限級数の部分和、すなわち、(7)をn=10項まで計算すれば、コンピュータで用いられるπの近似値に到達することがわかる。

 

なお、表中のABCは次の通り。

  

 


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

ネイピア数e=「ふなひとはちふたはち」の100桁計算に挑む [数値解析]

ネイピア数e=「ふなひとはちふたはち」の100桁計算に挑む

 

§1 ネイピア数eの100桁の近似値を求める

 

ネイピア数eの近似値は、つぎのマクローリン展開を利用することによって、簡単に計算できる。

  

したがって、n項で打ち切った場合の誤差は、

  

である。

0<θ<1なので、であるが、これをと考えると、n項で計算を打ち切った誤差は、

  

と見積もることができる。

したがって、小数点以下10桁まで正しく計算しようとすると、少なくともn=13項まで計算しなければならない。

理論上は、nを大きくすればするほどネイピア数eの真実の値に近づくはずだが、コンピュータは有限桁でしか計算をすることができないので、前回、表計算ソフトによる計算値を示したように、n≧16で計算してもこれ以上、精度の高いネイピア数eの近似値を求めることができない。

 

ということで、十進数1000桁の計算が可能な十進BASICを用いて、ネイピア数eの100桁までの近似値を求めてみることにした。

そして、これがその計算結果だにゃ。

 

e=2.71828 18284 59045 23536 02874 71352 66249 77572 47093  69995

    95749 66967 62772 40766 30353 54759 45713 82178  52516 64274

 

ウィキペディアにあるネイピア数eの1000桁表示の近似値と比較すればわかるけれど、100桁、正しく計算できているにゃ。

 

そして、この計算に使用したプログラムは、次の通り。

 

REM ネイピア数e100桁まで計算

OPTION ARITHMETIC DECIMAL_HIGH

 

LET N=80

LET A=1

LET S=1

FOR K=1 TO N

LET A=A/K

LET S = S+A

NEXT K

 

LET FMT$="#."&REPEAT$("#",100)

PRINT USING FMT$: s

END

 

こんな簡単なプログラムでこれほどの高精度計算ができてしまうのだから驚き。

十進BASICというソフトは、本当に、凄いよね~。

 

 

§2 ライプニッツの級数(グレゴリーの公式)を用いて円周率πの近似値を求める

 

関数は次のようにマクローリン展開できる。

  

また、

  

だから、

  

になる。

そして、(1)はtan xの逆関数tan⁻¹のマクローリン展開になる。

 

さらに、x=1のとき、(1)の級数

  

は収束するので、

  

となる。

この級数をライプニッツ級数グレゴリーの公式という。

 

このライプニッツ級数を使って、円周率πの近似値を求めようじゃないか。

 

表計算ソフトを使って求めようと一瞬考えたけれど、この級数の収束はとんでもなく遅いので、考え直し、十進BASICを使って、πの近似値を求めることにする。

 




計算回数と誤差を見ると、nを10倍にすると誤差が1/10になっていることがわかる。。つまり、この級数は1次収束なので収束がとんでもなく遅い。

このことは、円周率の次の桁の数字の正確な値を求めるための計算量は、およそ、その10倍の計算量になるということを意味する。さらに、1桁下の値を求めようとすると、さらにその10倍の計算量が必要で、2桁下の正確な値を求めるためには約100倍の計算量になってしまう。

そのくせ、n=10万にとっても小数点5桁目の値を正しく計算できない。

アルキメデスが紀元前3世紀に導いた

  

程度の精度を得るためには、n1000以上にする必要があるのだから、この計算法は、とてもじゃないけれど、採用できない。

 

なお、この計算に使用したプログラムは次の通り。

 

REM ライプニッツの級数を用いて円周率πの近似値を求める

INPUT n

LET t=1

LET s=0

FOR k=1 TO n

LET a=t/(2*k-1)

LET s=s+a

LET t=-1*t

NEXT k

PRINT "π(の近似値):",4*S

END

 

 


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

e=「ふなひとはちふたはち」を計算する [数値解析]

e=「ふなひとはちふたはち」を計算する

 

 

ネイピア数(自然底)の近似値を求めよと問われたら、多くのヒトがマクローリン展開

  

  

を用いてこの近似を求めるのではないか。

 

そこで、n=5として、計算してみると、

  

となり、小数点2桁まで正確に計算してくれる。

この時の打ち切り誤差は

  

程度と見積もることができる。

(正しくは、0.001615162であるが・・・)

 

n=6以上の計算は、さすがに手計算では辛いものがあるので、表計算ソフトで計算した結果は次の通り。

 

 

 

この計算結果を見ると、n=10程度にとれば、実用上、問題にならない精度で計算できることがわかる。

と同時に、n=16以上で、eの近似値を求めようとしても無意味なことがわかるに違いない。

n=16のときに誤差は0になっているが、eは無理数なので、表計算ソフトに出ているこの値が正しい訳ではない。

計算機特有の丸め誤差のために、これ以上という項をいくら増やして(この項の値はほとんど0)これを加えたとしても、計算結果に反映されないためだ。

 

さてさて、n=10まで計算するのは大変。

そこで、x=0.1n=5として、の近似値を求めるにする。

  

このときの打切誤差はおおよそ

  

だから、誤差は10⁻⁸オーダーだから、ほとんど無視できる程度。

そこで、先に求めたの近似値を10乗すると、

  

n=10まで計算したものよりも、このように計算したほうが正確に計算できる。

 

しかしながら、高校生はマクローリン展開を習っていない。したがって、このように計算することはできない。

そこで、高校で習う範囲で計算する方法を考えなければならない。

高校では、誤差のところで

  

という式は習うはずである。

そこで、とおき、x=0とすると、

  

という近似式を得ることができる。

この近似式の打切誤差は、だいたい、

  

程度だから、h=0.01にとれば、誤差はおよそ10⁻⁶オーダー、つまり、百万分の1になり、無視できるオーダーになる。

ということで、

  

そして、これを100回掛けると、

  

小数点以下4桁まで正確に計算することが出来た。

手回しのタイガー計算機と根性さえあれば、高校数学程度の知識でこれくらいの近似値を求めることができるに違いない。

 

 

 

大学院生のころ、配置された研究室には、これよりも旧式のタイガー計算機があった。これを初めて見たとき、凄いと感動し、使い方を少しだけ教えてもらった経験がある。そして「よく出来ているな〜」と感動した。

YouTubeの動画に出てくるタイガー計算機は、比較的、最近(1960年台)のモデルだと思う(笑)。

研究室にあったのは、1940〜1950年台のモデルだったように思う。

(画像元:https://goo.gl/YS2muZ

 

 

さらに、h=0.001にすれば、

  

これを1000回掛けると、

  

e=「ふなひとはちふたはち」を見事求めることが出来た!!

 

  

だから、9000回も掛け算をする回数を削減することに成功できたにゃ。

この差は大きいケロよ〜(笑)。

 

冗談はさておき、

コンピュータや電卓のない時代に、オイラーは

 

  e≒2.71828182845904523536028

 

という近似値を発表していたそうだ。

オイラーは(1)または(2)のマクローリン展開を用いてこの近似値を求めたのだろうと考えられているが、この近似値を出すためには少なくともn=23項まで計算する必要がある。

この計算で特に問題になるのが、

  

というアボガドロ定数に迫る大きな数の計算。

 (科学技術の計算に使用されるFortranC言語はこんな大きな整数を整数型として計算できない!!)

いやはや、いったい、どういう計算法で、オイラーさんは、この近似値を求めたのであろうか?

 


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

対流と拡散 第5回 一般公式 [数値解析]

対流と拡散 第5回 一般公式

 

J-node.png右図に示すように、距離δだけ離れた格子点ii+1を考える。

これらの格子点にある境界面を横切る全流束Jは次のように表すことができる。

ここで数P

であり、セルペクレ数である。

境界面におけるφの値はのある種の加重平均であり、にある数を掛けたものである。

したがって、次のように書くことができる。

ここで、αβPに依存する無次元数である。

よって、

と表すことができ、APの間には次の関係があることがわかる。

また、座標軸を逆にとれば、Pは−Pとなり、ABの役割が入れ替わる。したがって、A(P)B(P)は次式のような関係がある。

 

指数法から導かれるペクレ数Pに対するA(P)B(P)になる。

これを図に描くと次のようになっており、P=0に関して対称でになっていることがわかる。

 

 

P<0のとき、

となるので、Pのすべての値に対してA(P)は次のように書くことができる。

また、(2)式より

流束に関係する(1)を境界面ewに適用し、(4)と(5)を使えば、次の対流・拡散に関する一般公式を得ることができる。

ここで、

 

以下に代表的な方法のA(P)を示す。

 

 

 

 

 


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

対流と拡散 第4回 ハイブリッド法とべき乗法 [数値解析]

対流と拡散 第4回 ハイブリッド法とべき乗法

 

§1 ハイブリッド法

 

前回、定常1次元の対流・拡散をあらわす微分方程式

  hp-001.png

の厳密解を与える指数法を紹介した。

指数法とは、(1)を次のように離散化して、(1)の数値解を得る方法である。

  

(3)式の両辺をで割ると、次の式を得る。

  

横軸に、縦軸にをとり、(6)をグラフにすると、次のようになる。

 

 

このグラフには、

x→−∞のときの漸近線

  

x=0における接線

  

x→0のときの漸近線

  

も合せて破線として示してある。

この図から、(6)はこれらの3本の直線で近似できることがわかる。

すなわち、

  

これは

  

または、

  

と一本の式で表すことができる。

ここで、記号

  

は、ABCの最大値を表す。

数学の記号で書くと、

  

である。

 

この方法をハイブリッド法という。

ハイブリッド法は、のとき中心差分、のときは風上差分(上流差分)と一致する。

 

そして、ハイブリッド法による微分方程式(1)の離散化方程式は次のようになる。

  

 

 

§2 べき乗法

 

図を見ると、ハイブリッド法はの付近において厳密解とのズレが大きい。

そこで、次のように近似する。

のとき

  

のとき

  

のとき

  

 

のとき

  

この近似式に基づき、各に対してを次のようになる。

 

 

よく近似できていることがわかるだろう。

 

これをパタンカーのべき乗法という。

 

べき乗による(1)式の離散化方程式は次のようになる。

  

ここで、

  



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

対流と拡散 第3回 指数法 [数値解析]

対流と拡散 第3回 指数法

 

  

拡散係数Γが一定のとき(1)は解析的に解くことができ、境界条件が

x=0φ=φ₀

x=L

ならば、(1)の解は

  

になる。

ここでPは次式で定義されるペクレ数である。

  

 

ところで、

  

とおくと、(1)式は

  

となる。

これをコントロールボリュームで積分すると、

  

(1)式の厳密解のφ₀をそれぞれに、さらにLに置き換えて、代入すれば

  

ここで、

  

同様に、

  

(7)、(8)を(6)に代入すると、

  

これから、次の離散化方程式が得られる。

  

この離散化された方程式を用いて、(1)式の数値解を求める方法を指数法という。

 

拡散係数Γが一定のとき、丸め誤差などを無視すれば、指数法は、計算領域の分割数にの大小にかかわらず、定常一次元問題に関して、厳密解と一致する。

 

この指数法は、コンピュータの計算能力が低いとき、

 (1) 指数関数の計算に時間がかかる

 (2) 2次元、3次元問題、生成や消滅を表す生成項Sが存在するときに厳密でない

などの理由から、実際の数値計算の現場で使われることがなかった方法。

しかし、コンピュータの能力が飛躍的に向上した現在、(1)は指数法を使わない理由にならないよね(^^)

 

 


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

対流と拡散 第2回 [数値解析]

対流と拡散 第2回

 

質量流束ρuと拡散係数Γが一定の場合、

  

となるので、ⅰ次元定常の対流・拡散の輸送方程式

  

  

になる。

tai-kaku-fig-001.pngさらに密度ρと流速uを一定とし、とし、左辺の対流項に中心差分を用いると、(1)式は

  

となり、これから次の差分方程式を得ることができる。

  

ここで、

このことから、前回導いた(8)〜(11)式が中心差分と同等のものであることが理解できるだろう。

 

これに基づき、

  

を数値的に解いてみた。

計算領域[0,1]を10分割し、すなわち、δx=0.1で、計算している。

 

F=1D=0.1、すなわち、格子ペクレ数P

   

P=0.1のときの結果は次の通り。

 

Kekka-graph-001.png

 

この場合、厳密解とよく一致する。

F=5D=1とし、格子ペクレ数P=0.5、さらにF=10D=1、格子ペクレ数P=1のときの計算結果を以下に示す。

 

Kekka-graph-002.png

 

Kekka-graph-003.png

 

格子ペクレ数Pが大きくなるに従い、数値解と厳密解との差が大きくなってくる。

F=20D=1、格子ペクレ数=2のとき、対流項に中心差分を用いた数値解は、実質、意味を有さないものになってしまう。

 

Kekka-graph-004.png 

 

そして、格子ペクレ数Pを2より大きくすると、数値解は振動するようになる。

 

Kekka-graph-005.png

 

したがって、対流項に中心差分を用いる場合、格子ペクレ数Pを|P<2にとる必要がある。

 

対流項を風上差分(上流差分)を用いると、この振動は抑えられるが、中心差分の打切誤差が(Δx)³であるのに対し、風上差分の打切誤差は(Δx)²と、数値解と厳密解との誤差が大きいという欠点を有している。

 

 Kekka-graph-006.png

この計算に使用したプログラムを以下に示す。

 

parameter (n=10)
real t(0:n)
real ae(4), aw(4), ap(4)
real a(n),b(n),c(n),d(n)
real result(0:n,3)

t=0.
t(n)=1.0

dx=1./n
f=1.
gam=1.

de=gam/dx
pe=f/de
write(*,*) 'Cell Peclt Number', pe

 

! 中心差分   
ae(1)=de-0.5*f
aw(1)=de+0.5*f
ap(1)=ae(1)+aw(1)

! 風上差分
ae(2)=de+amax1(-f,0.)
aw(2)=de+amax1(f,0.)
ap(2)=ae(2)+aw(2)

! べき乗法
ae(3)=de*amax1(0.,(1.-0.1*abs(pe))**5)+amax1(-f,0.)
aw(3)=de*amax1(0.,(1.-0.1*abs(pe))**5)+amax1(f,0.)
ap(3)=ae(3)+aw(3)

! 指数法

if (pe.ne.0.) then
    ae(4)=de*abs(pe)/(exp(abs(pe))-1.)+amax1(-f,0.)
    aw(4)=de*abs(pe)/(exp(abs(pe))-1.)+amax1(f,0.)
else
    ae(4)=de+amax1(-f,0.)
    aw(4)=de+amax1(f,0.)
end if
ap(4)=ae(4)+aw(4)


do j=1, 3
    do i=1, n-1
        a(i)=-aw(j)
        b(i)=ap(j)
        c(i)=-ae(j)
        d(i)=0.
    end do
        d(1) = d(1)+aw(j)*t(0)
        d(n-1)= d(n-1)+ae(j)*t(n)

    call tdma(a,b,c,d,n-1) ! TDMAで連立方程式を解く

    do i=1,n-1
        t(i)=d(i)
        result(i,j) = d(i)
    end do
    result(0,j)=t(0)
    result(n,j)=t(n)
end do

do i=0,n
write(*,100) i*dx, result(i,1), result(i,2), result(i,3), exact(f/gam, i*dx)
end do

100 format(f9.5, 1x , 3(f10.7, 1x), f10.7)

end


!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do

! 後退代入
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

function exact(p,x)
exact = (exp(p*x)-1)/(exp(p)-1)
end

このプログラムには、中心差分よりも高精度、しかも、安定に計算できる、パタンカーの「べき乗法」による計算もできるようにしている。

 

 


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

対流と拡散 第1回 [数値解析]

対流と拡散 第1回

 

パタンカーのコントロールボリューム法(有限体積法の一種)に基づく輸送方程式の離散化について、何回か、話をすることにする。

 

対流・拡散を有する対流・拡散を伴う輸送現象の基礎方程式は、次の通り。

  taikaku-001.png

ここで、ρは密度、uは流体の速度、φは温度や濃度、エンタルピーなどの物理量で、Γは拡散係数。

1次元ならばn=1、2次元ならばn=2、3次元ならばn=3

(1)は質量保存則を表しており、連続の式と呼ばれる。

(2)式が物理量φの輸送方程式で、左辺第2項は対流項、右辺第1項は拡散項と呼ばれる。右辺第2項のSφの生成・消滅を表す。

 

輸送問題で最も単純で、基礎になる定常1次元の対流と拡散について考える。

このとき、(1)、(2)式は、それぞれ、

   

(3)式から

  

 

計算に使用する点(格子)の配置は次の通り。

そして、はそれぞれwWPEにおけるφの値であり、などは検査体積の境界面wと境界面eにおけるφなどの値を表すものとする。

tai-kaku-fig-001.png

 

(4)式を検査体積(コントロール・ボリューム)で積分すると、

  

差分法を用いて、

  

さらに、境界面eと境界面wにおけるφの値を

  

で近似すると、(5)式は

  

ここで、

  

とおくとことにする。

すると、(6)式は

  

ここで、

  

 

F=ρuだから、(3)式は

  

これを検査体積で積分すると、

  

よって、ならば、(11)式は

  

 

離散化方程式(8)〜(11)は、対流項に中心差分を用いて(4)式を離散化することによって得る離散化方程式と同等のものなので、以降、これを中心差分と呼ぶことにする。

 

さて、の場合について考える。

の値が与えられると、(8)式によっての値を求めることができる。

そこで、

  

の場合について計算してみる。

(9)〜(11)式より、

  

だから、(8)式より

  

となり、物理的に非現実的な答になってしまう。

 ――この問題では、計算領域内でφの生成・消滅がないので、の間になければおかしい――

こうした事態が発生するとき、離散化方程式の係数の少なくとも一つが負になっている。

F=ρu>0のとき、係数を負にするものは

  

なので、

  

でなければならない。

この時のペクレ数Pを計算すると、

  

このように、この計算方法は、ペクレ数P(格子レイノルズ数)が

  

のとき、物理的にあり得ない数値解を出す危険性があることが知られている。

たとえば、ねこ騙し数学の次の記事。

バーガース方程式の数値解法6 中心差分

 

 

こうした中心差分の欠点を補うために、次のように風上化が行われる。

  taikaku-002.png

ここで、

  

という記号をを導入すると、

  

で表すことができ、これを(5)の左辺に代入し離散化方程式を書き換えると、次のようになる。


   


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

ロジスティック方程式の数値解法2 [数値解析]

ロジスティック方程式の数値解法2

 

一般のロジスティック方程式は

  

である。

これから、

  

という差分方程式が得られる。

ここで、

  

とおくと、差分方程式(2)は

  

になる。

  

とおくと、(3)式は

  

になる。

 

さて、a>0f(x)の定義域を[0,1]とすると、

  

だから、f(x)x=1/2のときに最大値a/4を取ることがわかる。

したがって、0≦a≦1/4ならば、f[0,1]から[0,1]への写像になる。

これをロジスティック写像と呼ぶ。

ロジスティック写像の不動点をxとすると、

  

となるので、不動点はである。

 

x₀=01のときは、漸化式(3)より、

  

そこで、0<x₀<1とすると、0<a≦1のとき、n≧1に対して

  

になる。

 

数値計算の話なので、結論だけを言うと、

漸化式(3)で得られるは、n→∞のとき、

  

に収束し、3<a≦4のとき、収束しないことが知られている。

 

論より証拠ということで、x₀=0.6a=0.9として、表計算ソフトで計算した結果は次の通り。


x₀=0.6a=2のときの計算結果は次の通り。

 

0<a≦3のときの分岐図を以下に示す。横軸にはa、縦軸にはの極限値をとっている。

 

この図を見ると、a=1で極限値x0から1−1/aに分岐(?)していることがわかるだろう。

――a=3付近での収束は非常に遅い!!――

 

2≦a≦4の場合の分岐図は、次の通り。

3<a≦1+√6では、ある2つの値を交互に取るような振る舞いをする。

さらに、1+√6≦a≦3.57では、aの増加に伴い、4、8、16・・・と分岐してゆく。

 

そして、約3.57<a≦4ではカオス的になることが知られている。


ついでに、書いただけだにゃ。

ネムネコは、ただのネコだから、カオス理論、フラクタルや力学的不安定といった話は知らないケロ。

だから、これ以上の話を知りたいヒトは、カオス理論の本を買って勉強するなり、ウィキペディアの次の記事などを読むといいと思うにゃ。

 

https://goo.gl/mBzkEm

 

(非線形方程式の)力学的不安定の話に関しては、きっと、ddt³さんがしてくれるんじゃないか(^^

期待しているケロ!!

 


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