FORTRANプログラムのリファクタリングをしている。そんなことをするからにはデータフローもなんだかよくわからない悲惨なプログラムである。
で、ロジックを推測して論理的に書きなおすのであるが、プログラムを壊してしまっていないことを確認しながらバージョン管理をすすめていかないと険呑このうえない。 そこでアサーション（想定していない動作をしたら落ちるような確認）をほうぼうに入れていくのだが、なにせFORTRANだからその対象はきまって巨大な配列である。1行で確認するためには、定数データをプログラムに埋め込むわけにいかないから、外部にファイルで持つべきである。
で、こんなサブルーチンを書いてみた。
2013-01-31
2013-01-30
いろいろURLメモ（銚子積雪事例、CBS-MG14, RA-II, media type suffix）
ぜんぜんかんけいないのだがご容赦。
銚子積雪事例のＮＣＥＰ再解析データによる図など
http://taket.at.webry.info/201301/article_21.html
デイリーポータルZ：街路樹に名札がついている国、日本
http://portal.nifty.com/kiji/130125159288_1.htm
（桜が枯れなくなったのは、百年で地下水位が下がったということではないのか）
WMO/CBS-MG-14
http://www.wmo.int/pages/prog/www/WIS/wiswiki/tiki-index.php?page=CBS-MG-14
WMO/RA-II-15
http://raii-15.wmo.int/documents-english
RFC6838 "Media Type Specifications and Registration Procedures"
http://tools.ietf.org/html/rfc6838
メディアタイプのサフィックス +xml が一般化されて +zip +json とかが使えるようにサフィックスレジストリを立てるという。
初期登録設定値は RFC 6839 http://tools.ietf.org/html/rfc6839
Java の jar ファイルは application/java-archive+zip であるべきだったということになろう。
ちなみに、zip ファイルの内容物の URI スキームは jar: という
https://www.iana.org/assignments/uri-schemes/prov/jar
ところで +fastinfoset てのが ITU-T X.891 なんだが、EXI とは別のバイナリXMLなんだそうだ。
http://en.wikipedia.org/wiki/Fast_Infoset
銚子積雪事例のＮＣＥＰ再解析データによる図など
http://taket.at.webry.info/201301/article_21.html
デイリーポータルZ：街路樹に名札がついている国、日本
http://portal.nifty.com/kiji/130125159288_1.htm
（桜が枯れなくなったのは、百年で地下水位が下がったということではないのか）
WMO/CBS-MG-14
http://www.wmo.int/pages/prog/www/WIS/wiswiki/tiki-index.php?page=CBS-MG-14
WMO/RA-II-15
http://raii-15.wmo.int/documents-english
RFC6838 "Media Type Specifications and Registration Procedures"
http://tools.ietf.org/html/rfc6838
メディアタイプのサフィックス +xml が一般化されて +zip +json とかが使えるようにサフィックスレジストリを立てるという。
初期登録設定値は RFC 6839 http://tools.ietf.org/html/rfc6839
Java の jar ファイルは application/java-archive+zip であるべきだったということになろう。
ちなみに、zip ファイルの内容物の URI スキームは jar: という
https://www.iana.org/assignments/uri-schemes/prov/jar
ところで +fastinfoset てのが ITU-T X.891 なんだが、EXI とは別のバイナリXMLなんだそうだ。
http://en.wikipedia.org/wiki/Fast_Infoset
2013-01-22
豆：WebサーバのSSL証明書の有効期間を表示するコマンド
echo 'HEAD / HTTP/1.0' | \
openssl s_client -connect www.example.com:443 2>&1 | \
openssl x509 -text -noout | awk '/Validity/,/After/'
うまくいけばこんなふうに表示される。
Validity
Not Before: Jan 20 12:54:09 2013 GMT
Not After : Apr 23 11:57:43 2017 GMT
openssl s_client -connect www.example.com:443 2>&1 | \
openssl x509 -text -noout | awk '/Validity/,/After/'
うまくいけばこんなふうに表示される。
Validity
Not Before: Jan 20 12:54:09 2013 GMT
Not After : Apr 23 11:57:43 2017 GMT
2013-01-16
[気象庁防災XML] 流量調査
流量について話題になっていたのでちょっと調べてみた。
- 日通数は最小６００、最大１０００
- 日電文量は最小１０ＭＢ，最大２０ＭＢ
- 流量の多い日には通数に比例よりはやや多くなる（つまり荒天時の電文は通常より長い傾向があることを含意する）
2013-01-11
被積分関数をかえてみる
簡単な例なので被積分関数をいじってみる。dy/dx = f(x, y) = 2x というのは y に依存しないので簡単すぎたのだ。ならば、 f(x,y) = x + sqrt(y) とやっても解析的には同じことである
が、 y がずれると積分結果がずれていってしまうという意味でよい感じになる。
次の結果で上が４次ルンゲクッタ、下が前進差分。もちろんバイアスがないルンゲクッタで定数倍改善していることも注目に値するが、ステップを 1/10 にしたときに１桁以上精度が改善しているこ
とに注目される。
0.3000000000000000 0.9717281348583426 -0.0282718651416574
0.1000000000000000 0.9945397310971597 -0.0054602689028403
0.0300000000000000 0.9991023660366382 -0.0008976339633618
0.0100000000000000 0.9998272390531332 -0.0001727609468668
0.0030000000000000 0.9999716121300442 -0.0000283878699558
0.0010000000000000 0.9999945367413657 -0.0000054632586343
0.0003000000000000 0.9999991022947341 -0.0000008977052659
0.0001000000000000 0.9999998272366565 -0.0000001727633435
0.0000300000000000 0.9999999716130413 -0.0000000283869587
0.0000100000000000 0.9999999945402596 -0.0000000054597404
0.0000030000000000 0.9999999991117637 -0.0000000008882363
0.0000010000000000 0.9999999998104413 -0.0000000001895587
0.0000003000000000 1.0000000000950195 0.0000000000950195
0.0000001000000000 1.0000000004349120 0.0000000004349121
0.0000000300000000 0.9999999989783790 -0.0000000010216210
0.0000000100000000 0.9999999959029633 -0.0000000040970367
0.3000000000000000 0.5100000000000001 -0.4899999999999999
0.1000000000000000 0.8100000000000004 -0.1899999999999996
0.0300000000000000 0.9410999999999989 -0.0589000000000011
0.0100000000000000 0.9800999999999989 -0.0199000000000011
0.0030000000000000 0.9940109999999988 -0.0059890000000012
0.0010000000000000 0.9980009999999985 -0.0019990000000015
0.0003000000000000 0.9994001100000743 -0.0005998899999257
0.0001000000000000 0.9998000100001390 -0.0001999899998610
0.0000300000000000 0.9999400011009805 -0.0000599988990195
0.0000100000000000 0.9999800001035221 -0.0000199998964779
0.0000030000000000 0.9999940000204693 -0.0000059999795307
0.0000010000000000 0.9999979999842054 -0.0000020000157946
0.0000003000000000 0.9999994001235168 -0.0000005998764832
0.0000001000000000 0.9999998004403844 -0.0000001995596156
0.0000000300000000 0.9999999389792776 -0.0000000610207224
0.0000000100000000 0.9999999759031382 -0.0000000240968618
import java.text.DecimalFormat;
public strictfp class Z {
static double gradient(double x, double y) {
return x + Math.sqrt(y);
}
public static double rk4(double xmin, double xmax, double ymin,
double step) {
double dx = step;
double xstop = xmax - dx;
double inv6 = 1.0/6.0;
double k1, k2, k3, k4;
double x = xmin;
double y = ymin;
boolean continuep = true;
while (continuep) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
k1 = gradient(x, y);
k2 = gradient(x + 0.5 * dx, y + 0.5 * dx * k1);
k3 = gradient(x + 0.5 * dx, y + 0.5 * dx * k2);
k4 = gradient(x + dx, y + dx * k3);
y += dx * (k1 + 2.0 * (k2 + k3) + k4) * inv6;
x += dx;
}
return y;
}
public static double euler(double xmin, double xmax, double ymin,
double step) {
double dx = step;
double xstop = xmax - dx;
double dydx;
double x = xmin;
double y = ymin;
boolean continuep = true;
while (continuep) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
dydx = gradient(x, y);
y += dydx * dx;
x += dx;
}
return y;
}
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {3.0e-1, 1.0e-1, 3.0e-2, 1.0e-2, 3.0e-3, 1.0e-3, 3.0e-4,
1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
for (double step: dxs) {
double y = rk4(xmin, xmax, ymin, step);
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
System.out.println("");
for (double step: dxs) {
double y = euler(xmin, xmax, ymin, step);
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
が、 y がずれると積分結果がずれていってしまうという意味でよい感じになる。
次の結果で上が４次ルンゲクッタ、下が前進差分。もちろんバイアスがないルンゲクッタで定数倍改善していることも注目に値するが、ステップを 1/10 にしたときに１桁以上精度が改善しているこ
とに注目される。
0.3000000000000000 0.9717281348583426 -0.0282718651416574
0.1000000000000000 0.9945397310971597 -0.0054602689028403
0.0300000000000000 0.9991023660366382 -0.0008976339633618
0.0100000000000000 0.9998272390531332 -0.0001727609468668
0.0030000000000000 0.9999716121300442 -0.0000283878699558
0.0010000000000000 0.9999945367413657 -0.0000054632586343
0.0003000000000000 0.9999991022947341 -0.0000008977052659
0.0001000000000000 0.9999998272366565 -0.0000001727633435
0.0000300000000000 0.9999999716130413 -0.0000000283869587
0.0000100000000000 0.9999999945402596 -0.0000000054597404
0.0000030000000000 0.9999999991117637 -0.0000000008882363
0.0000010000000000 0.9999999998104413 -0.0000000001895587
0.0000003000000000 1.0000000000950195 0.0000000000950195
0.0000001000000000 1.0000000004349120 0.0000000004349121
0.0000000300000000 0.9999999989783790 -0.0000000010216210
0.0000000100000000 0.9999999959029633 -0.0000000040970367
0.3000000000000000 0.5100000000000001 -0.4899999999999999
0.1000000000000000 0.8100000000000004 -0.1899999999999996
0.0300000000000000 0.9410999999999989 -0.0589000000000011
0.0100000000000000 0.9800999999999989 -0.0199000000000011
0.0030000000000000 0.9940109999999988 -0.0059890000000012
0.0010000000000000 0.9980009999999985 -0.0019990000000015
0.0003000000000000 0.9994001100000743 -0.0005998899999257
0.0001000000000000 0.9998000100001390 -0.0001999899998610
0.0000300000000000 0.9999400011009805 -0.0000599988990195
0.0000100000000000 0.9999800001035221 -0.0000199998964779
0.0000030000000000 0.9999940000204693 -0.0000059999795307
0.0000010000000000 0.9999979999842054 -0.0000020000157946
0.0000003000000000 0.9999994001235168 -0.0000005998764832
0.0000001000000000 0.9999998004403844 -0.0000001995596156
0.0000000300000000 0.9999999389792776 -0.0000000610207224
0.0000000100000000 0.9999999759031382 -0.0000000240968618
import java.text.DecimalFormat;
public strictfp class Z {
static double gradient(double x, double y) {
return x + Math.sqrt(y);
}
public static double rk4(double xmin, double xmax, double ymin,
double step) {
double dx = step;
double xstop = xmax - dx;
double inv6 = 1.0/6.0;
double k1, k2, k3, k4;
double x = xmin;
double y = ymin;
boolean continuep = true;
while (continuep) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
k1 = gradient(x, y);
k2 = gradient(x + 0.5 * dx, y + 0.5 * dx * k1);
k3 = gradient(x + 0.5 * dx, y + 0.5 * dx * k2);
k4 = gradient(x + dx, y + dx * k3);
y += dx * (k1 + 2.0 * (k2 + k3) + k4) * inv6;
x += dx;
}
return y;
}
public static double euler(double xmin, double xmax, double ymin,
double step) {
double dx = step;
double xstop = xmax - dx;
double dydx;
double x = xmin;
double y = ymin;
boolean continuep = true;
while (continuep) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
dydx = gradient(x, y);
y += dydx * dx;
x += dx;
}
return y;
}
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {3.0e-1, 1.0e-1, 3.0e-2, 1.0e-2, 3.0e-3, 1.0e-3, 3.0e-4,
1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
for (double step: dxs) {
double y = rk4(xmin, xmax, ymin, step);
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
System.out.println("");
for (double step: dxs) {
double y = euler(xmin, xmax, ymin, step);
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
4次ルンゲクッタの威力（？）
さて、4次ルンゲクッタにしたところ、困ったことに誤差がなくなってしまった。
前進差分と違って、バイアスがないし、高次なものだから２次関数のシミュレートなんて厳密にできてしまう。刻みが 0.3 ですら解析解と全ビット同じ値をたたき出すのである。結果
として、ステップが小さくなるにつれて大きくなる情報落ちだけが誤差となる。
教材としては適当ではないが、情報落ちの可視化としては適当なのでメモしておく。
結果
ステップ 積分結果 誤差（対解析解）
0.3000000000000000 1.0000000000000000 0.0000000000000000
0.1000000000000000 1.0000000000000002 0.0000000000000002
0.0300000000000000 0.9999999999999991 -0.0000000000000009
0.0100000000000000 0.9999999999999990 -0.0000000000000010
0.0030000000000000 0.9999999999999990 -0.0000000000000010
0.0010000000000000 0.9999999999999988 -0.0000000000000012
0.0003000000000000 1.0000000000000695 0.0000000000000695
0.0001000000000000 1.0000000000001080 0.0000000000001079
0.0000300000000000 1.0000000000007920 0.0000000000007920
0.0000100000000000 1.0000000000031826 0.0000000000031826
0.0000030000000000 1.0000000000076747 0.0000000000076747
0.0000010000000000 0.9999999999833878 -0.0000000000166122
0.0000003000000000 1.0000000001106562 0.0000000001106562
0.0000001000000000 1.0000000003876117 0.0000000003876117
0.0000000300000000 0.9999999990878823 -0.0000000009121177
0.0000000100000000 0.9999999963603514 -0.0000000036396486
コード
import java.text.DecimalFormat;
public strictfp class Z {
static double gradient(double x, double y) {
return x * 2.0;
}
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {3.0e-1, 1.0e-1, 3.0e-2, 1.0e-2, 3.0e-3, 1.0e-3, 3.0e-4,
1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
int midlimit = 1;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
double x, y, k1, k2, k3, k4;
double inv6 = 1.0/6.0;
boolean continuep;
for (x = xmin, y = ymin, continuep = true; continuep; ) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
k1 = gradient(x, y);
k2 = gradient(x + 0.5 * dx, y + 0.5 * dx * k1);
k3 = gradient(x + 0.5 * dx, y + 0.5 * dx * k2);
k4 = gradient(x + dx, y + dx * k3);
y += dx * (k1 + 2.0 * (k2 + k3) + k4) * inv6;
x += dx;
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
前進差分と違って、バイアスがないし、高次なものだから２次関数のシミュレートなんて厳密にできてしまう。刻みが 0.3 ですら解析解と全ビット同じ値をたたき出すのである。結果
として、ステップが小さくなるにつれて大きくなる情報落ちだけが誤差となる。
教材としては適当ではないが、情報落ちの可視化としては適当なのでメモしておく。
結果
ステップ 積分結果 誤差（対解析解）
0.3000000000000000 1.0000000000000000 0.0000000000000000
0.1000000000000000 1.0000000000000002 0.0000000000000002
0.0300000000000000 0.9999999999999991 -0.0000000000000009
0.0100000000000000 0.9999999999999990 -0.0000000000000010
0.0030000000000000 0.9999999999999990 -0.0000000000000010
0.0010000000000000 0.9999999999999988 -0.0000000000000012
0.0003000000000000 1.0000000000000695 0.0000000000000695
0.0001000000000000 1.0000000000001080 0.0000000000001079
0.0000300000000000 1.0000000000007920 0.0000000000007920
0.0000100000000000 1.0000000000031826 0.0000000000031826
0.0000030000000000 1.0000000000076747 0.0000000000076747
0.0000010000000000 0.9999999999833878 -0.0000000000166122
0.0000003000000000 1.0000000001106562 0.0000000001106562
0.0000001000000000 1.0000000003876117 0.0000000003876117
0.0000000300000000 0.9999999990878823 -0.0000000009121177
0.0000000100000000 0.9999999963603514 -0.0000000036396486
コード
import java.text.DecimalFormat;
public strictfp class Z {
static double gradient(double x, double y) {
return x * 2.0;
}
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {3.0e-1, 1.0e-1, 3.0e-2, 1.0e-2, 3.0e-3, 1.0e-3, 3.0e-4,
1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
int midlimit = 1;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
double x, y, k1, k2, k3, k4;
double inv6 = 1.0/6.0;
boolean continuep;
for (x = xmin, y = ymin, continuep = true; continuep; ) {
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
k1 = gradient(x, y);
k2 = gradient(x + 0.5 * dx, y + 0.5 * dx * k1);
k3 = gradient(x + 0.5 * dx, y + 0.5 * dx * k2);
k4 = gradient(x + dx, y + dx * k3);
y += dx * (k1 + 2.0 * (k2 + k3) + k4) * inv6;
x += dx;
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
中間加算バッファによる情報落ち（ごめん桁落ちとは言わなかった）の軽減
さて、差分スキーム自体の話にいく前に、さっきの精度打ち止めの理由が情報落ちであることを確認しておく。
（用語訂正：さっきは桁落ちといったがふつうの用語ではなかった）。
情報落ちとは、指数部が大きく異なる浮動小数点数を加減算すると、増分の精度が数形式の精度より大きく劣ることである。
桁落ちとは、絶対値がほとんど同じで、なおかつ符号の反対な浮動小数点数を加算（符号が同じ数を減算といっても同じこと）すると、結果の相対精度が大きく劣ることで、ちょっと違う。
ステップが非常に小さい積分は、この情報落ちにぴったり当てはまる。
回避法は簡単で、積算時に中間積算を行って、あるていどまとまったところで積算変数への足しこみを行えばよい。次の例では1024回づつ中間積算することによって、ステップが 1e-9 でも精度打ち止めにならないよう
になった。（おそらくは精度を3桁改善できるだろうが、ものすごく時間がかかるので確認しない）
0.0001000000000000 0.9998999999999811 -0.0001000000000189
0.0000300000000000 0.9999700002000160 -0.0000299997999840
0.0000100000000000 0.9999900000000238 -0.0000099999999762
0.0000030000000000 0.9999970000020161 -0.0000029999979839
0.0000010000000000 0.9999989999999600 -0.0000010000000400
0.0000003000000000 0.9999996999999685 -0.0000003000000315
0.0000001000000000 0.9999998999999479 -0.0000001000000521
0.0000000300000000 0.9999999700006352 -0.0000000299993648
0.0000000100000000 0.9999999900010098 -0.0000000099989902
0.0000000030000000 0.9999999970078517 -0.0000000029921483
0.0000000010000000 0.9999999990116513 -0.0000000009883487
import java.text.DecimalFormat;
public strictfp class Z {
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8, 3.0e-9, 1.0e-9};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
int midlimit = 1024;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
double x, y, xmid, ymid, dydx;
boolean continuep;
int count;
for (x = xmin, y = ymin, xmid = 0.0, ymid = 0.0, continuep = true,
count = 0; continuep; ) {
dydx = (x + xmid) * 2.0;
if (x + xmid > xstop) {
dx = xmax - (x + xmid);
count = midlimit;
continuep = false;
}
ymid += dx * dydx;
xmid += dx;
if (++count > midlimit) {
x += xmid;
y += ymid;
xmid = ymid = 0.0;
count = 0;
}
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
（用語訂正：さっきは桁落ちといったがふつうの用語ではなかった）。
情報落ちとは、指数部が大きく異なる浮動小数点数を加減算すると、増分の精度が数形式の精度より大きく劣ることである。
桁落ちとは、絶対値がほとんど同じで、なおかつ符号の反対な浮動小数点数を加算（符号が同じ数を減算といっても同じこと）すると、結果の相対精度が大きく劣ることで、ちょっと違う。
ステップが非常に小さい積分は、この情報落ちにぴったり当てはまる。
回避法は簡単で、積算時に中間積算を行って、あるていどまとまったところで積算変数への足しこみを行えばよい。次の例では1024回づつ中間積算することによって、ステップが 1e-9 でも精度打ち止めにならないよう
になった。（おそらくは精度を3桁改善できるだろうが、ものすごく時間がかかるので確認しない）
0.0001000000000000 0.9998999999999811 -0.0001000000000189
0.0000300000000000 0.9999700002000160 -0.0000299997999840
0.0000100000000000 0.9999900000000238 -0.0000099999999762
0.0000030000000000 0.9999970000020161 -0.0000029999979839
0.0000010000000000 0.9999989999999600 -0.0000010000000400
0.0000003000000000 0.9999996999999685 -0.0000003000000315
0.0000001000000000 0.9999998999999479 -0.0000001000000521
0.0000000300000000 0.9999999700006352 -0.0000000299993648
0.0000000100000000 0.9999999900010098 -0.0000000099989902
0.0000000030000000 0.9999999970078517 -0.0000000029921483
0.0000000010000000 0.9999999990116513 -0.0000000009883487
import java.text.DecimalFormat;
public strictfp class Z {
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8, 3.0e-9, 1.0e-9};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
int midlimit = 1024;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
double x, y, xmid, ymid, dydx;
boolean continuep;
int count;
for (x = xmin, y = ymin, xmid = 0.0, ymid = 0.0, continuep = true,
count = 0; continuep; ) {
dydx = (x + xmid) * 2.0;
if (x + xmid > xstop) {
dx = xmax - (x + xmid);
count = midlimit;
continuep = false;
}
ymid += dx * dydx;
xmid += dx;
if (++count > midlimit) {
x += xmid;
y += ymid;
xmid = ymid = 0.0;
count = 0;
}
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
差分スキームの誤差と情報落ち [件名訂正]
差分スキームの誤差は O(Δx) のようにステップの冪であらわされる。つまり、アルゴリズムをなにも改良しなくても、ステップを小さくすれば誤差を小さくできる。しかしそれは無限精度の実数で計算ができる場合で、現実の浮動小数点演算においては別種の誤差もあるので無限高精度の積分を達成することはできない。ステップが小さいのに漫然と加算をしていれば桁落ち[訂正：情報落ち]がまず思いつく脅威である。
あまりいい例ではないかもしれないが、わかりやすいので、前進差分で dy/dx = 2x
を x: [0, 1] で積分してみた。もちろん解析解は ｙ = x^2 で答は 1 である。
出力をみると、ステップ 1e-8 程度で精度は頭打ちになっている。
import java.text.DecimalFormat;
public strictfp class Z {
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8, 3.0e-9, 1.0e-9};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
double x, y, dydx;
boolean continuep;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
for (x = xmin, y = ymin, continuep = true; continuep; x += dx) {
dydx = x * 2.0;
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
y += dx * dydx;
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
/* output
0.0001000000000000 0.9999000000001081 -0.0000999999998919
0.0000300000000000 0.9999700002007929 -0.0000299997992071
0.0000100000000000 0.9999900000031827 -0.0000099999968173
0.0000030000000000 0.9999970000096747 -0.0000029999903253
0.0000010000000000 0.9999989999833878 -0.0000010000166122
0.0000003000000000 0.9999997001106760 -0.0000002998893240
0.0000001000000000 0.9999999003876101 -0.0000000996123899
0.0000000300000000 0.9999999690878816 -0.0000000309121184
0.0000000100000000 0.9999999863603514 -0.0000000136396486
0.0000000030000000 0.9999999888525976 -0.0000000111474024
0.0000000010000000 1.0000000151388213 0.0000000151388213
*/
あまりいい例ではないかもしれないが、わかりやすいので、前進差分で dy/dx = 2x
を x: [0, 1] で積分してみた。もちろん解析解は ｙ = x^2 で答は 1 である。
出力をみると、ステップ 1e-8 程度で精度は頭打ちになっている。
import java.text.DecimalFormat;
public strictfp class Z {
public static void main(String[] args) {
DecimalFormat df = new DecimalFormat(" 0.0;-0.0");
df.setMinimumIntegerDigits(1);
df.setMaximumIntegerDigits(1);
df.setMinimumFractionDigits(16);
df.setMaximumFractionDigits(16);
double dxs[] = {1.0e-4, 3.0e-5, 1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7,
1.0e-7, 3.0e-8, 1.0e-8, 3.0e-9, 1.0e-9};
double xmin = 0.0;
double xmax = 1.0;
double ymin = 0.0;
double x, y, dydx;
boolean continuep;
for (double step: dxs) {
double dx = step;
double xstop = xmax - dx;
for (x = xmin, y = ymin, continuep = true; continuep; x += dx) {
dydx = x * 2.0;
if (x > xstop) {
dx = xmax - x;
continuep = false;
}
y += dx * dydx;
}
System.out.println(
df.format(step)
+ " " + df.format(y)
+ " " + df.format(y - (xmax * xmax - xmin * xmin))
);
}
}
}
/* output
0.0001000000000000 0.9999000000001081 -0.0000999999998919
0.0000300000000000 0.9999700002007929 -0.0000299997992071
0.0000100000000000 0.9999900000031827 -0.0000099999968173
0.0000030000000000 0.9999970000096747 -0.0000029999903253
0.0000010000000000 0.9999989999833878 -0.0000010000166122
0.0000003000000000 0.9999997001106760 -0.0000002998893240
0.0000001000000000 0.9999999003876101 -0.0000000996123899
0.0000000300000000 0.9999999690878816 -0.0000000309121184
0.0000000100000000 0.9999999863603514 -0.0000000136396486
0.0000000030000000 0.9999999888525976 -0.0000000111474024
0.0000000010000000 1.0000000151388213 0.0000000151388213
*/
2013-01-08
[気象庁防災XML] Title と InfoKind の関係
辞書やフォーマットの解説（ http://xml.kishou.go.jp/tec_material.html で得られる）を読む限り、データ型のサブクラスに相当する物は jmx_ib:InfoKind タグで指示されるようである。しか
し、全内容出力スタイルシートのカタログ（
http://xml.kishou.go.jp/jmaxml_20120615_Verification_StyleSheet_list.pdf ）では
InfoKind を用いずに jmx:Title で振り分けている。どういうことか、ここ２週間ほど入電した電文をデコードして得た (InfoKind, Title) 関係とつきあわせてみた：
・ある Title では同一 InfoKind しか出てこない。
・ある InfoKind は複数の Title で使われうる。
・同一 InfoKind であっても、Title によって異なる XSL を割り当てられている例がある。
・つまり、振り分けは jmx:Title によらざるを得ない。
InfoKind の使い道がわからないが、まあ整理学上のものというべきだろう。なぜこうなったかというと、XSLはデータ型というより所管課によく対応して作られていることはちょっと対象を知っている人ならばお見通しであろう。人だよ人。
InfoKind Title XSL
全般海上警報 全般海上警報（臨時） 120329_umikeiho1.xsl
全般海上警報 全般海上警報（定時） 120329_umikeiho1.xsl
台風解析・予報情報（３日予報） 台風解析・予報情報（３日予報）
120329_kfxc80.xsl
台風解析・予報情報（５日予報） 台風解析・予報情報（５日予報）
120329_kfxc80.xsl
同一現象用平文情報 府県気象情報 120329_ippanho.xsl
同一現象用平文情報 全般台風情報（定型） 120329_ippanho.xsl
同一現象用平文情報 地方気象情報 120329_ippanho.xsl
同一現象用平文情報 地方潮位情報 120329_choui1.xsl
同一現象用平文情報 全般気象情報 120329_ippanho.xsl
同一現象用平文情報 府県潮位情報 120329_choui1.xsl
同一現象用平文情報 全般潮位情報 120329_choui1.xsl
噴火に関する火山観測報 噴火に関する火山観測報 120615_kazan.xsl
地方海上予報 地方海上予報 120329_chihoumiyoho.xsl
地方海上警報 地方海上警報 120329_chihoumikeiho.xsl
地震情報 震源・震度に関する情報
120615_jishin_tsunami_tokai_decode_all.xsl
天候情報 全般天候情報 120329_TenkoJohoHtmlPlain.xsl
天候情報 府県天候情報 120329_TenkoJohoHtmlPlain.xsl
天候情報 地方天候情報 120329_TenkoJohoHtmlPlain.xsl
季節予報 全般１か月予報 120329_KisetuYohoHtmlPlain.xsl
季節予報 地方３か月予報 120329_KisetuYohoHtmlPlain.xsl
季節予報 全般３か月予報 120329_KisetuYohoHtmlPlain.xsl
季節予報 地方１か月予報 120329_KisetuYohoHtmlPlain.xsl
平文情報 府県天気概況 120329_ippanho.xsl
平文情報 全般週間天気予報 120329_ippanho.xsl
平文情報 地方週間天気予報 120329_ippanho.xsl
府県天気予報 府県天気予報 120329_ippanho.xsl
府県週間天気予報 府県週間天気予報 120329_shukan.xsl
東海地震関連情報 東海地震観測情報
120615_jishin_tsunami_tokai_decode_all.xsl
気象警報・注意報 気象警報・注意報 120329_kei_all_line_div.xsl
火山の状況に関する解説情報 火山の状況に関する解説情報
120615_kazan.xsl
特殊気象報 特殊気象報 120329_tokusyu.xsl
特殊気象報 季節観測 120329_tokusyu.xsl
環境気象情報 紫外線観測データ 120329_tokusyu.xsl
生物季節観測報告気象報 生物季節観測 120329_Sbt.xsl
異常天候早期警戒情報 異常天候早期警戒情報 120329_SoukeiHtmlPlain.xsl
竜巻注意情報 竜巻注意情報 120329_tatsumaki.xsl
震度速報 震度速報 120615_jishin_tsunami_tokai_decode_all.xsl
震源速報 震源に関する情報
120615_jishin_tsunami_tokai_decode_all.xsl
し、全内容出力スタイルシートのカタログ（
http://xml.kishou.go.jp/jmaxml_20120615_Verification_StyleSheet_list.pdf ）では
InfoKind を用いずに jmx:Title で振り分けている。どういうことか、ここ２週間ほど入電した電文をデコードして得た (InfoKind, Title) 関係とつきあわせてみた：
・ある Title では同一 InfoKind しか出てこない。
・ある InfoKind は複数の Title で使われうる。
・同一 InfoKind であっても、Title によって異なる XSL を割り当てられている例がある。
・つまり、振り分けは jmx:Title によらざるを得ない。
InfoKind の使い道がわからないが、まあ整理学上のものというべきだろう。なぜこうなったかというと、XSLはデータ型というより所管課によく対応して作られていることはちょっと対象を知っている人ならばお見通しであろう。人だよ人。
InfoKind Title XSL
全般海上警報 全般海上警報（臨時） 120329_umikeiho1.xsl
全般海上警報 全般海上警報（定時） 120329_umikeiho1.xsl
台風解析・予報情報（３日予報） 台風解析・予報情報（３日予報）
120329_kfxc80.xsl
台風解析・予報情報（５日予報） 台風解析・予報情報（５日予報）
120329_kfxc80.xsl
同一現象用平文情報 府県気象情報 120329_ippanho.xsl
同一現象用平文情報 全般台風情報（定型） 120329_ippanho.xsl
同一現象用平文情報 地方気象情報 120329_ippanho.xsl
同一現象用平文情報 地方潮位情報 120329_choui1.xsl
同一現象用平文情報 全般気象情報 120329_ippanho.xsl
同一現象用平文情報 府県潮位情報 120329_choui1.xsl
同一現象用平文情報 全般潮位情報 120329_choui1.xsl
噴火に関する火山観測報 噴火に関する火山観測報 120615_kazan.xsl
地方海上予報 地方海上予報 120329_chihoumiyoho.xsl
地方海上警報 地方海上警報 120329_chihoumikeiho.xsl
地震情報 震源・震度に関する情報
120615_jishin_tsunami_tokai_decode_all.xsl
天候情報 全般天候情報 120329_TenkoJohoHtmlPlain.xsl
天候情報 府県天候情報 120329_TenkoJohoHtmlPlain.xsl
天候情報 地方天候情報 120329_TenkoJohoHtmlPlain.xsl
季節予報 全般１か月予報 120329_KisetuYohoHtmlPlain.xsl
季節予報 地方３か月予報 120329_KisetuYohoHtmlPlain.xsl
季節予報 全般３か月予報 120329_KisetuYohoHtmlPlain.xsl
季節予報 地方１か月予報 120329_KisetuYohoHtmlPlain.xsl
平文情報 府県天気概況 120329_ippanho.xsl
平文情報 全般週間天気予報 120329_ippanho.xsl
平文情報 地方週間天気予報 120329_ippanho.xsl
府県天気予報 府県天気予報 120329_ippanho.xsl
府県週間天気予報 府県週間天気予報 120329_shukan.xsl
東海地震関連情報 東海地震観測情報
120615_jishin_tsunami_tokai_decode_all.xsl
気象警報・注意報 気象警報・注意報 120329_kei_all_line_div.xsl
火山の状況に関する解説情報 火山の状況に関する解説情報
120615_kazan.xsl
特殊気象報 特殊気象報 120329_tokusyu.xsl
特殊気象報 季節観測 120329_tokusyu.xsl
環境気象情報 紫外線観測データ 120329_tokusyu.xsl
生物季節観測報告気象報 生物季節観測 120329_Sbt.xsl
異常天候早期警戒情報 異常天候早期警戒情報 120329_SoukeiHtmlPlain.xsl
竜巻注意情報 竜巻注意情報 120329_tatsumaki.xsl
震度速報 震度速報 120615_jishin_tsunami_tokai_decode_all.xsl
震源速報 震源に関する情報
120615_jishin_tsunami_tokai_decode_all.xsl
登録:
投稿 (Atom)