So-net無料ブログ作成
ひとこと言わねば ブログトップ
前の10件 | -

予定を少し変更!! [ひとこと言わねば]

今日、ddt³さんから、興味深い記事を頂いたので、当初の予定を変更し、頂いたddt³さんの記事を紹介します。
明日、掲載分の数学の記事を書いていなかったので、助かったにゃ。まさに、天佑だケロよ!!


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

ネムネコ、独り言をいう・・・ [ひとこと言わねば]

あっ、そうか。陽解法の方がクランク・ニコルソン法よりも精度が良くなるように思えるのは、厳密解である

この(フーリエ)級数の計算に問題があるのか。

このプログラムでは、n=110までをとって厳密解としている。

このため、例えば、t=0のとき、x=1で3、x=24にならず、小数点以下4位で打切誤差が入るんだケロ。

n=30くらいにすると、3.000013.99999になる。

こうすると、少し違った印象をもったのかもしれない。

 

そうは言っても、フーリエ級数の打切誤差は1/10000程度だからね〜。グラフにしたとき、この差がグラフに反映されるとも思えない。

そうだとすると、理由は不明ながら、Δt=0.1Δx=0.4とし、例の問題を解いたとき、陽解法がもっとも精度よく計算できるということになるのだろうか・・・。




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

ネムネコ、さらに混乱する!! [ひとこと言わねば]

1次元の熱伝導方程式の数値解法の最後の記事として、ルンゲ・クッタ法を用いた解法を紹介しようと、プログラムを作ってみた。
計算は発散するし、精度は最も簡単な陽解法より悪いわと散々な結果だったケロ。

ルンゲ・クッタの計算段数を増やせば増やすほど、精度が悪くなるんだにゃ。1次の時、つまり、オイラー法で陽解法より、2次のルンゲ・クッタの精度が悪く、一般的に使われる4次のルンゲ・クッタはさらに精度が悪化する。不思議だケロね〜。
オレ、計算法を間違っているんかな〜。

どこが悪いのか、何故、精度が悪化するのか、理由はわからないのだけれど、せっかくプログラムを作ったのだから、明後日、数学の記事として、作ったプログラムを公開するにゃ。
だから、どこか間違っているところがあったら、教えて欲しいにゃ。


熱や流体力学の数値計算の本には、これらの計算でルンゲ・クッタが使われることがあるみたいなことが少し書いてあるんだけれどね。


慣れないことは、するもんじゃないケロね〜!!

