WScript.Echo(3 + 5) WScript.Echo(3 - 5) WScript.Echo(3 * 5) WScript.Echo(Math.pow(3, 5)) WScript.Echo(5 / 3) WScript.Echo(parseInt(5 / 3)) WScript.Echo(5 % 3) WScript.StdOut.Write(3 * 5 + "\n") WScript.StdOut.WriteLine(3 * 5)
var i = 3 * 5; WScript.Echo("3 * 5 = " + i);
for (var i = 1; i < 10; i++) { WScript.StdOut.Write(i + ", "); } WScript.StdOut.WriteLine();
for (var i = 1; i < 10; i++) { if (i % 3 == 0) { WScript.StdOut.Write(i + ", "); } } WScript.StdOut.WriteLine();
var sum = 0; for (var i = 1; i < 100; i++) { if (i % 3 == 0) { sum += i; } } WScript.Echo(sum);
Array.prototype.range = function(start, end) { var array = new Array(); var j = 0; for (var i = start; i <= end; i++) { array[j++] = i; } return array; } Array.prototype.map = function(func) { for (var i = 0; i < this.length; i++) { func(this[i]); } } var array = new Array(). range(1, 9). map(function (n) { WScript.Echo(n); });
Array.prototype.range = function(start, end) { var array = new Array(); var j = 0; for (var i = start; i <= end; i++) { array[j++] = i; } return array; } Array.prototype.map = function(func) { for (var i = 0; i < this.length; i++) { func(this[i]); } } Array.prototype.filter = function(func) { var array = new Array(); var j = 0; for (var i = 0; i < this.length; i++) { if (func(this[i])) array[j++] = this[i]; } return array; } var array = new Array(). range(1, 9). filter(function (n) { return (n % 3 == 0); }). map(function (n) { WScript.Echo(n); });
Array.prototype.range = function(start, end) { var array = new Array(); var j = 0; for (var i = start; i <= end; i++) { array[j++] = i; } return array; } Array.prototype.filter = function(func) { var array = new Array(); var j = 0; for (var i = 0; i < this.length; i++) { if (func(this[i])) array[j++] = this[i]; } return array; } Array.prototype.sum = function() { var res = 0 for (var i = 0; i < this.length; i++) { res += this[i]; } return res; } var sum = new Array(). range(1, 99). filter(function (n) { return (n % 3 == 0); }). sum(); WScript.Echo(sum);
// 初項:a, 公差:a で 上限:lim の数列の総和を返す関数 function sn(a, lim) { var n = parseInt(lim / a) // 項数:n = 上限:lim / 公差:a var l = n * a // 末項:l = 項数:n * 公差:a return (a + l) * n / 2 // 総和:sn = (初項:a + 末項:l) * 項数:n / 2 } // 3 の倍数の合計 WScript.Echo(sn(3, 999))
// 10000 までの 自然数の和 // 項目数 n = 10000 var n = 10000 WScript.Echo(n * (n + 1) / 2)
// 10000 までの 偶数の和 // 項目数 n = 5000 var n = 10000 / 2 WScript.Echo(n * (n + 1))
// 10000 までの 奇数の和 // 項目数 n = 5000 var n = 10000 / 2 WScript.Echo(Math.pow(n, 2))
// 1000 までの 自然数の2乗の和 var n = 1000 WScript.Echo(n * (n + 1) * (2 * n + 1) / 6)
// 100 までの 自然数の3乗の和 var n = 100 WScript.Echo(Math.pow(n, 2) * (Math.pow((n + 1), 2)) / 4)
// 初項 2, 公比 3, 項数 10 の等比数列の和 var n = 10 var a = 2 var r = 3 WScript.Echo((a * (Math.pow(r, n) - 1)) / (r - 1))
var a = 5 // 初項 5 var d = 3 // 公差 3 var n = 10 // 項数 10 var p = 1 // 積 for (var i = 1; i <= n; i++) { var m = a + (d * (i - 1)) p *= m } WScript.Echo(p)
// 初項 5, 公差 3, 項数 10 の数列の積を表示する WScript.Echo(product(5, 3, 10)) function product(m, d, n) { if (n == 0) return 1 else return m * product(m + d, d, n - 1) }
// 階乗を求める関数 function Fact(n) { if (n <= 1) return 1 else return n * Fact(n - 1) } // 10の階乗 WScript.Echo(Fact(10)) WScript.Echo(10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1)
// 下降階乗冪 function FallingFact(x, n) { if (n <= 1) return x else return x * FallingFact(x - 1, n - 1) } // 10 から 6 までの 総乗 WScript.Echo(FallingFact(10, 5)) WScript.Echo(10 * 9 * 8 * 7 * 6)
// 上昇階乗冪 function RisingFact(x, n) { if (n <= 1) return x else return x * RisingFact(x + 1, n - 1) } // 10 から 14 までの 総乗 WScript.Echo(RisingFact(10, 5)) WScript.Echo(10 * 11 * 12 * 13 * 14)
// 階乗 function Fact(n) { if (n <= 1) return 1 else return n * Fact(n - 1) } // 下降階乗冪 function FallingFact(x, n) { if (n <= 1) return x else return x * FallingFact(x - 1, n - 1) } // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) var n = 10 var r = 5 WScript.Echo(Fact(n) / Fact(n - r)) WScript.Echo(FallingFact(n, r))
// 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) var n = 10 var r = 5 WScript.Echo(Math.pow(n, r))
// 組合せ function Comb(n, r) { if (r == 0 || r == n) return 1 else if (r == 1) return n else return Comb(n - 1, r - 1) + Comb(n - 1, r) } // 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) var n = 10 var r = 5 WScript.Echo(Comb(n, r))
// 組合せ function Comb(n, r) { if (r == 0 || r == n) return 1 else if (r == 1) return n else return Comb(n - 1, r - 1) + Comb(n - 1, r) } // 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) var n = 10 var r = 5 WScript.Echo(Comb(n + r - 1, r))
for (var degree = 0; degree <= 360; degree += 15) { if (degree % 30 == 0 || degree % 45 == 0) { var radian = degree * Math.PI / 180.0; // 自作の正弦関数 var d1 = mySin(radian, 1, false, radian, 1.0, radian); // 標準の正弦関数 var d2 = Math.sin(radian); // 標準関数との差異 WScript.StdOut.Write((" " + degree ).slice( -3) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } } // 自作の正弦関数 function mySin(x, n, nega, numerator, denominator, y) { var m = 2 * n; denominator = denominator * (m + 1) * m; numerator = numerator * x * x; var a = numerator / denominator; // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) return y; else return y + mySin(x, ++n, !nega, numerator, denominator, nega ? a : -a); }
for (var degree = 0; degree <= 360; degree += 15) { if (degree % 30 == 0 || degree % 45 == 0) { var radian = degree * Math.PI / 180.0; // 自作の余弦関数 var d1 = myCos(radian, 1, false, 1.0, 1.0, 1.0); // 標準の余弦関数 var d2 = Math.cos(radian); // 標準関数との差異 WScript.StdOut.Write((" " + degree ).slice( -3) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } } // 自作の余弦関数 function myCos(x, n, nega, numerator, denominator, y) { var m = 2 * n; denominator = denominator * m * (m - 1); numerator = numerator * x * x; var a = numerator / denominator; // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) return y; else return y + myCos(x, ++n, !nega, numerator, denominator, nega ? a : -a); }
for (var degree = -90; degree <= 90; degree += 15) { if ((degree + 90) % 180 != 0) { var radian = degree * Math.PI / 180.0; var x2 = radian * radian; // 自作の正接関数 var d1 = myTan(radian, x2, 15, 0.0); // 15:必要な精度が得られる十分大きな奇数 // 標準の正接関数 var d2 = Math.tan(radian); // 標準関数との差異 WScript.StdOut.Write((" " + degree ).slice( -3) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } } // 自作の正接関数 function myTan(x, x2, n, t) { t = x2 / (n - t); n -= 2; if (n <= 1) return x / (1 - t); else return myTan(x, x2, n, t); }
for (var i = -10; i <= 10; i++) { var x = i / 4.0; // 標準の指数関数 var d1 = Math.exp(x); // 自作の指数関数 var d2 = myExp(x, 1, 1.0, 1.0, 1.0); // 標準関数との差異 WScript.StdOut.Write((" " + x.toFixed(2) ).slice(-5) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } // 自作の指数関数 function myExp(x, n, numerator, denominator, y) { denominator = denominator * n; numerator = numerator * x; var a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) return y; else return y + myExp(x, ++n, numerator, denominator, a); }
for (var i = -10; i <= 10; i++) { var x = i / 4.0; // 標準の指数関数 var d1 = Math.exp(x); // 自作の指数関数 var x2 = x * x; var d2 = myExp(x, x2, 30, 0.0); // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる // 標準関数との差異 WScript.StdOut.Write((" " + x.toFixed(2) ).slice(-5) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } // 自作の指数関数 function myExp(x, x2, n, t) { t = x2 / (n + t); n -= 4; if (n < 6) return 1 + ((2 * x) / (2 - x + t)); else return myExp(x, x2, n, t); }
for (var i = 1; i <= 20; i++) { var x = i / 5.0; // 標準の対数関数 var d1 = Math.log(x); // 自作の対数関数 var x2 = (x - 1) / (x + 1); var d2 = 2 * myLog(x2, x2, 1.0, x2); // 標準関数との差異 WScript.StdOut.Write((" " + x.toFixed(2) ).slice(-5) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } // 自作の対数関数 function myLog(x2, numerator, denominator, y) { denominator = denominator + 2; numerator = numerator * x2 * x2; var a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) return y; else return y + myLog(x2, numerator, denominator, a); }
for (var i = 1; i <= 20; i++) { var x = i / 5.0; // 標準の対数関数 var d1 = Math.log(x); // 自作の対数関数 var d2 = myLog(x - 1, 27, 0.0); // 27:必要な精度が得られる十分大きな奇数 // 標準関数との差異 WScript.StdOut.Write((" " + x.toFixed(2) ).slice(-5) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-13) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-13) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } // 自作の対数関数 function myLog(x, n, t) { var n2 = n; var x2 = x; if (n > 3) { if (n % 2 == 0) n2 = 2; x2 = x * parseInt(n / 2); } t = x2 / (n2 + t); if (n <= 2) return x / (1 + t); else return myLog(x, --n, t); }
for (var x = -10; x <= 10; x++) { // 自作の双曲線正弦関数 var d1 = mySinh(x, 1, x, 1.0, x); // 標準の双曲線正弦関数はない var d2 = (Math.exp(x) - Math.exp(-x)) / 2; // 標準関数との差異 WScript.StdOut.Write((" " + x ).slice( -3) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-17) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-17) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } // 自作の双曲線正弦関数 function mySinh(x, n, numerator, denominator, y) { var m = 2 * n; denominator = denominator * (m + 1) * m; numerator = numerator * x * x; var a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) return y; else return y + mySinh(x, ++n, numerator, denominator, a); }
for (var x = -10; x <= 10; x++) { // 自作の双曲線余弦関数 var d1 = myCosh(x, 1, 1.0, 1.0, 1.0); // 標準の双曲線余弦関数はない var d2 = (Math.exp(x) + Math.exp(-x)) / 2; // 標準関数との差異 WScript.StdOut.Write((" " + x ).slice( -3) + " : "); WScript.StdOut.Write((" " + d1.toFixed(10) ).slice(-17) + " - "); WScript.StdOut.Write((" " + d2.toFixed(10) ).slice(-17) + " = "); WScript.StdOut.Write((" " + (d1 - d2).toFixed(10)).slice(-13) + "\n" ); } // 自作の双曲線余弦関数 function myCosh(x, n, numerator, denominator, y) { var m = 2 * n; denominator = denominator * m * (m - 1); numerator = numerator * x * x; var a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) return y; else return y + myCosh(x, ++n, numerator, denominator, a); }
var a = 0 var b = 1 // 台形則で積分 var n = 2; for (var j = 1; j <= 10; j++) { var h = (b - a) / n var s = 0 var x = a for (var i = 1; i <= n - 1; i++) { x += h s += f(x) } s = h * ((f(a) + f(b)) / 2 + s) n *= 2 // 結果を π と比較 WScript.StdOut.Write((" " + j ).slice( -2) + " : "); WScript.StdOut.Write((" " + s.toFixed(10) ).slice(-13) + ", "); WScript.StdOut.Write((" " + (s - Math.PI).toFixed(10)).slice(-13) + "\n" ); } function f(x) { return 4 / (1 + x * x) }
var a = 0 var b = 1 // 中点則で積分 var n = 2; for (var j = 1; j <= 10; j++) { var h = (b - a) / n var s = 0 var x = a + (h / 2) for (var i = 1; i <= n; i++) { s += f(x) x += h } s *= h n *= 2 // 結果を π と比較 WScript.StdOut.Write((" " + j ).slice( -2) + " : "); WScript.StdOut.Write((" " + s.toFixed(10) ).slice(-13) + ", "); WScript.StdOut.Write((" " + (s - Math.PI).toFixed(10)).slice(-13) + "\n" ); } function f(x) { return 4 / (1 + x * x) }
var a = 0 var b = 1 // Simpson則で積分 var n = 2 for (var j = 1; j <= 5; j++) { var h = (b - a) / n var s2 = 0 var s4 = 0 var x = a + h for (var i = 1; i <= n / 2; i++) { s4 += f(x) x += h s2 += f(x) x += h } s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 var s = (s2 + s4) * h / 3 n *= 2 // 結果を π と比較 WScript.StdOut.Write((" " + j ).slice( -2) + " : "); WScript.StdOut.Write((" " + s.toFixed(10) ).slice(-13) + ", "); WScript.StdOut.Write((" " + (s - Math.PI).toFixed(10)).slice(-13) + "\n" ); } function f(x) { return 4 / (1 + x * x) }
var a = 0 var b = 1 var t = [] // 台形則で積分 var n = 2; for (var i = 1; i <= 6; i++) { var h = (b - a) / n var s = 0 var x = a for (var j = 1; j <= n - 1; j++) { x += h s += f(x) } // 結果を保存 t[i] = [] t[i][1] = h * ((f(a) + f(b)) / 2 + s) n *= 2 } // Richardsonの補外法 n = 4 for (var j = 2; j <= 6; j++) { for (var i = j; i <= 6; i++) { t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1) if (i == j) { // 結果を π と比較 WScript.StdOut.Write((" " + j ).slice( -2) + " : "); WScript.StdOut.Write((" " + t[i][j].toFixed(10) ).slice(-13) + ", "); WScript.StdOut.Write((" " + (t[i][j] - Math.PI).toFixed(10)).slice(-13) + "\n" ); } } n *= 4 } function f(x) { return 4 / (1 + x * x) }
// データ点の数 var N = 7 var x = [] var y = [] // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (var i = 0; i < N; i++) { var d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 0.5刻みで 与えられていない値を補間 for (var i = 0; i <= 18; i++) { var d1 = i * 0.5 - 4.5 var d2 = f(d1) var d3 = lagrange(d1, x, y) // 元の関数と比較 WScript.StdOut.Write((" " + d1.toFixed(2) ).slice(-5) + "\t") WScript.StdOut.Write((" " + d2.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + d3.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + (d2 - d3).toFixed(5)).slice(-8) + "\n") } // 元の関数 function f(x) { return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Lagrange (ラグランジュ) 補間 function lagrange(d, x, y) { var sum = 0; for (var i = 0; i < N; i++) { prod = y[i] for (var j = 0; j < N; j++) { if (j != i) prod *= (d - x[j]) / (x[i] - x[j]) } sum += prod } return sum }
// データ点の数 var N = 7 var x = [] var y = [] // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (var i = 0; i < N; i++) { var d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 0.5刻みで 与えられていない値を補間 for (var i = 0; i <= 18; i++) { var d1 = i * 0.5 - 4.5 var d2 = f(d1) var d3 = neville(d1, x, y) // 元の関数と比較 WScript.StdOut.Write((" " + d1.toFixed(2) ).slice(-5) + "\t") WScript.StdOut.Write((" " + d2.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + d3.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + (d2 - d3).toFixed(5)).slice(-8) + "\n") } // 元の関数 function f(x) { return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2); } // Neville (ネヴィル) 補間 function neville(d, x, y) { w = [] w[0] = [] for (var i = 0; i < N; i++) w[0][i] = y[i] for (var j = 1; j < N; j++) { w[j] = [] for (var i = 0; i < N - j; i++) w[j][i] = w[j-1][i+1] + (w[j-1][i+1] - w[j-1][i]) * (d - x[i+j]) / (x[i+j] - x[i]) } return w[N-1][0] }
// データ点の数 var N = 7 var x = [] var y = [] // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (var i = 0; i < N; i++) { var d = i * 1.5 - 4.5 x[i] = d y[i] = f(d) } // 差分商の表を作る var d = [] d[0] = [] for (var j = 0; j < N; j++) d[0][j] = y[j] for (var i = 1; i < N; i++) { d[i] = [] for (var j = 0; j < N - i; j++) d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j]) } // n階差分商 var a = [] for (var j = 0; j < N; j++) a[j] = d[j][0] // 0.5刻みで 与えられていない値を補間 for (var i = 0; i <= 18; i++) { var d1 = i * 0.5 - 4.5 var d2 = f(d1) var d3 = newton(d1, x, a) // 元の関数と比較 WScript.StdOut.Write((" " + d1.toFixed(2) ).slice(-5) + "\t") WScript.StdOut.Write((" " + d2.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + d3.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + (d2 - d3).toFixed(5)).slice(-8) + "\n") } // 元の関数 function f(x) { return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Newton (ニュートン) 補間 function newton(d, x, a) { var sum = a[0] for (var i = 1; i < N; i++) { var prod = a[i] for (var j = 0; j < i; j++) prod *= (d - x[j]) sum += prod } return sum }
// データ点の数 var N = 7 var Nx2 = 14 var x = [] var y = [] var yd = [] // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (var i = 0; i < N; i++) { var d1 = i * 1.5 - 4.5 x[i] = d1 y[i] = f(d1) yd[i] = fd(d1) } // 差分商の表を作る var z = [] var d = [] d[0] = [] for (var i = 0; i < Nx2; i++) { var j = parseInt(i / 2) z[i] = x[j] d[0][i] = y[j] } for (var i = 1; i < Nx2; i++) { d[i] = [] for (var j = 0; j < Nx2 - i; j++) { if (i == 1 && j % 2 == 0) d[i][j] = yd[parseInt(j / 2)] else d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]) } } // n階差分商 var a = [] for (var j = 0; j < Nx2; j++) a[j] = d[j][0] // 0.5刻みで 与えられていない値を補間 for (var i = 0; i <= 18; i++) { var d1 = i * 0.5 - 4.5 var d2 = f(d1) var d3 = hermite(d1, z, a) // 元の関数と比較 WScript.StdOut.Write((" " + d1.toFixed(2) ).slice(-5) + "\t") WScript.StdOut.Write((" " + d2.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + d3.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + (d2 - d3).toFixed(5)).slice(-8) + "\n") } // 元の関数 function f(x) { return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // 導関数 function fd(x) { return 1 - Math.pow(x,2) / 2 + Math.pow(x,4) / (4 * 3 * 2) } // Hermite (エルミート) 補間 function hermite(d, z, a) { var sum = a[0] for (var i = 1; i < Nx2; i++) { var prod = a[i] for (var j = 0; j < i; j++) prod *= (d - z[j]) sum += prod } return sum }
// データ点の数 var N = 7 var x = [] var y = [] // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (var i = 0; i < N; i++) { var d1 = i * 1.5 - 4.5 x[i] = d1 y[i] = f(d1) } // 3項方程式の係数の表を作る var a = [] var b = [] var c = [] var d = [] for (var i = 1; i < N - 1; i++) { a[i] = x[i] - x[i-1] b[i] = 2.0 * (x[i+1] - x[i-1]) c[i] = x[i+1] - x[i] d[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / (x[i] - x[i-1])) } // 3項方程式を解く (ト−マス法) var g = [] var s = [] g[1] = b[1] s[1] = d[1] for (var i = 2; i < N - 1; i++) { g[i] = b[i] - a[i] * c[i-1] / g[i-1] s[i] = d[i] - a[i] * s[i-1] / g[i-1] } var z = [] z[0] = 0 z[N-1] = 0 z[N-2] = s[N-2] / g[N-2] for (var i = N - 3; i >= 1; i--) z[i] = (s[i] - c[i] * z[i+1]) / g[i] // 0.5刻みで 与えられていない値を補間 for (var i = 0; i <= 18; i++) { var d1 = i * 0.5 - 4.5 var d2 = f(d1) var d3 = spline(d1, x, y, z) // 元の関数と比較 WScript.StdOut.Write((" " + d1.toFixed(2) ).slice(-5) + "\t") WScript.StdOut.Write((" " + d2.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + d3.toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + (d2 - d3).toFixed(5)).slice(-8) + "\n") } // 元の関数 function f(x) { return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Spline (スプライン) 補間 function spline(d, x, y, z) { // 補間関数値がどの区間にあるか var k = -1 for (var i = 1; i < N; i++) { if (d <= x[i]) { k = i - 1 break } } if (k < 0) k = N - 1 var d1 = x[k+1] - d var d2 = d - x[k] var d3 = x[k+1] - x[k] return (z[k] * Math.pow(d1,3) + z[k+1] * Math.pow(d2,3)) / (6.0 * d3) + (y[k] / d3 - z[k] * d3 / 6.0) * d1 + (y[k+1] / d3 - z[k+1] * d3 / 6.0) * d2 }
// 重力加速度 var g = -9.8 // 空気抵抗係数 var k = -0.01 // 時間間隔(秒) var h = 0.01 // 角度 var degree = 45 var radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 var v = parseInt(250 * 1000 / 3600) // 水平方向の速度 var vx = [] vx[0] = v * Math.cos(radian) // 鉛直方向の速度 var vy = [] vy[0] = v * Math.sin(radian) // 経過秒数 var t = 0.0 // 位置 var x = 0.0 var y = 0.0 // Euler法 for (var i = 1; y >= 0.0; i++) { // 経過秒数 t = i * h // 位置 x += h * vx[0] y += h * vy[0] WScript.StdOut.Write((" " + t.toFixed(2) ).slice(-4) + "\t") WScript.StdOut.Write((" " + vx[0].toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + vy[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + x.toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + y.toFixed(5) ).slice(-8) + "\n") // 速度 vx[1] = vx[0] + h * fx(vx[0], vy[0]) vy[1] = vy[0] + h * fy(vx[0], vy[0]) vx[0] = vx[1] vy[0] = vy[1] } // 空気抵抗による水平方向の減速分 function fx(vx, vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 function fy(vx, vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy) }
// 重力加速度 var g = -9.8 // 空気抵抗係数 var k = -0.01 // 時間間隔(秒) var h = 0.01 // 角度 var degree = 45 var radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 var v = parseInt(250 * 1000 / 3600) // 水平方向の速度 var vx = [] vx[0] = v * Math.cos(radian) // 鉛直方向の速度 var vy = [] vy[0] = v * Math.sin(radian) // 経過秒数 var t = 0.0 // 位置 var x = [] x[0] = 0.0 var y = [] y[0] = 0.0 // Heun法 for (var i = 1; y[0] >= 0.0; i++) { // 経過秒数 t = i * h // 位置・速度 x[1] = x[0] + h * vx[0] y[1] = y[0] + h * vy[0] vx[1] = vx[0] + h * fx(vx[0], vy[0]) vy[1] = vy[0] + h * fy(vx[0], vy[0]) x[2] = x[0] + h * ( vx[0] + vx[1] ) / 2 y[2] = y[0] + h * ( vy[0] + vy[1] ) / 2 vx[2] = vx[0] + h * (fx(vx[0], vy[0]) + fx(vx[1], vy[1])) / 2 vy[2] = vy[0] + h * (fy(vx[0], vy[0]) + fy(vx[1], vy[1])) / 2 x[0] = x[2] y[0] = y[2] vx[0] = vx[2] vy[0] = vy[2] WScript.StdOut.Write((" " + t.toFixed(2) ).slice(-4) + "\t") WScript.StdOut.Write((" " + vx[0].toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + vy[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + x[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + y[0].toFixed(5) ).slice(-8) + "\n") } // 空気抵抗による水平方向の減速分 function fx(vx, vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 function fy(vx, vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy) }
// 重力加速度 var g = -9.8 // 空気抵抗係数 var k = -0.01 // 時間間隔(秒) var h = 0.01 // 角度 var degree = 45 var radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 var v = parseInt(250 * 1000 / 3600) // 水平方向の速度 var vx = [] vx[0] = v * Math.cos(radian) // 鉛直方向の速度 var vy = [] vy[0] = v * Math.sin(radian) // 経過秒数 var t = 0.0 // 位置 var x = [] x[0] = 0.0 var y = [] y[0] = 0.0 // 中点法 for (var i = 1; y[0] >= 0.0; i++) { // 経過秒数 t = i * h // 位置・速度 vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) var wx = vx[0] + vx[1] / 2.0 var wy = vy[0] + vy[1] / 2.0 vx[0] = vx[0] + h * fx(wx, wy) vy[0] = vy[0] + h * fy(wx, wy) x[0] = x[0] + h * wx y[0] = y[0] + h * wy WScript.StdOut.Write((" " + t.toFixed(2) ).slice(-4) + "\t") WScript.StdOut.Write((" " + vx[0].toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + vy[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + x[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + y[0].toFixed(5) ).slice(-8) + "\n") } // 空気抵抗による水平方向の減速分 function fx(vx, vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 function fy(vx, vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy) }
// 重力加速度 var g = -9.8 // 空気抵抗係数 var k = -0.01 // 時間間隔(秒) var h = 0.01 // 角度 var degree = 45 var radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 var v = parseInt(250 * 1000 / 3600) // 水平方向の速度 var vx = [] vx[0] = v * Math.cos(radian) // 鉛直方向の速度 var vy = [] vy[0] = v * Math.sin(radian) // 経過秒数 var t = 0.0 // 位置 var x = [] x[0] = 0.0 var y = [] y[0] = 0.0 // Runge-Kutta法 for (var i = 1; y[0] >= 0.0; i++) { // 経過秒数 t = i * h // 位置・速度 x[1] = h * vx[0] y[1] = h * vy[0] vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) var wx = vx[0] + vx[1] / 2 var wy = vy[0] + vy[1] / 2 x[2] = h * wx y[2] = h * wy vx[2] = h * fx(wx, wy) vy[2] = h * fy(wx, wy) wx = vx[0] + vx[2] / 2 wy = vy[0] + vy[2] / 2 x[3] = h * wx y[3] = h * wy vx[3] = h * fx(wx, wy) vy[3] = h * fy(wx, wy) wx = vx[0] + vx[3] wy = vy[0] + vy[3] x[4] = h * wx y[4] = h * wy vx[4] = h * fx(wx, wy) vy[4] = h * fy(wx, wy) x[0] += ( x[1] + x[2] * 2 + x[3] * 2 + x[4]) / 6 y[0] += ( y[1] + y[2] * 2 + y[3] * 2 + y[4]) / 6 vx[0] += (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6 vy[0] += (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6 WScript.StdOut.Write((" " + t.toFixed(2) ).slice(-4) + "\t") WScript.StdOut.Write((" " + vx[0].toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + vy[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + x[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + y[0].toFixed(5) ).slice(-8) + "\n") } // 空気抵抗による水平方向の減速分 function fx(vx, vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 function fy(vx, vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy) }
// 重力加速度 var g = -9.8 // 空気抵抗係数 var k = -0.01 // 時間間隔(秒) var h = 0.01 // 角度 var degree = 45 var radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 var v = parseInt(250 * 1000 / 3600) // 水平方向の速度 var vx = [] vx[0] = v * Math.cos(radian) // 鉛直方向の速度 var vy = [] vy[0] = v * Math.sin(radian) // 経過秒数 var t = 0.0 // 位置 var x = [] x[0] = 0.0 var y = [] y[0] = 0.0 // Runge-Kutta-Gill法 for (var i = 1; y[0] >= 0.0; i++) { // 経過秒数 t = i * h // 位置・速度 x[1] = h * vx[0] y[1] = h * vy[0] vx[1] = h * fx(vx[0], vy[0]) vy[1] = h * fy(vx[0], vy[0]) var wx = vx[0] + vx[1] / 2 var wy = vy[0] + vy[1] / 2 x[2] = h * wx y[2] = h * wy vx[2] = h * fx(wx, wy) vy[2] = h * fy(wx, wy) wx = vx[0] + vx[1] * ((Math.sqrt(2.0) - 1) / 2) + vx[2] * (1 - 1 / Math.sqrt(2.0)) wy = vy[0] + vy[1] * ((Math.sqrt(2.0) - 1) / 2) + vy[2] * (1 - 1 / Math.sqrt(2.0)) x[3] = h * wx y[3] = h * wy vx[3] = h * fx(wx, wy) vy[3] = h * fy(wx, wy) wx = vx[0] - vx[2] / Math.sqrt(2.0) + vx[3] * (1 + 1 / Math.sqrt(2.0)) wy = vy[0] - vy[2] / Math.sqrt(2.0) + vy[3] * (1 + 1 / Math.sqrt(2.0)) x[4] = h * wx y[4] = h * wy vx[4] = h * fx(wx, wy) vy[4] = h * fy(wx, wy) x[0] += ( x[1] + x[2] * (2 - Math.sqrt(2.0)) + x[3] * (2 + Math.sqrt(2.0)) + x[4]) / 6 y[0] += ( y[1] + y[2] * (2 - Math.sqrt(2.0)) + y[3] * (2 + Math.sqrt(2.0)) + y[4]) / 6 vx[0] += (vx[1] + vx[2] * (2 - Math.sqrt(2.0)) + vx[3] * (2 + Math.sqrt(2.0)) + vx[4]) / 6 vy[0] += (vy[1] + vy[2] * (2 - Math.sqrt(2.0)) + vy[3] * (2 + Math.sqrt(2.0)) + vy[4]) / 6 WScript.StdOut.Write((" " + t.toFixed(2) ).slice(-4) + "\t") WScript.StdOut.Write((" " + vx[0].toFixed(5) ).slice(-8) + "\t") WScript.StdOut.Write((" " + vy[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + x[0].toFixed(5) ).slice(-9) + "\t") WScript.StdOut.Write((" " + y[0].toFixed(5) ).slice(-8) + "\n") } // 空気抵抗による水平方向の減速分 function fx(vx, vy) { return k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 function fy(vx, vy) { return g + (k * Math.sqrt(vx * vx + vy * vy) * vy) }
var a = 1 var b = 2 WScript.StdOut.Write((" " + bisection(a, b).toFixed(10) ).slice(-12) + "\n") function bisection(a, b) { var c while (true) { // 区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2 WScript.StdOut.Write((" " + c.toFixed(10) ).slice(-12) + "\t") WScript.StdOut.Write((" " + (c - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n") var fc = f(c) if (Math.abs(fc) < 0.0000000001) break if (fc < 0){ // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c } } return c } function f(x) { return x * x - 2 }
var a = 1 var b = 2 WScript.StdOut.Write((" " + falseposition(a, b).toFixed(10) ).slice(-12) + "\n") function falseposition(a, b) { var c while (true) { // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c = (a * f(b) - b * f(a)) / (f(b) - f(a)) WScript.StdOut.Write((" " + c.toFixed(10) ).slice(-12) + "\t") WScript.StdOut.Write((" " + (c - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n") var fc = f(c) if (Math.abs(fc) < 0.0000000001) break if (fc < 0){ // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c } } return c } function f(x) { return x * x - 2 }
var x = 1 WScript.StdOut.Write((" " + iterative(x).toFixed(10) ).slice(-12) + "\n") function iterative(x0) { var x1 while (true) { x1 = g(x0) WScript.StdOut.Write((" " + x1.toFixed(10) ).slice(-12) + "\t") WScript.StdOut.Write((" " + (x1 - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n") if (Math.abs(x1 - x0) < 0.0000000001) break x0 = x1 } return x1 } function g(x) { return (x / 2) + (1 / x) }
var x = 2 WScript.StdOut.Write((" " + newton(x).toFixed(10) ).slice(-12) + "\n") function newton(x0) { var x1 while (true) { x1 = x0 - (f0(x0) / f1(x0)) WScript.StdOut.Write((" " + x1.toFixed(10) ).slice(-12) + "\t") WScript.StdOut.Write((" " + (x1 - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n") if (Math.abs(x1 - x0) < 0.0000000001) break x0 = x1 } return x1 } function f0(x) { return x * x - 2 } function f1(x) { return 2 * x }
var x = 2 WScript.StdOut.Write((" " + bailey(x).toFixed(10) ).slice(-12) + "\n") function bailey(x0) { var x1 while (true) { x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))) WScript.StdOut.Write((" " + x1.toFixed(10) ).slice(-12) + "\t") WScript.StdOut.Write((" " + (x1 - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n") if (Math.abs(x1 - x0) < 0.0000000001) break x0 = x1 } return x1 } function f0(x) { return x * x - 2 } function f1(x) { return 2 * x } function f2(x) { return 2 }
var x0 = 1 var x1 = 2 WScript.StdOut.Write((" " + secant(x0, x1).toFixed(10) ).slice(-12) + "\n") function secant(x0, x1) { var x2 while (true) { x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) WScript.StdOut.Write((" " + x2.toFixed(10) ).slice(-12) + "\t") WScript.StdOut.Write((" " + (x2 - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n") if (Math.abs(x2 - x1) < 0.0000000001) break x0 = x1 x1 = x2 } return x2 } function f(x) { return x * x - 2 }
var N = 4 var a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]] var b = [20,16,8,17] var c = [0,0,0,0] // ヤコビの反復法 jacobi(a,b,c) WScript.StdOut.WriteLine("X") disp_vector(c) // ヤコビの反復法 function jacobi(a, b, x0) { while (true) { x1 = [] var finish = true for (i = 0; i < N; i++) { x1[i] = 0 for (j = 0; j < N; j++) if (j != i) x1[i] += a[i][j] * x0[j] x1[i] = (b[i] - x1[i]) / a[i][i] if (Math.abs(x1[i] - x0[i]) > 0.0000000001) finish = false } for (i = 0; i < N; i++) x0[i] = x1[i] if (finish) return disp_vector(x0) } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() }
var N = 4 var a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]] var b = [20,16,8,17] var c = [0,0,0,0] // ガウス・ザイデル法 gauss(a,b,c) WScript.StdOut.WriteLine("X") disp_vector(c) // ガウス・ザイデル法 function gauss(a, b, x0) { while (true) { var finish = true for (i = 0; i < N; i++) { var x1 = 0 for (j = 0; j < N; j++) if (j != i) x1 += a[i][j] * x0[j] x1 = (b[i] - x1) / a[i][i] if (Math.abs(x1 - x0[i]) > 0.0000000001) finish = false x0[i] = x1 } if (finish) return disp_vector(x0) } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() }
var N = 4 var a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]] var b = [8,17,20,16] // ピボット選択 pivoting(a,b) WScript.StdOut.WriteLine("pivoting") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 前進消去 forward_elimination(a,b) WScript.StdOut.WriteLine("forward elimination") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 後退代入 backward_substitution(a,b) WScript.StdOut.WriteLine("X") disp_vector(b) // 前進消去 function forward_elimination(a, b) { for (pivot = 0; pivot < N - 1; pivot++) { for (row = pivot + 1; row < N; row++) { var s = a[row][pivot] / a[pivot][pivot] for (col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s } } } // 後退代入 function backward_substitution(a, b) { for (row = N - 1; row >= 0; row--) { for (col = N - 1; col > row; col--) b[row] -= a[row][col] * b[col] b[row] /= a[row][row] } } // ピボット選択 function pivoting(a, b) { for (pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val = 0 for (row = pivot; row < N; row++) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp for (col = 0; col < N; col++) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + matrix[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
var N = 4 var a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]] var b = [8,17,20,16] // ピボット選択 pivoting(a,b) WScript.StdOut.WriteLine("pivoting") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 前進消去 forward_elimination(a,b) WScript.StdOut.WriteLine("forward elimination") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 後退代入 backward_substitution(a,b) WScript.StdOut.WriteLine("X") disp_vector(b) // 前進消去 function forward_elimination(a, b) { for (pivot = 0; pivot < N; pivot++) { for (row = 0; row < N; row++) { if (row == pivot) continue var s = a[row][pivot] / a[pivot][pivot] for (col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s b[row] -= b[pivot] * s } } } // 後退代入 function backward_substitution(a, b) { for (pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot][pivot] } // ピボット選択 function pivoting(a, b) { for (pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val = 0 for (row = pivot; row < N; row++) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp for (col = 0; col < N; col++) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 状態を確認 function disp_progress(a, b) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + a[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine((" " + b[i].toFixed(10)).slice(-14)) } WScript.StdOut.WriteLine() } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + matrix[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
var N = 4 var a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]] var b = [8,17,20,16] // ピボット選択 pivoting(a,b) WScript.StdOut.WriteLine("pivoting") WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // LU分解 // 前進消去 forward_elimination(a,b) WScript.StdOut.WriteLine("LU") disp_matrix(a) // Ly=b から y を求める (前進代入) var y = [0,0,0,0] forward_substitution(a,b,y) WScript.StdOut.WriteLine("Y") disp_vector(y) // Ux=y から x を求める (後退代入) var x = [0,0,0,0] backward_substitution(a,y,x) WScript.StdOut.WriteLine("X") disp_vector(x) // 前進消去 function forward_elimination(a, b) { for (pivot = 0; pivot < N - 1; pivot++) { for (row = pivot + 1; row < N; row++) { var s = a[row][pivot] / a[pivot][pivot] for (col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s // これが 上三角行列 a[row][pivot] = s // これが 下三角行列 // b[row] -= b[pivot] * s // この値は変更しない } } } // 前進代入 function forward_substitution(a, b, y) { for (row = 0; row < N; row++) { for (col = 0; col < row; col++) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // 後退代入 function backward_substitution(a, y, x) { for (row = N - 1; row >= 0; row--) { for (col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // ピボット選択 function pivoting(a, b) { for (pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す var max_row = pivot var max_val = 0 for (row = pivot; row < N; row++) { if (Math.abs(a[row][pivot]) > max_val) { // 一番値が大きい行 max_val = Math.abs(a[row][pivot]) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp for (col = 0; col < N; col++) { tmp = a[max_row][col] a[max_row][col] = a[pivot][col] a[pivot][col] = tmp } tmp = b[max_row] b[max_row] = b[pivot] b[pivot] = tmp } } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + matrix[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
var N = 4 var a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]] var b = [34,68,96,125] WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 前進消去 forward_elimination(a,b) WScript.StdOut.WriteLine("LL^T") disp_matrix(a) // Ly=b から y を求める (前進代入) var y = [0,0,0,0] forward_substitution(a,b,y) WScript.StdOut.WriteLine("Y") disp_vector(y) // L^Tx=y から x を求める (後退代入) var x = [0,0,0,0] backward_substitution(a,y,x) WScript.StdOut.WriteLine("X") disp_vector(x) // 前進消去 function forward_elimination(a, b) { for (pivot = 0; pivot < N; pivot++) { var s = 0 for (col = 0; col < pivot; col++) s += a[pivot][col] * a[pivot][col] // ここで根号の中が負の値になると計算できない! a[pivot][pivot] = Math.sqrt(a[pivot][pivot] - s) for (row = pivot + 1; row < N; row++) { s = 0 for (col = 0; col < pivot; col++) s += a[row][col] * a[pivot][col] a[row][pivot] = (a[row][pivot] - s) / a[pivot][pivot] a[pivot][row] = a[row][pivot] } } } // 前進代入 function forward_substitution(a, b, y) { for (row = 0; row < N; row++) { for (col = 0; col < row; col++) b[row] -= a[row][col] * y[col] y[row] = b[row] / a[row][row] } } // 後退代入 function backward_substitution(a, y, x) { for (row = N - 1; row >= 0; row--) { for (col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + matrix[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
var N = 4 var a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]] var b = [34,68,96,125] WScript.StdOut.WriteLine("A") disp_matrix(a) WScript.StdOut.WriteLine("B") disp_vector(b) WScript.StdOut.WriteLine() // 前進消去 forward_elimination(a,b) WScript.StdOut.WriteLine("LDL^T") disp_matrix(a) // Ly=b から y を求める (前進代入) var y = [0,0,0,0] forward_substitution(a,b,y) WScript.StdOut.WriteLine("Y") disp_vector(y) // DL^Tx=y から x を求める (後退代入) var x = [0,0,0,0] backward_substitution(a,y,x) WScript.StdOut.WriteLine("X") disp_vector(x) // 前進消去 function forward_elimination(a, b) { for (pivot = 0; pivot < N; pivot++) { var s // pivot < k の場合 for (col = 0; col < pivot; col++) { s = a[pivot][col] for (k = 0; k < col; k++) s -= a[pivot][k] * a[col][k] * a[k][k] a[pivot][col] = s / a[col][col] a[col][pivot] = a[pivot][col] } // pivot == k の場合 s = a[pivot][pivot] for (k = 0; k < pivot; k++) s -= a[pivot][k] * a[pivot][k] * a[k][k] a[pivot][pivot] = s } } // 前進代入 function forward_substitution(a, b, y) { for (row = 0; row < N; row++) { for (col = 0; col < row; col++) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // 後退代入 function backward_substitution(a, y, x) { for (row = N - 1; row >= 0; row--) { for (col = N - 1; col > row; col--) y[row] -= a[row][col] * a[row][row] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (i = 0; i < N; i++) { for (j = 0; j < N; j++) WScript.StdOut.Write((" " + matrix[i][j].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }
var N = 4 // ベキ乗法で最大固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var x = [1.0 ,0.0 ,0.0 ,0.0] // ベキ乗法 var lambda = power(a, x) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") WScript.StdOut.WriteLine((" " + lambda.toFixed(10)).slice(-14)) WScript.StdOut.WriteLine("eigenvector") disp_vector(x) } // ベキ乗法 function power(a, x0) { var lambda = 0.0 // 正規化 (ベクトル x0 の長さを1にする) normarize(x0) var e0 = 0.0 for (var i = 0; i < N; i++) e0 += x0[i] for (var k = 1; k <= 200; k++) { // 1次元配列を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_vector(x0) // 行列の積 x1 = A × x0 var x1 = [0.0 ,0.0 ,0.0 ,0.0] for (var i = 0; i < N; i++) for (var j = 0; j < N; j++) x1[i] += a[i][j] * x0[j] // 内積 var p0 = 0.0 var p1 = 0.0 for (var i = 0; i < N; i++) { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p0 / p1 // 正規化 (ベクトル x1 の長さを1にする) normarize(x1) // 収束判定 var e1 = 0.0 for (var i = 0; i < N; i++) e1 += x1[i] if (Math.abs(e0 - e1) < 0.00000000001) break for (var i = 0; i < N; i++) x0[i] = x1[i] e0 = e1 } return lambda } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 正規化 (ベクトルの長さを1にする) function normarize(x) { var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] s = Math.sqrt(s) for (var i = 0; i < N; i++) x[i] /= s }
var N = 4 // 逆ベキ乗法で最小固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var x = [1.0 ,0.0 ,0.0 ,0.0] // LU分解 forward_elimination(a) // 逆ベキ乗法 var lambda = inverse(a, x) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") WScript.StdOut.WriteLine((" " + lambda.toFixed(10)).slice(-14)) WScript.StdOut.WriteLine("eigenvector") disp_vector(x) } // 逆ベキ乗法 function inverse(a, x0) { var lambda = 0.0 // 正規化 (ベクトル x0 の長さを1にする) normarize(x0) var e0 = 0.0 for (var i = 0; i < N; i++) e0 += x0[i] for (var k = 1; k <= 200; k++) { // 1次元配列を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_vector(x0) // Ly=b から y を求める (前進代入) var b = [0,0,0,0] var y = [0,0,0,0] for (var i = 0; i < N; i++) b[i] = x0[i] forward_substitution(a,y,b) // Ux=y から x を求める (後退代入) var x1 = [0,0,0,0] backward_substitution(a,x1,y) // 内積 var p0 = 0.0 var p1 = 0.0 for (var i = 0; i < N; i++) { p0 += x1[i] * x1[i] p1 += x1[i] * x0[i] } // 固有値 lambda = p1 / p0 // 正規化 (ベクトル x1 の長さを1にする) normarize(x1) // 収束判定 var e1 = 0.0 for (var i = 0; i < N; i++) e1 += x1[i] if (Math.abs(e0 - e1) < 0.00000000001) break for (var i = 0; i < N; i++) x0[i] = x1[i] e0 = e1 } return lambda } // 前進消去 function forward_elimination(a) { for (var pivot = 0; pivot < N - 1; pivot++) { for (var row = pivot + 1; row < N; row++) { var s = a[row][pivot] / a[pivot][pivot] for (var col = pivot; col < N; col++) a[row][col] -= a[pivot][col] * s // これが 上三角行列 a[row][pivot] = s // これが 下三角行列 } } } // 前進代入 function forward_substitution(a, y, b) { for (var row = 0; row < N; row++) { for (var col = 0; col < row; col++) b[row] -= a[row][col] * y[col] y[row] = b[row] } } // 後退代入 function backward_substitution(a, x, y) { for (var row = N - 1; row >= 0; row--) { for (var col = N - 1; col > row; col--) y[row] -= a[row][col] * x[col] x[row] = y[row] / a[row][row] } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 正規化 (ベクトルの長さを1にする) function normarize(x) { var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] s = Math.sqrt(s) for (var i = 0; i < N; i++) x[i] /= s }
var N = 4 // LR分解で固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var l = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] var u = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] for (var k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u) // 行列の積 multiply(u, l, a) // 対角要素を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_eigenvalue(a) // 収束判定 var e = 0.0 for (var i = 1; i < N; i++) for (var j = 0; j < i; j++) e += Math.abs(a[i][j]) if (e < 0.00000000001) break } WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_eigenvalue(a) } // LU分解 function decomp(a, l, u) { for (var i = 0; i < N; i++) { for (var j = 0; j < N; j++) { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for (var j = 0; j < N; j++) u[0][j] = a[0][j] for (var i = 1; i < N; i++) { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for (var i = 1; i < N; i++) { l[i][i] = 1.0 var t = a[i][i] for (var k = 0; k <= i; k++) t -= l[i][k] * u[k][i] u[i][i] = t for (var j = i + 1; j < N; j++) { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for (var k = 0; k <= i; k++) t -= l[j][k] * u[k][i] l[j][i] = t / u[i][i] t = a[i][j] for (var k = 0; k <= i; k++) t -= l[i][k] * u[k][j] u[i][j] = t } } } // 行列の積 function multiply(a, b, c) { for (var i = 0; i < N; i++) { for (var j = 0; j < N; j++) { var s = 0.0 for (var k = 0; k < N; k++) s += a[i][k] * b[k][j] c[i][j] = s } } } // 対角要素を表示 function disp_eigenvalue(a) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + a[i][i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() }
var N = 4 // QR分解で固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var q = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] var r = [[0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0], [0.0 ,0.0 ,0.0 ,0.0]] for (var k = 1; k <= 200; k++) { // QR分解 decomp(a, q, r) // 行列の積 multiply(r, q, a) // 対角要素を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_eigenvalue(a) // 収束判定 var e = 0.0 for (var i = 1; i < N; i++) for (var j = 0; j < i; j++) e += Math.abs(a[i][j]) if (e < 0.00000000001) break } WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_eigenvalue(a) } // QR分解 function decomp(a, q, r) { var x = [0.0 ,0.0 ,0.0 ,0.0] for (var k = 0; k < N; k++) { for (var i = 0; i < N; i++) x[i] = a[i][k] for (var j = 0; j < k; j++) { var t = 0.0 for (var i = 0; i < N; i++) t += a[i][k] * q[i][j] r[j][k] = t r[k][j] = 0.0 for (var i = 0; i < N; i++) x[i] -= t * q[i][j] } var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] r[k][k] = Math.sqrt(s) for (var i = 0; i < N; i++) q[i][k] = x[i] / r[k][k] } } // 行列の積 function multiply(a, b, c) { for (var i = 0; i < N; i++) { for (var j = 0; j < N; j++) { var s = 0.0 for (var k = 0; k < N; k++) s += a[i][k] * b[k][j] c[i][j] = s } } } // 対角要素を表示 function disp_eigenvalue(a) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + a[i][i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() }
var N = 4 // ヤコビ法で固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var v = [[1.0 ,0.0 ,0.0 ,0.0], [0.0 ,1.0 ,0.0 ,0.0], [0.0 ,0.0 ,1.0 ,0.0], [0.0 ,0.0 ,0.0 ,1.0]] // ヤコビ法 jacobi(a, v) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_eigenvalue(a) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvector") disp_eigenvector(v) } // ヤコビ法 function jacobi(a, v) { for (var k = 1; k <= 100; k++) { // 最大値を探す var p = 0 var q = 0 var max_val = 0.0 for (var i = 0; i < N; i++) { for (var j = i + 1; j < N; j++) { if (max_val < Math.abs(a[i][j])) { max_val = Math.abs(a[i][j]) p = i q = j } } } // θ を求める var t = 0.0 if (Math.abs(a[p][p] - a[q][q]) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 if (a[p][p] < 0) t = -t } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0 } // θ を使って 行列 U を作成し、A = U^t × A × U var c = Math.cos(t) var s = Math.sin(t) // U^t × A var t1 = 0.0 var t2 = 0.0 for (var i = 0; i < N; i++) { t1 = a[p][i] * c + a[q][i] * s t2 = -a[p][i] * s + a[q][i] * c a[p][i] = t1 a[q][i] = t2 // 固有ベクトル t1 = v[p][i] * c + v[q][i] * s t2 = -v[p][i] * s + v[q][i] * c v[p][i] = t1 v[q][i] = t2 } // A × U for (var i = 0; i < N; i++) { t1 = a[i][p] * c + a[i][q] * s t2 = -a[i][p] * s + a[i][q] * c a[i][p] = t1 a[i][q] = t2 } // 対角要素を表示 WScript.StdOut.Write((" " + k).slice( -3) + "\t") disp_eigenvalue(a) // 収束判定 if (max_val < 0.00000000001) break } } // 対角要素を表示 function disp_eigenvalue(a) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + a[i][i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 固有ベクトルを表示 function disp_eigenvector(matrix) { for (var i = 0; i < N; i++) { var row = [0.0 ,0.0 ,0.0 ,0.0] for (var j = 0; j < N; j++) row[j] = matrix[i][j] normarize(row) disp_vector(row) } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 正規化 (ベクトルの長さを1にする) function normarize(x) { var s = 0.0 for (var i = 0; i < N; i++) s += x[i] * x[i] s = Math.sqrt(s) for (var i = 0; i < N; i++) x[i] /= s }
var N = 4 // ハウスホルダー変換とQR分解で固有値を求める main() function main() { var a = [[5.0, 4.0, 1.0, 1.0], [4.0, 5.0, 1.0, 1.0], [1.0, 1.0, 4.0, 2.0], [1.0, 1.0, 2.0, 4.0]] var d = [1.0 ,0.0 ,0.0 ,0.0] var e = [1.0 ,0.0 ,0.0 ,0.0] // ハウスホルダー変換 tridiagonalize(a, d, e) // QR分解 decomp(a, d, e) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvalue") disp_vector(d) WScript.StdOut.WriteLine() WScript.StdOut.WriteLine("eigenvector") disp_matrix(a) } // ハウスホルダー変換 function tridiagonalize(a, d, e) { var v = [0.0 ,0.0 ,0.0 ,0.0] for (var k = 0; k < N - 2; k++) { var w = [0.0 ,0.0 ,0.0 ,0.0] d[k] = a[k][k] var t = 0.0 for (var i = k + 1; i < N; i++) { w[i] = a[i][k] t += w[i] * w[i] } var s = Math.sqrt(t) if (w[k + 1] < 0) s = -s if (Math.abs(s) < 0.00000000001) e[k + 1] = 0.0 else { e[k + 1] = -s w[k + 1] += s s = 1 / Math.sqrt(w[k + 1] * s) for (var i = k + 1; i < N; i++) w[i] *= s for (var i = k + 1; i < N; i++) { s = 0.0 for (var j = k + 1; j < N; j++) { if (j <= i) s += a[i][j] * w[j] else s += a[j][i] * w[j] } v[i] = s } s = 0.0 for (var i = k + 1; i < N; i++) s += w[i] * v[i] s /= 2.0 for (var i = k + 1; i < N; i++) v[i] -= s * w[i] for (var i = k + 1; i < N; i++) for (var j = k + 1; j <= i; j++) a[i][j] -= w[i] * v[j] + w[j] * v[i] // 次の行は固有ベクトルを求めないなら不要 for (var i = k + 1; i < N; i++) a[i][k] = w[i] } } d[N - 2] = a[N - 2][N - 2] d[N - 1] = a[N - 1][N - 1] e[0] = 0.0 e[N - 1] = a[N - 1][N - 2] // 次の行は固有ベクトルを求めないなら不要 for (var k = N - 1; k >= 0; k--) { var w = [0.0 ,0.0 ,0.0 ,0.0] if (k < N - 2) { for (var i = k + 1; i < N; i++) w[i] = a[i][k] for (var i = k + 1; i < N; i++) { var s = 0.0 for (var j = k + 1; j < N; j++) s += a[i][j] * w[j] v[i] = s } for (var i = k + 1; i < N; i++) for (var j = k + 1; j < N; j++) a[i][j] -= v[i] * w[j] } for (var i = 0; i < N; i++) a[i][k] = 0.0 a[k][k] = 1.0 } } // QR分解 function decomp(a, d, e) { e[0] = 1.0 var h = N - 1 while (Math.abs(e[h]) < 0.00000000001) h-- while (h > 0) { e[0] = 0.0 var l = h - 1 while (Math.abs(e[l]) >= 0.00000000001) l-- for (var j = 1; j <= 100; j++) { var w = (d[h - 1] - d[h]) / 2.0 var s = Math.sqrt(w * w + e[h] * e[h]) if (w < 0.0) s = -s var x = d[l] - d[h] + e[h] * e[h] / (w + s) var y = e[l + 1] var z = 0.0 var t = 0.0 var u = 0.0 for (var k = l; k < h; k++) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / Math.sqrt(t * t + 1.0) if (t < 0) s = -s u = t * s } w = d[k] - d[k + 1] t = (w * s + 2 * u * e[k + 1]) * s d[k ] = d[k ] - t d[k + 1] = d[k + 1] + t e[k ] = u * e[k] - s * z e[k + 1] = e[k + 1] * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for (var i = 0; i < N; i++) { x = a[k ][i] y = a[k + 1][i] a[k ][i] = u * x - s * y a[k + 1][i] = s * x + u * y } if (k < N - 2) { x = e[k + 1] y = -s * e[k + 2] z = y e[k + 2] = u * e[k + 2] } } WScript.StdOut.Write((" " + j).slice( -3) + "\t") disp_vector(d) // 収束判定 if (Math.abs(e[h]) < 0.00000000001) break } e[0] = 1.0 while (Math.abs(e[h]) < 0.00000000001) h-- } // 次の行は固有ベクトルを求めないなら不要 for (var k = 0; k < N - 1; k++) { var l = k for (var i = k + 1; i < N; i++) if (d[i] > d[l]) l = i var t = d[k] d[k] = d[l] d[l] = t for (var i = 0; i < N; i++) { t = a[k][i] a[k][i] = a[l][i] a[l][i] = t } } } // 1次元配列を表示 function disp_vector(row) { for (var i = 0; i < N; i++) WScript.StdOut.Write((" " + row[i].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } // 2次元配列を表示 function disp_matrix(matrix) { for (var row = 0; row < N; row++) { for (var col = 0; col < N; col++) WScript.StdOut.Write((" " + matrix[row][col].toFixed(10)).slice(-14) + "\t") WScript.StdOut.WriteLine() } }