第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」でコンパイルしています.
結果は次のようになりました.(全体をプロットしたものは最後にあります.)
python | julia | fortran | |
---|---|---|---|
sec [s] | 3774.5655 | 33.1704 | 50.7343 |
pythonではグローバル変数で行っているとはいえ,少しのチューニングで動的型付けであるjuliaは100倍以上早くなっている.
fortranは静的型付けであるためチューニングをしてないが,十分な速度がでている.
やはりpythonはデータ処理用で,数値解析には向かないみたいです.
というわけで,次回以降はますます格子が増え計算が重くなるので,juliaやfortranを使おうと思います.
よろしくお願いいたします.
ここまで見てくださり,誠にありがとうございます.
by hide
*おまけ
2 compariosnの条件でステップ数を増やした結果を添付しておきます.やはりビジュアルに強いやつが欲しかったからです.
Awesome post! Keep up the great work! 🙂
Great content! Super high-quality! Keep it up! 🙂
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.
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.
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.
Ahaa, its good conversation regarding this paragraph here at this webpage,
I have read all that, so at this time me also commenting here.
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
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!
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.
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!
As the admin of this web page is working, no question very rapidly it will be well-known, due to its feature contents.
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!
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.
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.
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!
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!!
Very neat article post.Really looking forward to read more. Cool.
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!
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!
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
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!