プログラム、1箇所間違っていたケロ(^^ゞ
式が長かったから、間違いに気づかなかったにゃ。
安定性は、陽解法よりもいくらか緩和されるようだけれど、思ったほどの精度は出ないね。

でも、これでまともな記事を書ける。良かったにゃ。

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

当惑するネムネコ [ひとこと言わねば]

当惑するネムネコ

 

明後日23日の記事のために

  

を、次の条件で数値的に解いてみた。

  

ちなみに、κ=1/2

 

Δt=0.1Δx=0.4として、陽解法、陰解法、明日やる予定のクランク・ニコルソン法で解いてみたんだケロ。

そうしたら、最も精度が良いはずのクランク・ニコルソン法よりも単純な陽解法の方が精度よく計算できてしまったケロ。

経験上、陽解法の方が(純)陰解法よりも厳密解に近い近似解を求めることは知っていた。

しかし、この問題の場合、何故、理論上最も高精度のクランク・ニコルソン法より陽解法のほうが精度よく計算できたか、その理由がわからないにゃ。

 

Explicit.png

Crank-Nicolson method.png

上のグラフが陽解法、下のグラフがクランク・ニコルソン法。その差は微妙だけれど、心持ち陽解法のほうが精度よく計算できている。

連立方程式を解く過程で、丸め誤差や桁落ち、情報落ちなどが起きてこのようなことが起きた(・・?

だから、倍精度で計算すると、結果が変わる(・・?

 

予想外の結果にただただ当惑するネムネコでありました。




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

ネムネコ、プログラムで大きな間違いをおかす!! [ひとこと言わねば]

ずっと前に、ブログで紹介した非定常の1次元熱伝導方程式を解くプログラムに大きな間違いがあったにゃ。
なので、コッソリと問題箇所を修正する(^^ゞ

境界条件が0だから、この間違いがあったとしても、計算結果は正しいのだけれど、境界条件が0じゃないと、修正前のプログラムはおかしな答えを出す。これは係数の符号の正負を勘違いしていたために起きたミス。たまには、こういううっかりミスもおかすにゃ。

紙と鉛筆をいっさい使わず、勘と見込み(思い込みか?)だけでプログラムをつくり、しかも、1つの例で正しい答えが出れば、このプログラムには間違いはないと考えれるのだから、間違えることもあるわな。アタリマエのことなのだけれど、もう少し真面目にプログラムのテストはしないといけないね〜。

この間違ったプログラムとは別に新たにプログラムを作ったにゃ。コチラは、境界条件が0でない場合でもちゃんと計算してくれる。明日、数学の記事と共に、このプログラムを公開するにゃ。
この記事を読むと、何故、このプログラムで計算できるのか、その仕組み、そして、プログラムで何をやっているかが分かると思うにゃ。
さらに、明後日には、別の解法を用いた、精度の良い計算法とプログラムを公開するにゃ。

この2つのプログラムは間違っていないと思うけれど、ネムネコが恥を晒し続けないために、間違っているところがあったら、ネムネコに教えて欲しいにゃ。



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

「陰的オイラー法」の記事で使用したスプレッドシート [ひとこと言わねば]

例によりまして、「陰的オイラー法」という記事で用いたスプレッドシート公開のお知らせ。



上記リンクをクリックすると、スプレッドシートがダウンロードされるにゃ。

「ネムネコはこのスプレッドシートに何かウィルスを仕込んでいるに違いない!危ないケロ」というヒトは、ここ↓を見るといいにゃ。
https://docs.google.com/spreadsheets/d/e/2PACX-1vR0UKslY05jzOGMZeMK2THHGxlEGN9-PzJirmX6LVfSi3DhBPU0Le-Y7uApx3j3Qq2lH9XfSmvXuM3K/pubhtml

計算はしてくれないけれど、大きな画面で計算結果とグラフを見ることができるケロよ。

これも前回のスプレッドシート同様に、λ、hのセル、すなわち、I2セルとJ2セルの値を変えると、この値で自動的に計算し、その結果をグラフにしてくれるにゃ。ただ、Googleの表計算ソフトを一回経由しているので、グラフはオリジナルのものと違うものになる。だから、使っている表計算ソフトのグラフ機能を用いて、グラフを作ったほうがいいと思うにゃ。



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

ネムネコ、たまげる! そして、自画自賛する!! [ひとこと言わねば]

「オイラー法の誤差の内訳」で書いた打切誤差が具体的に何を意味するのか分からず――「積分の理論的な誤差なのか、そうではなく、オイラー法を用いて計算した近似値に基づく積分の誤差なのか分からなかった」と書いても、何のことか、誰もわからないに違いない(^^ゞ――、電卓でちょこっと計算したみたケロ。

λ=1h=0.1のとき、

  

の厳密解だから、

  

そして、

  

だから、(1)で計算した誤差は、−0.012718406になる。

 

 

この表の打切誤差を見ると、ドンピシャ、この値が出ている。後追い的にこのことが確かめられて、驚いたケロよ。

この表の打切誤差は、相対誤差から伝播誤差を差し引いて求めたんだけどね。

あまりに理論通りなので、驚きだケロよ。


つまり、「オイラー法の誤差の内訳」におけるネムネコの計算は、理論的なところを含めて、全て正しかったということになるにゃ。

嘘じゃないにゃ。


と、相対誤差と伝播誤差のさから打切誤差を求めているにゃ。この計算に使用したスプレッドシートを公開し、誰でもダウンロード可能でこのことを確かめられるようになっているから、嘘はつけないにゃ。


さすが、ネムネコだにゃ!!






今回は、まじまじと見せつけたに違いないにゃ。

打切誤差と伝播誤差についてここまで詳しく、そして、具体的な数値まで出して紹介しているブログ、ネットの記事(日本語版)はないと思うにゃ。
したがって、あの記事を読めたヒトは幸せだと思うにゃ。
だって、どこにも出てねえんだから!!


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

「オイラー法の誤差の内訳」で使用したスプレッドシートを公開したにゃ [ひとこと言わねば]

「オイラー法の誤差の内訳」で使用したスプレッドシートは、
エクセル版
https://docs.google.com/spreadsheets/d/e/2PACX-1vSXCnqWUxyG_3iwBoC_Kg2DkIb9SiRoIr2wxcg7oX8jLUXZVgYC7r-8csKkJu-w19NQngJYhoUgFZUw/pub?output=xlsx
LibreOffice Calc版
https://docs.google.com/spreadsheets/d/e/2PACX-1vSXCnqWUxyG_3iwBoC_Kg2DkIb9SiRoIr2wxcg7oX8jLUXZVgYC7r-8csKkJu-w19NQngJYhoUgFZUw/pub?output=ods
に公開してあるにゃ。
物好きなヒトは、上記アドレスをクリックすると、このスプレッドシートをダウンロードすることができるケロ。

G2セルと、F2セルの値を変えると、それに応じて、自動的に計算してくれるというすぐれものだにゃ。しかも、グラフまで書いてくれるにゃ。
とは言え、Googleさんのスプレッドシートを経由なので、グラフには難ありなので、グラフの部分はスプレッドシートをダウンロードしたあと、使っている表計算ソフトのグラフ機能を使ったほうがいいにゃ。

「ネムネコはこのスプレッドシートに何かマクロを仕込んでいるかもしれない。危ないケロ」と思うヒトは、
https://docs.google.com/spreadsheets/d/e/2PACX-1vSXCnqWUxyG_3iwBoC_Kg2DkIb9SiRoIr2wxcg7oX8jLUXZVgYC7r-8csKkJu-w19NQngJYhoUgFZUw/pubhtml

を見るといいにゃ。λ=−1、h=0.1のときの計算結果が出ているにゃ。


「ねこ騙し数学」は、計算に使用したプログラム、スプレッドシートは、原則として、すべて公開するにゃ。何も隠さないというのがこのブログの方針だにゃ。


自画自賛、自尊だけじゃないケロ、このブログは。ネムネコは常に訪問者の便宜と幸せを考えているにゃ。



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

ひょっとして、平行平板ポアズイユ流れの熱伝達のプログラム(バージョンアップ)が欲しいケロか? [ひとこと言わねば]

今日の未明だかに書いた「平行平板ポアズイユ流れの熱伝達」に関するアクセスが思いの外、多いケロ。

ひょっとして、平行平板ポアズイユ流れの熱伝達のプログラム(バージョンアップ)が欲しいケロか?

しょうがないにゃ。計算に使用したプログラムを紹介するにゃ。


! U∂T/∂X=Γ(∂²T/∂X²+∂²T/∂Y²)を解くプログラム
parameter (m=10,n=5)
real T(0:m,0:n)
real a(5,m,n)
real u(n), tb(m),dnu(m)

limit = 1000 ! 反復回数の上限

! 温度Tの初期化
do i=0,m
    do j=0,n
        t(i,j)=0.
    end do
end do

! 壁の温度T=1を設定
do i=1,m
    t(i,0)=1.
end do

! ΔX、ΔYを計算
dx=10./m ; dy = 0.5/n
dx2=dx*dx; dy2=dy*dy

! 計算条件 Pr=0.7、Re=100
pr=0.7 ; re=100
gam=1./(pr*re)

! 流速Uの設定
do j = 1,n
    y=dy*j
    u(j)=6.0*y*(1-y)
end do

! 連立方程式の係数を計算
do i=1,m
    do j=1,n
        a(1,i,j)=gam/dx2+u(j)/dx
        a(2,i,j)=gam/dx2
        a(3,i,j)=gam/dy2
        a(4,i,j)=gam/dy2
    end do
end do

! 出口の条件d²T/dX²=0を係数に反映
do j=1,n
    a(1,m,j)=u(j)/dx
    a(2,m,j)=0.
end do

! 格子点(i,j)におけるT(i,j)の係数を計算
do i=1,m
    do j=1,n
        a(5,i,j)=a(1,i,j)+a(2,i,j)+a(3,i,j)+a(4,i,j)
    end do
end do

! ガウス・ザイデル反復法で連立方程式を解く
do k=1,limit
    res_max = 0.
    do i=1,m-1
        do j=1,n
            if (j.ne.n) then
                res =a(1,i,j)*t(i-1,j)+a(2,i,j)*t(i+1,j)+a(3,i,j)*t(i,j-1)+a(4,i,j)*t(i,j+1) &
                 -a(5,i,j)*t(i,j)
            else
                res =a(1,i,n)*t(i-1,n)+a(2,i,n)*t(i+1,n)+2*a(3,i,n)*t(i,n-1) &
                 -a(5,i,j)*t(i,j)
            end if
            res_max=amax1(abs(res),res_max)
            t(i,j)=t(i,j)+res/a(5,i,j)
        end do
    end do
    do j=1,n
        if (j.ne.n) then
            res = a(1,m,j)*t(m-1,j)+a(3,m,j)*t(m,j-1)+a(4,m,j)*t(m,j+1) - a(5,m,j)*t(m,j)
        else
            res =a(1,m,n)*t(m-1,n)+2*a(3,m,n)*t(m,n-1) - a(5,m,n)*t(m,j)
        end if
        res_max=amax1(abs(res),res_max)
        t(m,j)=t(m,j)+res/a(5,m,j)
    end do
    if (res_max .lt. 1.0e-5) exit ! 収束条件を満たせば反復計算を抜け出す
end do

! 混合平均温度 Tb=∫UTdy/∫udYを計算
do i=1,m
    sum = 0.5*u(n)*t(i,n)
    do j=1,n-1
        sum = sum + u(j)*t(i,j)
    end do
    tb(i)=2*sum*dy
end do

write(6,*) '*** result***'
write(6,110) 'x  / y  ',(dy*j,j=1,n),'Tb', 'Nu'
do i=1,m
    dtdy= (4*t(i,1)-t(i,2)-3*t(i,0))/(2*dy) ! 壁の温度勾配
    dnu(i) = -dtdy/(1.0-tb(i)) ! ヌセルト数Nuを計算
    write(6,100) dx*i,(t(i,j),j=1,n),tb(i),dnu(i)
end do

100 format(7(f8.5,1x),f8.5)
110 format(a8,1x,5(f8.5,1x),a8,1x,a8)
end

収束を早めるために、ガウス・ザイデルじゃなくSOR法を使おうとしたら、何故か、計算が発散したり、収束がガウス・ザイデル法よりも遅くなったりと、予期せぬ事態に出食わした。でもまぁ〜、それっぽい答えが出ているから、たぶん、プログラムは間違っていないと思うにゃ。




紹介したFortranプログラムをコンパイルし、実行すると、上の計算結果が得られるはずだにゃ。そして、m=100、n=10として、ヌセルト数Nu(熱伝達係数を無次元化したもの)を計算したものが、先の記事のこのグラフだにゃ。



さらに、この動画を埋め込んでおこう!!



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

平行平板ポアズイユ流れの熱伝達がさらに進化 [ひとこと言わねば]

平行平板ポアズイユ流れの熱伝達の方程式に

X方向の熱伝導の項を考慮した

 

にして計算してみたにゃ。

ここで、

で、(2)の境界条件は以下の通り。

 

を加えても、流入口X=0付近の値が変わるだけで、(1)、(2)で解こうがあまり変わらないということは知っていたのだけれど、実際に計算してみて、予想以上にこの差が小さくてショックを受けたケロ。

 

(2)を解くためのプログラムを作って、損こいたにゃ。

 

ちなみに、Pr=0.7Re=100での計算結果は次の通り。

 

 

Parabolicと書いてあるのが(1)、Ellipticと書いてあるのが(2)だにゃ。

Nuは、熱伝達係数を無次元化したもので、熱伝達係数と考えていい。


nice!(2)  コメント(0) 
前の10件 | - ひとこと言わねば ブログトップ