Python で1次元の熱方程式を2通りの方法で解いてみた

両端が断熱的な導線の熱伝導を記述する, 次の斉次 Neumann 境界付きの熱方程式の初期値境界値問題

\begin{equation}
\left\{
\begin{split}
&\partial_t u = \partial_x^2 u,\quad &t>0 , x\in (0,\pi),\\
&\partial_x  u(t,0) =\partial_x u(t,\pi) = 0 ,\quad& t>0,\\
&u(0,x) = f(x),\quad & x\in (0,\pi)
\end{split}
\right.
\end{equation}
を考えます*1.
ここで,  \partial_t  =\frac{\partial}{ \partial t},  \partial_x  =\frac{\partial}{ \partial x} を表します.
前進差分法とスペクトル法を用いる方法の2通りでの数値計算例のデモを作成しました. *2
コードがメインで, 数値解析の手法については多くを説明しません.
意味が通らない箇所があれば気軽にコメントをよろしくお願いします.

*1: u の第1引数を時間変数, 第2引数を空間変数としています.

*2:有限要素法でもやりたかった

続きを読む

不偏標本分散の一致性を示してみる

本記事では省略されがちであったり, 正規分布であることを仮定してあったりすることの多い, 不偏標本分散の母分散の推定量としての一致性を証明していきます.

続きを読む

Rにたくさんさいころを振らせて、さいころの分散を求めてみる

前回記事Rにさいころをたくさん振らせて、さいころの平均を求めてみる - わかばめにっきでは、Rの擬似乱数を使って標本平均が母平均の一致推定量であることを示しました。今回は後で述べる不偏標本分散という量が母分散の不偏推定量であることを調べていきます。

続きを読む

Rにさいころをたくさん振らせて、さいころの平均を求めてみる

最近は東京大学出版会の統計学入門を用いて、統計学に再入門しています。標本分散の定義に  n-1 が出てくるのはなぜだっけ?というところからよくわかっていなかったので、何回かに分けてその理由をメモしていきます。
推定とはなにかというところから出発します。若干迂遠に感じられるかもしれませんが、間違い探しだと思ってお付き合いください。

      • -

確率的に変動する要素がある場合、わたしたちはその変動をできるだけ正確に見積もりたいのに対して、厳密に調べつくすことは現実的ではないという問題に遭遇します。
例えばさいころを振ることを考えても、実際にはわたしたちは有限回しかさいころを振るという行動はできません。
したがって、「理想的なさいころであればそれぞれの目が出る確率は  \frac16 になるはずだ」ということに対して、手元にあるさいころの背後にある確率分布を精密に調べることは不可能です。
そうであれば、実験、調査、測定などによって得られたデータ(これを標本と呼ぶ)のどのような性質を見れば、確率分布の持つ平均分散などといった特徴(これを母数と呼ぶ)を調べられるでしょうか?また、性質の見方には指標があるでしょうか?ということを考える必要が生じます。
この記事群では一致推定量不偏推定量という観点から、母平均と母分散を"推定"することを考えていきます。*1

*1:実は個人的に計算でつまづいたところがあったので、自分用メモの側面もあります。

続きを読む

Python で波動方程式の数値計算と動画 gif の書き出しをやらせてみよう

2次元正方形領域における波動方程式の初期値境界値問題を考えます。
こんな感じで波が伝播していく様子を気軽にシミュレーションできます。
f:id:wakabame:20180307201626g:plain

続きを読む

プログラミングコンテストにおける Python での標準入出力のまとめ

入力

1行に1要素の場合

N = int(input())
  • 文字列
S = input()

beta 版のコードテスト環境だと改行を含む場合もあるので注意しよう。
無印のコードテストとテスト時には改行は含まれない。
f:id:wakabame:20180304121614p:plain
気になる場合は input().rstrip() とすればよい。

1行にスペース区切りで複数要素の場合

  • 2つの数字を別の変数 N, M に
N, M = map(int, input().split())
  • 2つの文字列を別の変数 S, T に
S, T = map(str, input().split())
  • 複数の数字をリスト A に
A = list(map(int, input().split()))

複数行に1要素がある場合

  • 複数の数字をリスト A に
A = [int(input()) for _ in range(N)]
  • 複数の文字列をリスト A に
A = [input() for _ in range(N)]

beta版のコードテスト環境下では input().rstrip() とすればよい。

複数行にスペース区切りで複数要素の場合

  • 複数の数字をリスト A に
A = [list(map(int, input().split())) for _ in range(N)]

A[i][j] は  i+1 行目の j+1 文字目を表す

出力

1要素の出力

  • 文字列または数字 A を改行つきで表示
print(A)

複数要素の出力

  • リスト A を1要素ごとに改行して出力
[print(A[i]) for i in range(N)]
  • リスト A を空白区切りで出力
print(' '.join(map(str, B)))

クォーテーション内の半角スペースが挿入される。
別の文字(例えばカンマ)で区切りたい場合はそこを変更すればよい。

ABC088 を Julia で解いてみる

先ほどまで開催されていた AtCoder の ABC088 に参加してみました。
問題や解説についてはリンク先をご覧ください。
abc088.contest.atcoder.jp
最近は数値計算をしてみたいなと思うところがあり、Julia の勉強を始めていました。
というわけで、使い慣れてきた Python と Julia の二つで解いていきました。

続きを読む