工学系大学院生のブログ

2人の大学院生による雑記ブログ

第7-5回 二次元拡散方程式[python, julia, fortran]



第7-5回では, 今までのまとめとして, 二次元拡散方程式をコードに落とし込んでいきます.



今回は「コードが回っているかの確認」の バリデーションと「python, julia fortranの差」の2つを行いました.




1 validation

$$\frac{\partial}{\partial t}(\rho C T)= \frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x} \right ) + \frac{\partial}{\partial y} \left ( \lambda \frac{\partial T}{\partial y} \right ) $$


計算条件として密度\(\rho\),比熱\(C\),熱伝導率\(\lambda\)一定条件を課します.


すると拡散係数\(D\)を用いて方程式は次のように書けます.


$$\frac{\partial T}{\partial t}=\frac{\lambda}{\rho C} \left( \frac{\partial T}{\partial x^2}+\frac{\partial T}{\partial y^2} \right)  $$

$$\frac{\partial T}{\partial t}=D \left( \frac{\partial T}{\partial x^2}+\frac{\partial T}{\partial y^2} \right)  $$




さて,下記の条件でpythonを用いて計算を行いました.


行った内容は「1次元の厳密解」「1次元の陽解法」「1次元の陰解法」「2次元の陽解法」「2次元の陰解法」の5つになります.


時間ステップ数,時間刻み幅はそれぞれの条件で行っています.


初期温度:\( T(x) = \left\{ \begin{array}{c} 2x \ \ \ \ \ (0\le x \le 1/2) \\ 2(1-x) \ \ \ \ \ (1/2 \le x \le 1) \end{array} \right.\)

境界条件:\( T(0,t)=T(L,t)=0\)

拡散係数:\(D=1.0\)

x方向のセル数:\(nx0=100\)

y方向のセル数:\(ny0=100 \)

x方向の刻み幅:\(dx=0.01\)

y方向の刻み幅:\(dy=0.01\)

⇒上記の条件から横幅は次の通り決まる.
\(L_x=1.0, L_y=1.0\)




上記の条件における1次元の厳密解は

$$T(x,t)=\sum ^{\inf}_{n=1} \frac{8D}{n^2 \pi ^2} sin(\frac{n \pi}{2}) sin(\frac{n \pi x}{L}) e^{(-D \frac{n^2 \pi ^2}{L^2}t)}$$


\(D=1.0, L=1.0\)より

$$T(x,t)=\sum ^{\inf}_{n=1} \frac{8}{n^2 \pi ^2} sin(\frac{n \pi}{2}) sin(n \pi x) e^{(-n^2 \pi ^2t)}$$

$$ \simeq \frac{8}{\pi ^2} \left[ sin(\pi x) e^{-\pi ^2 t}- \frac{1}{3^2} sin(3 \pi x) e^{-3^2 \pi ^2 t} + \frac{1}{5^2} sin(5 \pi x) e^{-5^2 \pi ^2 t} \right]$$


厳密解の算出は面倒なので,結果だけ書いておきました.



コードが多すぎるので下記Githubの[validation]内にコードを置いています.

https://github.com/hide-dog/2d-heat-equation



結果は次のようになりました.

2dでも計算が上手くいっているのではないのでしょうか.


今回は全ての条件において,同様の大きさの\(\Delta x, \Delta y \)を用いており,精度改善のため\(\Delta t\)をそれぞれ変更しています.


二次元の計算が目に見えるほどの誤差を持つのは\(\Delta x, \Delta y, \Delta t\)がまだ粗いからと考えています.気になる方がいましたら計算時間は増えますが,\(\Delta t\)を小さくして見て下さい.


簡単に誤差が小さくなります.







2 comparison

python, julia, fortranの比較を行います.



というのも今回のような計算では,pythonだとおもく,時間ステップ数を増やすとものすごく時間がかかったからです.





条件は次のように設定しました.

時間ステップ数:\(nt=1000\)

時間刻み幅:\(dt=0.01\)

初期温度:\(T(x,y,0)=293 (20℃)\)

境界条件: \( \left\{ \begin{array}{c} \frac{\partial T(x,0,t)}{\partial y}= \frac{\partial T(x,L_y,t)}{\partial y} =0\ \\ T(0,y,t)=573 (300℃) ,\ \ \ \ \ T(L_x,y,t)=293 (20℃) \end{array} \right.\)

アルミナの熱伝導:\(\lambda=32.0\)

アルミナの密度:\(\rho=3.94e3\)

アルミナの比熱:\(c=0.78e3\)

x方向のセル数:\(nx0=100\)

y方向のセル数:\(ny0=100 \)

x方向の刻み幅:\(dx=0.001\)

y方向の刻み幅:\(dy=0.001\)


行ったコードは「python」「julia」「fortran」の3つになります.



コードが多すぎるので下記Githubの[comparison]内にコードを置いてます.

https://github.com/hide-dog/2d-heat-equation



juliaはグローバル変数が極力減るようにチューニングしています.


fortranはコンパイルの際に最適化するようにオプション「-O4」でコンパイルしています.





結果は次のようになりました.(全体をプロットしたものは最後にあります.)

 pythonjuliafortran
sec [s]3774.565533.170450.7343


pythonではグローバル変数で行っているとはいえ,少しのチューニングで動的型付けであるjuliaは100倍以上早くなっている.


fortranは静的型付けであるためチューニングをしてないが,十分な速度がでている.



やはりpythonはデータ処理用で,数値解析には向かないみたいです.




というわけで,次回以降はますます格子が増え計算が重くなるので,juliaやfortranを使おうと思います.


よろしくお願いいたします.



ここまで見てくださり,誠にありがとうございます.


by hide



*おまけ


2 compariosnの条件でステップ数を増やした結果を添付しておきます.やはりビジュアルに強いやつが欲しかったからです.

このエントリーをはてなブックマークに追加

“第7-5回 二次元拡散方程式[python, julia, fortran]” への21件のフィードバック

  1. Awesome post! Keep up the great work! 🙂

  2. AffiliateLabz より:

    Great content! Super high-quality! Keep it up! 🙂

  3. Lasonya Mosher より:

    Hi there. I found your website by the use of Google even as searching for a comparable matter, your site got here up. It looks great. I have bookmarked it in my google bookmarks to come back then.

  4. Shelba Baynard より:

    Hey there. I discovered your blog via Google at the same time as searching for a similar subject, your site got here up. It seems to be great. I’ve bookmarked it in my google bookmarks to come back then.

  5. webhosting より:

    Hmm it seems like your website ate my first comment (it was extremely long) so I guess I’ll
    just sum it up what I submitted and say, I’m thoroughly enjoying your blog.
    I as well am an aspiring blog writer but I’m still new to
    the whole thing. Do you have any helpful hints for newbie blog writers?
    I’d certainly appreciate it.

  6. web hosting より:

    Ahaa, its good conversation regarding this paragraph here at this webpage,
    I have read all that, so at this time me also commenting here.

  7. What’s up all, here every person is sharing these experience, thus it’s fastidious
    to read this blog, and I used to go to see
    this website daily. adreamoftrains web hosting

  8. Cooper より:

    Hi! This is my first comment here so I just
    wanted to give a quick shout out and say I genuinely enjoy reading your articles.
    Can you suggest any other blogs/websites/forums that go over
    the same subjects? Many thanks!

  9. Rosia Lichty より:

    Hey there, You’ve performed a fantastic job. I will definitely digg it and personally suggest to my friends. I am sure they’ll be benefited from this site.

  10. Guy より:

    This is the perfect web site for anybody who really wants to understand this
    topic. You know a whole lot its almost hard to argue with you (not that I actually would want to…HaHa).
    You definitely put a new spin on a subject that’s been written about for years.

    Excellent stuff, just wonderful!

  11. Valorie より:

    As the admin of this web page is working, no question very rapidly it will be well-known, due to its feature contents.

  12. cheap flights より:

    It’s appropriate time to make some plans for the future and
    it is time to be happy. I’ve read this post and if I could I wish to suggest you few
    interesting things or tips. Perhaps you can write next articles referring to
    this article. I desire to read even more things about it!

  13. cheap flights より:

    I’ve been surfing on-line more than three hours lately, yet I by no means
    discovered any interesting article like yours. It’s beautiful value sufficient for me.
    Personally, if all webmasters and bloggers made excellent content
    as you did, the web will probably be a lot more helpful than ever before.

  14. cheap flights より:

    Great post. I was checking continuously this blog and I’m impressed!
    Extremely useful info specifically the ultimate part 🙂 I care for such info a lot.
    I used to be looking for this certain info for a long time.
    Thank you and good luck.

  15. Bettie より:

    Hi there! I could have sworn I’ve been to this web site before but after looking at some
    of the articles I realized it’s new to me. Nonetheless, I’m definitely happy I
    stumbled upon it and I’ll be bookmarking it and checking back often!

  16. Logan より:

    Hi, i think that i saw you visited my site thus
    i came to “return the favor”.I’m trying to find things
    to enhance my website!I suppose its ok to use a few of your ideas!!

  17. Very neat article post.Really looking forward to read more. Cool.

  18. Emilie より:

    Hey! I’m at work surfing around your blog from my new apple iphone!
    Just wanted to say I love reading through your blog and look forward to all your posts!
    Carry on the great work!

  19. Nannie より:

    Hi there, I found your web site by way of Google whilst looking for a
    comparable subject, your web site got here up, it appears to be
    like good. I’ve bookmarked it in my google bookmarks.
    Hi there, just become alert to your weblog thru Google, and found that it’s really informative.
    I am going to watch out for brussels. I will be grateful in the event you proceed this
    in future. Lots of people shall be benefited out
    of your writing. Cheers!

  20. 0mniartist より:

    It’s going to be ending of mine day, however before finish I am reading this wonderful paragraph to increase my know-how.

    0mniartist asmr

  21. graliontorile より:

    Hi there, just became alert to your blog through Google, and found that it’s really informative. I am gonna watch out for brussels. I will appreciate if you continue this in future. Many people will be benefited from your writing. Cheers!

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です