public class CS0101 { public static void Main() { System.Console.WriteLine(3 + 5); System.Console.WriteLine(3 - 5); System.Console.WriteLine(3 * 5); System.Console.WriteLine(System.Math.Pow(3, 5)); System.Console.WriteLine(5 / 3); System.Console.WriteLine(5.0 / 3); System.Console.WriteLine(5 / 3.0); System.Console.WriteLine(5 % 3); System.Console.Write(3 * 5 + "\n"); System.Console.WriteLine(string.Format("{0,3:D}", 3 * 5)); System.Console.WriteLine(string.Format("{0,23:F20}", 5 / 3.0)); } }
public class CS0102 { public static void Main() { int i = 3 * 5; System.Console.WriteLine("3 * 5 = " + i); System.Console.WriteLine(string.Format("3 * 5 = {0}", i)); } }
public class CS0103 { public static void Main() { for (int i = 1; i < 10; i++) { System.Console.Write(i + ", "); } System.Console.WriteLine(); } }
public class CS0104 { public static void Main() { for (int i = 1; i < 10; i++) { if (i % 3 == 0) { System.Console.Write(i + ", "); } } System.Console.WriteLine(); } }
public class CS0105 { public static void Main() { int sum = 0; for (int i = 1; i < 100; i++) { if (i % 3 == 0) { sum += i; } } System.Console.WriteLine(sum); } }
using System.Linq; using System.Collections.Generic; public class CS0201 { public static void Main() { IEnumerable<int> squares = Enumerable.Range(1, 9); foreach (int num in squares) { System.Console.Write(num + ", "); } System.Console.WriteLine(); } }
using System.Linq; using System.Collections.Generic; public class CS0202 { public static void Main() { IEnumerable<int> squares = Enumerable.Range(1, 9).Where(num => num % 3 == 0); foreach (int num in squares) { System.Console.Write(num + ", "); } System.Console.WriteLine(); } }
using System.Linq; public class CS0203 { public static void Main() { int sum = Enumerable.Range(1, 99).Where(num => (num % 3 == 0)).Sum(); System.Console.WriteLine(sum); sum = Enumerable.Range(1, 99).Where(num => (num % 3 == 0)).Aggregate((x, n) => x + n); System.Console.WriteLine(sum); } }
public class CS0301 { public static void Main() { // 3 の倍数の合計 System.Console.WriteLine( sn(3, 999) ); } // 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 private static int sn(int a, int lim) { int n = lim / a; // 項数:n = 上限:lim / 公差:a int l = n * a; // 末項:l = 項数:n * 公差:a return (a + l) * n / 2; // 総和:sn = (初項:a + 末項:l) * 項数:n / 2 } }
public class CS0302 { public static void Main() { // 10000 までの 自然数の和 // 項目数 n = 10000 int n = 10000; System.Console.WriteLine( n * (n + 1) / 2 ); } }
public class CS0303 { public static void Main() { // 10000 までの 偶数の和 // 項目数 n = 5000 int n = 10000 / 2; System.Console.WriteLine( n * (n + 1) ); } }
public class CS0304 { public static void Main() { // 10000 までの 奇数の和 // 項目数 n = 5000 int n = 10000 / 2; System.Console.WriteLine( System.Math.Pow(n, 2) ); } }
public class CS0305 { public static void Main() { // 1000 までの 自然数の2乗の和 int n = 1000; System.Console.WriteLine( n * (n + 1) * (2 * n + 1) / 6 ); } }
public class CS0306 { public static void Main() { // 100 までの 自然数の3乗の和 int n = 100; System.Console.WriteLine( System.Math.Pow(n, 2) * System.Math.Pow((n + 1), 2) / 4 ); } }
public class CS0307 { public static void Main() { // 初項 2, 公比 3, 項数 10 の等比数列の和 int n = 10; int a = 2; int r = 3; System.Console.WriteLine( (a * (System.Math.Pow(r, n) - 1)) / (r - 1) ); } }
using System; public class CS0401 { public static void Main() { int a = 5; // 初項 5 int d = 3; // 公差 3 int n = 10; // 項数 10 long p = 1; // 積 for (int i = 1; i <= n; i++) { int m = a + (d * (i - 1)); p *= m; } Console.WriteLine(p); } }
using System; public class CS0402 { public static void Main() { // 初項 5, 公差 3, 項数 10 の数列の積を表示する Console.WriteLine(product(5, 3, 10)); } private static long product(int m, int d, int n) { if (n == 0) return 1; else return m * product(m + d, d, n - 1); } }
using System; public class CS0403 { // 階乗を求める関数 private static int Fact(int n) { if (n <= 1) return 1; else return n * Fact(n - 1); } public static void Main() { // 10の階乗 Console.WriteLine(Fact(10)); Console.WriteLine(10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1); } }
using System; public class CS0404 { // 下降階乗冪 private static int FallingFact(int x, int n) { if (n <= 1) return x; else return x * FallingFact(x - 1, n - 1); } public static void Main() { // 10 から 6 までの 総乗 Console.WriteLine(FallingFact(10, 5)); Console.WriteLine(10 * 9 * 8 * 7 * 6); } }
using System; public class CS0405 { // 上昇階乗冪 private static int RisingFact(int x, int n) { if (n <= 1) return x; else return x * RisingFact(x + 1, n - 1); } public static void Main() { // 10 から 14 までの 総乗 Console.WriteLine(RisingFact(10, 5)); Console.WriteLine(10 * 11 * 12 * 13 * 14); } }
using System; public class CS0406 { // 階乗 private static int Fact(int n) { if (n <= 1) return 1; else return n * Fact(n - 1); } // 下降階乗冪 private static int FallingFact(int x, int n) { if (n <= 1) return x; else return x * FallingFact(x - 1, n - 1); } public static void Main() { // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) int n = 10; int r = 5; Console.WriteLine(Fact(n) / Fact(n - r)); Console.WriteLine(FallingFact(n, r)); } }
using System; public class CS0407 { public static void Main() { // 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) int n = 10; int r = 5; Console.WriteLine(Math.Pow(n, r)); } }
using System; public class CS040101 { // 組合せ private static int Comb(int n, int 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); } public static void Main() { // 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) int n = 10; int r = 5; Console.WriteLine(Comb(n, r)); } }
using System; public class CS0409 { // 組合せ private static int Comb(int n, int 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); } public static void Main() { // 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) int n = 10; int r = 5; Console.WriteLine(Comb(n + r - 1, r)); } }
using System; public class CS0501 { public static void Main() { for (int degree = 0; degree <= 360; degree += 15) { if (degree % 30 == 0 || degree % 45 == 0) { double radian = degree * Math.PI / 180.0; // 自作の正弦関数 double d1 = mySin(radian, 1, false, radian, 1.0, radian); // 標準の正弦関数 double d2 = Math.Sin(radian); // 標準関数との差異 Console.WriteLine(string.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)); } } } // 自作の正弦関数 private static double mySin(double x, int n, bool nega, double numerator, double denominator, double y) { int m = 2 * n; denominator = denominator * (m + 1) * m; numerator = numerator * x * x; double a = numerator / denominator; // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) return y; else return y + mySin(x, ++n, !nega, numerator, denominator, nega ? a : -a); } }
using System; public class CS0502 { public static void Main() { for (int degree = 0; degree <= 360; degree += 15) { if (degree % 30 == 0 || degree % 45 == 0) { double radian = degree * Math.PI / 180.0; // 自作の余弦関数 double d1 = myCos(radian, 1, false, 1.0, 1.0, 1.0); // 標準の余弦関数 double d2 = Math.Cos(radian); // 標準関数との差異 Console.WriteLine(string.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)); } } } // 自作の余弦関数 private static double myCos(double x, int n, bool nega, double numerator, double denominator, double y) { int m = 2 * n; denominator = denominator * m * (m - 1); numerator = numerator * x * x; double a = numerator / denominator; // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) return y; else return y + myCos(x, ++n, !nega, numerator, denominator, nega ? a : -a); } }
using System; public class CS0503 { public static void Main() { for (int degree = -90; degree <= 90; degree += 15) { if ((degree + 90) % 180 != 0) { double radian = degree * Math.PI / 180.0; double x2 = radian * radian; // 自作の正接関数 double d1 = myTan(radian, x2, 15, 0.0); // 15:必要な精度が得られる十分大きな奇数 // 標準の正接関数 double d2 = Math.Tan(radian); // 標準関数との差異 Console.WriteLine(string.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)); } } } // 自作の正接関数 private static double myTan(double x, double x2, int n, double t) { t = x2 / (n - t); n -= 2; if (n <= 1) return x / (1 - t); else return myTan(x, x2, n, t); } }
using System; public class CS0504 { public static void Main() { for (int i = -10; i <= 10; i++) { double x = i / 4.0; // 標準の指数関数 double d1 = Math.Exp(x); // 自作の指数関数 double d2 = myExp(x, 1, 1.0, 1.0, 1.0); // 標準関数との差異 Console.WriteLine(string.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)); } } // 自作の指数関数 private static double myExp(double x, int n, double numerator, double denominator, double y) { denominator = denominator * n; numerator = numerator * x; double a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.Abs(a) <= 0.00000000001) return y; else return y + myExp(x, ++n, numerator, denominator, a); } }
using System; public class CS0505 { public static void Main() { for (int i = -10; i <= 10; i++) { double x = i / 4.0; // 標準の指数関数 double d1 = Math.Exp(x); // 自作の指数関数 double x2 = x * x; double d2 = myExp(x, x2, 30, 0.0); // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる // 標準関数との差異 Console.WriteLine(string.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)); } } // 自作の指数関数 private static double myExp(double x, double x2, int n, double t) { t = x2 / (n + t); n -= 4; if (n < 6) return 1 + ((2 * x) / (2 - x + t)); else return myExp(x, x2, n, t); } }
using System; public class CS0506 { public static void Main() { for (int i = 1; i <= 20; i++) { double x = i / 5.0; // 標準の対数関数 double d1 = Math.Log(x); // 自作の対数関数 double x2 = (x - 1) / (x + 1); double d2 = 2 * myLog(x2, x2, 1.0, x2); // 標準関数との差異 Console.WriteLine(string.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)); } } // 自作の対数関数 private static double myLog(double x2, double numerator, double denominator, double y) { denominator = denominator + 2; numerator = numerator * x2 * x2; double a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.Abs(a) <= 0.00000000001) return y; else return y + myLog(x2, numerator, denominator, a); } }
using System; public class CS0507 { public static void Main() { for (int i = 1; i <= 20; i++) { double x = i / 5.0; // 標準の対数関数 double d1 = Math.Log(x); // 自作の対数関数 double d2 = myLog(x - 1, 27, 0.0); // 27:必要な精度が得られる十分大きな奇数 // 標準関数との差異 Console.WriteLine(string.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)); } } // 自作の対数関数 private static double myLog(double x, int n, double t) { int n2 = n; double x2 = x; if (n > 3) { if (n % 2 == 0) n2 = 2; x2 = x * (n / 2); } t = x2 / (n2 + t); if (n <= 2) return x / (1 + t); else return myLog(x, --n, t); } }
using System; public class CS0508 { public static void Main() { for (int x = -10; x <= 10; x++) { // 自作の双曲線正弦関数 double d1 = mySinh(x, 1, x, 1.0, x); // 標準の双曲線正弦関数 double d2 = Math.Sinh(x); // 標準関数との差異 Console.WriteLine(string.Format("{0,3:D} : {1,17:F10} - {2,17:F10} = {3,13:F10}", x, d1, d2, d1 - d2)); } } // 自作の双曲線正弦関数 private static double mySinh(double x, int n, double numerator, double denominator, double y) { int m = 2 * n; denominator = denominator * (m + 1) * m; numerator = numerator * x * x; double a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.Abs(a) <= 0.00000000001) return y; else return y + mySinh(x, ++n, numerator, denominator, a); } }
using System; public class CS0509 { public static void Main() { for (int x = -10; x <= 10; x++) { // 自作の双曲線余弦関数 double d1 = myCosh(x, 1, 1.0, 1.0, 1.0); // 標準の双曲線余弦関数 double d2 = Math.Cosh(x); // 標準関数との差異 Console.WriteLine(string.Format("{0,3:D} : {1,17:F10} - {2,17:F10} = {3,13:F10}", x, d1, d2, d1 - d2)); } } // 自作の双曲線余弦関数 private static double myCosh(double x, int n, double numerator, double denominator, double y) { int m = 2 * n; denominator = denominator * m * (m - 1); numerator = numerator * x * x; double a = numerator / denominator; // 十分な精度になったら処理を抜ける if (Math.Abs(a) <= 0.00000000001) return y; else return y + myCosh(x, ++n, numerator, denominator, a); } }
using System; public class CS0510 { public static void Main() { //for (int degree = -45; degree <= 45; degree += 15) for (int max = 1000; max <= 40000; max += 1000) { double radian = 1; // degree * Math.PI / 180.0; // 自作の逆正接関数 double d1 = myAtan(radian, false, radian, 1, radian, max); // 標準の逆正接関数 double d2 = Math.Atan(radian); // 標準関数との差異 //Console.WriteLine(string.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)); Console.WriteLine(string.Format("{0,5:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", max, d1 * 4, d2 * 4, (d1 * 4) - (d2 * 4))); } } // 自作の逆正接関数 private static double myAtan(double x, bool nega, double numerator, int denominator, double y, int max) { denominator = denominator + 2; numerator = numerator * x * x; double a = numerator / denominator; // 十分な精度になったら処理を抜ける //if (Math.Abs(a) <= 0.00000000001) if (denominator >= max) return y; else return y + myAtan(x, !nega, numerator, denominator, nega ? a : -a, max); } } /* 1000 : 3.1395926556 - 3.1415926536 = -0.0019999980 2000 : 3.1405926538 - 3.1415926536 = -0.0009999998 3000 : 3.1409259870 - 3.1415926536 = -0.0006666666 4000 : 3.1410926536 - 3.1415926536 = -0.0005000000 5000 : 3.1411926536 - 3.1415926536 = -0.0004000000 6000 : 3.1412593203 - 3.1415926536 = -0.0003333333 7000 : 3.1413069393 - 3.1415926536 = -0.0002857143 8000 : 3.1413426536 - 3.1415926536 = -0.0002500000 9000 : 3.1413704314 - 3.1415926536 = -0.0002222222 10000 : 3.1413926536 - 3.1415926536 = -0.0002000000 11000 : 3.1414108354 - 3.1415926536 = -0.0001818182 12000 : 3.1414259869 - 3.1415926536 = -0.0001666667 13000 : 3.1414388074 - 3.1415926536 = -0.0001538462 14000 : 3.1414497964 - 3.1415926536 = -0.0001428571 15000 : 3.1414593203 - 3.1415926536 = -0.0001333333 16000 : 3.1414676536 - 3.1415926536 = -0.0001250000 17000 : 3.1414750065 - 3.1415926536 = -0.0001176471 18000 : 3.1414815425 - 3.1415926536 = -0.0001111111 19000 : 3.1414873904 - 3.1415926536 = -0.0001052632 20000 : 3.1414926536 - 3.1415926536 = -0.0001000000 21000 : 3.1414974155 - 3.1415926536 = -0.0000952381 22000 : 3.1415017445 - 3.1415926536 = -0.0000909091 23000 : 3.1415056971 - 3.1415926536 = -0.0000869565 24000 : 3.1415093203 - 3.1415926536 = -0.0000833333 25000 : 3.1415126536 - 3.1415926536 = -0.0000800000 26000 : 3.1415157305 - 3.1415926536 = -0.0000769231 27000 : 3.1415185795 - 3.1415926536 = -0.0000740741 28000 : 3.1415212250 - 3.1415926536 = -0.0000714286 29000 : 3.1415236881 - 3.1415926536 = -0.0000689655 30000 : 3.1415259869 - 3.1415926536 = -0.0000666667 31000 : 3.1415281375 - 3.1415926536 = -0.0000645161 32000 : 3.1415301536 - 3.1415926536 = -0.0000625000 33000 : 3.1415320475 - 3.1415926536 = -0.0000606061 34000 : 3.1415338301 - 3.1415926536 = -0.0000588235 35000 : 3.1415355107 - 3.1415926536 = -0.0000571429 36000 : 3.1415370980 - 3.1415926536 = -0.0000555556 37000 : 3.1415385995 - 3.1415926536 = -0.0000540541 38000 : 3.1415400220 - 3.1415926536 = -0.0000526316 39000 : 3.1415413715 - 3.1415926536 = -0.0000512821 40000 : 3.1415426536 - 3.1415926536 = -0.0000500000 */
using System; public class CS0511 { public static void Main() { for (int degree = -45; degree <= 45; degree += 15) { double radian = degree * Math.PI / 180.0; double x2 = radian * radian; // 自作の逆正接関数 double d1 = myAtan(radian, x2, 23, 0.0); // 23:必要な精度が得られる十分大きな奇数 // 標準の逆正接関数 double d2 = Math.Atan(radian); // 標準関数との差異 Console.WriteLine(string.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)); } } // 自作の逆正接関数 private static double myAtan(double x, double x2, int n, double t) { int m = n / 2; t = (m * m * x2) / (n + t); n -= 2; if (n <= 1) return x / (1 + t); else return myAtan(x, x2, n, t); } }
using System; public class CS0512 { public static void Main() { for (int i = 11; i <= 31; i += 2) { double radian = 1; // degree * Math.PI / 180.0; double x2 = radian * radian; // 自作の逆正接関数 double d1 = myAtan(radian, x2, i, 0.0); // i:必要な精度が得られる十分大きな奇数 // 標準の逆正接関数 //double d2 = Math.Atan(radian); // 標準関数との差異 Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", i, d1 * 4, d1 * 4 - Math.PI)); //Console.WriteLine(string.Format("{0,2D} : {1,13:F10}", i, d1)); //Console.WriteLine(string.Format("{0,2:D}", i)); //Console.WriteLine(string.Format("{0,13:F10}, {1,13:F10}, {2,13:F10}", d1 * 4, d2 * 4, Math.PI)); } } // 自作の逆正接関数 private static double myAtan(double x, double x2, int n, double t) { int m = n / 2; t = (m * m * x2) / (n + t); n -= 2; if (n <= 1) return x / (1 + t); else return myAtan(x, x2, n, t); } }
using System; public class CS0601 { private static double f(double x) { return 4 / (1 + x * x); } public static void Main() { const double a = 0; const double b = 1; // 台形則で積分 int n = 2; for (int j = 1; j <= 10; j++) { double h = (b - a) / n; double s = 0; double x = a; for (int i = 1; i <= n - 1; i++) { x += h; s += f(x); } s = h * ((f(a) + f(b)) / 2 + s); n *= 2; // 結果を π と比較 Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)); } } }
using System; public class CS0602 { private static double f(double x) { return 4 / (1 + x * x); } public static void Main() { const double a = 0; const double b = 1; // 中点則で積分 int n = 2; for (int j = 1; j <= 10; j++) { double h = (b - a) / n; double s = 0; double x = a + (h / 2); for (int i = 1; i <= n; i++) { s += f(x); x += h; } s *= h; n *= 2; // 結果を π と比較 Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)); } } }
using System; public class CS0603 { private static double f(double x) { return 4 / (1 + x * x); } public static void Main() { const double a = 0; const double b = 1; // Simpson則で積分 int n = 2; for (int j = 1; j <= 5; j++) { double h = (b - a) / n; double s2 = 0; double s4 = 0; double x = a + h; for (int 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; double s = (s2 + s4) * h / 3; n *= 2; // 結果を π と比較 Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)); } } }
using System; public class CS0604 { private static double f(double x) { return 4 / (1 + x * x); } public static void Main() { const double a = 0; const double b = 1; double[,] t = new double[7,7]; // 台形則で積分 int n = 2; for (int i = 1; i <= 6; i++) { double h = (b - a) / n; double s = 0; double x = a; for (int j = 1; j <= n - 1; j++) { x += h; s += f(x); } // 結果を保存 t[i,1] = h * ((f(a) + f(b)) / 2 + s); n *= 2; } // Richardsonの補外法 n = 4; for (int j = 2; j <= 6; j++) { for (int 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) { // 結果を π と比較 Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, t[i,j], t[i,j] - Math.PI)); } } n *= 4; } } }
using System; public class CS0701 { // データ点の数 private const int N = 7; public static void Main() { double[] x = new double[N]; double[] y = new double[N]; // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (int i = 0; i < N; i++) { double d = i * 1.5 - 4.5; x[i] = d; y[i] = f(d); } // 0.5刻みで 与えられていない値を補間 for (int i = 0; i <= 18; i++) { double d = i * 0.5 - 4.5; double d1 = f(d); double d2 = lagrange(d, x, y); // 元の関数と比較 Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d, d1, d2, d1 - d2)); } } // 元の関数 private static double f(double x) { return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2); } // Lagrange (ラグランジュ) 補間 private static double lagrange(double d, double[] x, double[] y) { double sum = 0; for (int i = 0; i < N; i++) { double prod = y[i]; for (int j = 0; j < N; j++) { if (j != i) prod *= (d - x[j]) / (x[i] - x[j]); } sum += prod; } return sum; } }
using System; public class CS0702 { // データ点の数 private const int N = 7; public static void Main() { double[] x = new double[N]; double[] y = new double[N]; // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (int i = 0; i < N; i++) { double d = i * 1.5 - 4.5; x[i] = d; y[i] = f(d); } // 0.5刻みで 与えられていない値を補間 for (int i = 0; i <= 18; i++) { double d = i * 0.5 - 4.5; double d1 = f(d); double d2 = neville(d, x, y); // 元の関数と比較 Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d, d1, d2, d1 - d2)); } } // 元の関数 private static double f(double x) { return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2); } // Neville (ネヴィル) 補間 private static double neville(double d, double[] x, double[] y) { double[,] w = new double[N,N]; for (int i = 0; i < N; i++) w[0,i] = y[i]; for (int j = 1; j < N; j++) { for (int 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]; } }
using System; public class CS0703 { // データ点の数 private const int N = 7; public static void Main() { double[] x = new double[N]; double[] y = new double[N]; // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (int i = 0; i < N; i++) { double d1 = i * 1.5 - 4.5; x[i] = d1; y[i] = f(d1); } // 差分商の表を作る double[,] d = new double[N,N]; for (int j = 0; j < N; j++) d[0,j] = y[j]; for (int i = 1; i < N; i++) { for (int j = 0; j < N - i; j++) d[i,j] = (d[i-1,j+1] - d[i-1,j]) / (x[j+i] - x[j]); } // n階差分商 double[] a = new double[N]; for (int j = 0; j < N; j++) a[j] = d[j,0]; // 0.5刻みで 与えられていない値を補間 for (int i = 0; i <= 18; i++) { double d1 = i * 0.5 - 4.5; double d2 = f(d1); double d3 = newton(d1, x, a); // 元の関数と比較 Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d1, d2, d3, d2 - d3)); } } // 元の関数 private static double f(double x) { return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2); } // Newton (ニュートン) 補間 private static double newton(double d, double[] x, double[] a) { double sum = a[0]; for (int i = 1; i < N; i++) { double prod = a[i]; for (int j = 0; j < i; j++) prod *= (d - x[j]); sum += prod; } return sum; } }
using System; public class CS0704 { // データ点の数 private const int N = 7; private const int Nx2 = 14; public static void Main() { double[] x = new double[N]; double[] y = new double[N]; double[] yd = new double[N]; // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (int i = 0; i < N; i++) { double d1 = i * 1.5 - 4.5; x[i] = d1; y[i] = f(d1); yd[i] = fd(d1); } // 差分商の表を作る double[] z = new double[Nx2]; double[,] d = new double[Nx2,Nx2]; for (int i = 0; i < Nx2; i++) { int j = i / 2; z[i] = x[j]; d[0,i] = y[j]; } for (int i = 1; i < Nx2; i++) { for (int j = 0; j < Nx2 - i; j++) { if (i == 1 && j % 2 == 0) d[i,j] = yd[j / 2]; else d[i,j] = (d[i-1,j+1] - d[i-1,j]) / (z[j+i] - z[j]); } } // n階差分商 double[] a = new double[Nx2]; for (int j = 0; j < Nx2; j++) a[j] = d[j,0]; // 0.5刻みで 与えられていない値を補間 for (int i = 0; i <= 18; i++) { double d1 = i * 0.5 - 4.5; double d2 = f(d1); double d3 = hermite(d1, z, a); // 元の関数と比較 Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d1, d2, d3, d2 - d3)); } } // 元の関数 private static double f(double x) { return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2); } // 導関数 private static double fd(double x) { return 1 - Math.Pow(x,2) / 2 + Math.Pow(x,4) / (4 * 3 * 2); } // Hermite (エルミート) 補間 private static double hermite(double d, double[] z, double[] a) { double sum = a[0]; for (int i = 1; i < Nx2; i++) { double prod = a[i]; for (int j = 0; j < i; j++) prod *= (d - z[j]); sum += prod; } return sum; } }
using System; public class CS0705 { // データ点の数 private const int N = 7; public static void Main() { double[] x = new double[N]; double[] y = new double[N]; // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for (int i = 0; i < N; i++) { double d1 = i * 1.5 - 4.5; x[i] = d1; y[i] = f(d1); } // 3項方程式の係数の表を作る double[] a = new double[N]; double[] b = new double[N]; double[] c = new double[N]; double[] d = new double[N]; for (int 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項方程式を解く (ト−マス法) double[] g = new double[N]; double[] s = new double[N]; g[1] = b[1]; s[1] = d[1]; for (int 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]; } double[] z = new double[N]; z[0] = 0; z[N-1] = 0; z[N-2] = s[N-2] / g[N-2]; for (int i = N - 3; i >= 1; i--) { z[i] = (s[i] - c[i] * z[i+1]) / g[i]; } // 0.5刻みで 与えられていない値を補間 for (int i = 0; i <= 18; i++) { double d1 = i * 0.5 - 4.5; double d2 = f(d1); double d3 = spline(d1, x, y, z); // 元の関数と比較 Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d1, d2, d3, d2 - d3)); } } // 元の関数 private static double f(double x) { return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2); } // Spline (スプライン) 補間 private static double spline(double d, double[] x, double[] y, double[] z) { // 補間関数値がどの区間にあるか int k = -1; for (int i = 1; i < N; i++) { if (d <= x[i]) { k = i - 1; break; } } if (k < 0) k = N - 1; double d1 = x[k+1] - d; double d2 = d - x[k]; double 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; } }
using System; public class CS0801 { // 重力加速度 private const double g = -9.8; // 空気抵抗係数 private const double k = -0.01; // 時間間隔(秒) private const double h = 0.01; public static void Main() { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[2]; vx[0] = v * Math.Cos(radian); // 鉛直方向の速度 double[] vy = new double[2]; vy[0] = v * Math.Sin(radian); // 経過秒数 double t = 0.0; // 位置 double x = 0.0; double y = 0.0; // Euler法 for (int i = 1; y >= 0.0; i++) { // 経過秒数 t = i * h; // 位置 x += h * vx[0]; y += h * vy[0]; Console.WriteLine(string.Format("{0,4:F2}\t{1,8:F5}\t{2,9:F5}\t{3,9:F5}\t{4,8:F5}", t, vx[0], vy[0], x, y)); // 速度 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]; } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.Sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy); } }
using System; public class CS0802 { // 重力加速度 private const double g = -9.8; // 空気抵抗係数 private const double k = -0.01; // 時間間隔(秒) private const double h = 0.01; public static void Main() { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[3]; vx[0] = v * Math.Cos(radian); // 鉛直方向の速度 double[] vy = new double[3]; vy[0] = v * Math.Sin(radian); // 経過秒数 double t = 0.0; // 位置 double[] x = new double[3]; x[0] = 0.0; double[] y = new double[3]; y[0] = 0.0; // Heun法 for (int 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]; Console.WriteLine(string.Format("{0,4:F2}\t{1,8:F5}\t{2,9:F5}\t{3,9:F5}\t{4,8:F5}", t, vx[0], vy[0], x[0], y[0])); } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.Sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy); } }
using System; public class CS0803 { // 重力加速度 private const double g = -9.8; // 空気抵抗係数 private const double k = -0.01; // 時間間隔(秒) private const double h = 0.01; public static void Main() { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[2]; vx[0] = v * Math.Cos(radian); // 鉛直方向の速度 double[] vy = new double[2]; vy[0] = v * Math.Sin(radian); // 経過秒数 double t = 0.0; // 位置 double[] x = new double[2]; x[0] = 0.0; double[] y = new double[2]; y[0] = 0.0; // 中点法 for (int 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]); double wx = vx[0] + vx[1] / 2; double wy = vy[0] + vy[1] / 2; 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; Console.WriteLine(string.Format("{0,4:F2}\t{1,8:F5}\t{2,9:F5}\t{3,9:F5}\t{4,8:F5}", t, vx[0], vy[0], x[0], y[0])); } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.Sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy); } }
using System; public class CS0804 { // 重力加速度 private const double g = -9.8; // 空気抵抗係数 private const double k = -0.01; // 時間間隔(秒) private const double h = 0.01; public static void Main() { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[5]; vx[0] = v * Math.Cos(radian); // 鉛直方向の速度 double[] vy = new double[5]; vy[0] = v * Math.Sin(radian); // 経過秒数 double t = 0.0; // 位置 double[] x = new double[5]; x[0] = 0.0; double[] y = new double[5]; y[0] = 0.0; // Runge-Kutta法 for (int 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]); double wx = vx[0] + vx[1] / 2; double 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; Console.WriteLine(string.Format("{0,4:F2}\t{1,8:F5}\t{2,9:F5}\t{3,9:F5}\t{4,8:F5}", t, vx[0], vy[0], x[0], y[0])); } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.Sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy); } }
using System; public class CS0805 { // 重力加速度 private const double g = -9.8; // 空気抵抗係数 private const double k = -0.01; // 時間間隔(秒) private const double h = 0.01; public static void Main() { // 角度 double degree = 45; double radian = degree * Math.PI / 180.0; // 初速 250 km/h -> 秒速に変換 double v = 250 * 1000 / 3600; // 水平方向の速度 double[] vx = new double[5]; vx[0] = v * Math.Cos(radian); // 鉛直方向の速度 double[] vy = new double[5]; vy[0] = v * Math.Sin(radian); // 経過秒数 double t = 0.0; // 位置 double[] x = new double[5]; x[0] = 0.0; double[] y = new double[5]; y[0] = 0.0; // Runge-Kutta-Gill法 for (int 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]); double wx = vx[0] + vx[1] / 2; double 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; Console.WriteLine(string.Format("{0,4:F2}\t{1,8:F5}\t{2,9:F5}\t{3,9:F5}\t{4,8:F5}", t, vx[0], vy[0], x[0], y[0])); } } // 空気抵抗による水平方向の減速分 private static double fx(double vx, double vy) { return k * Math.Sqrt(vx * vx + vy * vy) * vx; } // 重力と空気抵抗による鉛直方向の減速分 private static double fy(double vx, double vy) { return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy); } }
using System; public class CS0901 { public static void Main() { double a = 1; double b = 2; Console.WriteLine(string.Format("{0,12:F10}", bisection(a, b))); } private static double bisection(double a, double b) { double c; while (true) { // 区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2; Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", c, c - Math.Sqrt(2))); double 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; } private static double f(double x) { return x * x - 2; } }
using System; public class CS0902 { public static void Main() { double a = 1; double b = 2; Console.WriteLine(string.Format("{0,12:F10}", falseposition(a, b))); } private static double falseposition(double a, double b) { double c; while (true) { // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c c = (a * f(b) - b * f(a)) / (f(b) - f(a)); Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", c, c - Math.Sqrt(2))); double 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; } private static double f(double x) { return x * x - 2; } }
using System; public class CS0903 { public static void Main() { double x = 1; Console.WriteLine(string.Format("{0,12:F10}", iterative(x))); } private static double iterative(double x0) { double x1; while (true) { x1 = g(x0); Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", x1, x1 - Math.Sqrt(2))); if (Math.Abs(x1 - x0) < 0.0000000001) break; x0 = x1; } return x1; } private static double g(double x) { return (x / 2) + (1 / x); } }
using System; public class CS0904 { public static void Main() { double x = 2; Console.WriteLine(string.Format("{0,12:F10}", newton(x))); } private static double newton(double x0) { double x1; while (true) { x1 = x0 - (f0(x0) / f1(x0)); Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", x1, x1 - Math.Sqrt(2))); if (Math.Abs(x1 - x0) < 0.0000000001) break; x0 = x1; } return x1; } private static double f0(double x) { return x * x - 2; } private static double f1(double x) { return 2 * x; } }
using System; public class CS0905 { public static void Main() { double x = 2; Console.WriteLine(string.Format("{0,12:F10}", bailey(x))); } private static double bailey(double x0) { double x1; while (true) { x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))); Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", x1, x1 - Math.Sqrt(2))); if (Math.Abs(x1 - x0) < 0.0000000001) break; x0 = x1; } return x1; } private static double f0(double x) { return x * x - 2; } private static double f1(double x) { return 2 * x; } private static double f2(double x) { return 2; } }
using System; public class CS0906 { public static void Main() { double x0 = 1; double x1 = 2; Console.WriteLine(string.Format("{0,12:F10}", secant(x0, x1))); } private static double secant(double x0, double x1) { double x2; while (true) { x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)); Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", x2, x2 - Math.Sqrt(2))); if (Math.Abs(x2 - x1) < 0.0000000001) break; x0 = x1; x1 = x2; } return x2; } private static double f(double x) { return x * x - 2; } }
using System; public class CS1001 { private const int N = 4; public static void Main() { double[,] a = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}}; double[] b = {20,16,8,17}; double[] c = {0,0,0,0}; // ヤコビの反復法 jacobi(a,b,c); Console.WriteLine("X"); disp_vector(c); } // ヤコビの反復法 private static void jacobi(double[,] a, double[] b, double[] x0) { while (true) { double[] x1 = new double[N]; bool finish = true; for (int i = 0; i < N; i++) { x1[i] = 0; for (int 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 (int i = 0; i < N; i++) x0[i] = x1[i]; if (finish) return; disp_vector(x0); } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } }
using System; public class CS1002 { private const int N = 4; public static void Main() { double[,] a = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}}; double[] b = {20,16,8,17}; double[] c = {0,0,0,0}; // ガウス・ザイデル法 gauss(a,b,c); Console.WriteLine("X"); disp_vector(c); } // ガウス・ザイデル法 private static void gauss(double[,] a, double[] b, double[] x0) { while (true) { bool finish = true; for (int i = 0; i < N; i++) { double x1 = 0; for (int 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次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } }
using System; public class CS1003 { private const int N = 4; // ガウスの消去法 public static void Main() { double[,] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double[] b = {8,17,20,16}; // ピボット選択 pivoting(a,b); Console.WriteLine("pivoting"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 前進消去 forward_elimination(a,b); Console.WriteLine("forward elimination"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 後退代入 backward_substitution(a,b); Console.WriteLine("X"); disp_vector(b); } // 前進消去 private static void forward_elimination(double[,] a, double[] b) { for (int pivot = 0; pivot < N - 1; pivot++) { for (int row = pivot + 1; row < N; row++) { double s = a[row,pivot] / a[pivot,pivot]; for (int col = pivot; col < N; col++) a[row,col] -= a[pivot,col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 private static void backward_substitution(double[,] a, double[] b) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) b[row] -= a[row,col] * b[col]; b[row] /= a[row,row]; } } // ピボット選択 private static void pivoting(double[,] a, double[] b) { for(int pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (int 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) { double tmp; for (int 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次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { int i = 0; foreach (double col in matrix) { Console.Write(string.Format("{0,14:F10}\t", col)); if (++i >= N) { i = 0; Console.WriteLine(); } } } }
using System; public class CS1004 { private const int N = 4; // ガウス・ジョルダン法 public static void Main() { double[,] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double[] b = {8,17,20,16}; // ピボット選択 pivoting(a,b); Console.WriteLine("pivoting"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 前進消去 forward_elimination(a,b); Console.WriteLine("forward elimination"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 後退代入 backward_substitution(a,b); Console.WriteLine("X"); disp_vector(b); } // 前進消去 private static void forward_elimination(double[,] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) { for (int row = 0; row < N; row++) { if (row == pivot) continue; double s = a[row,pivot] / a[pivot,pivot]; for (int col = pivot; col < N; col++) a[row,col] -= a[pivot,col] * s; b[row] -= b[pivot] * s; } } } // 後退代入 private static void backward_substitution(double[,] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) b[pivot] /= a[pivot,pivot]; } // ピボット選択 private static void pivoting(double[,] a, double[] b) { for(int pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (int 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) { double tmp; for (int 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次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { int i = 0; foreach (double col in matrix) { Console.Write(string.Format("{0,14:F10}\t", col)); if (++i >= N) { i = 0; Console.WriteLine(); } } } }
using System; public class CS1005 { private const int N = 4; // LU分解 public static void Main() { double[,] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}; double[] b = {8,17,20,16}; // ピボット選択 pivoting(a,b); Console.WriteLine("pivoting"); Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 前進消去 forward_elimination(a,b); Console.WriteLine("LU"); disp_matrix(a); // Ly=b から y を求める (前進代入) double[] y = {0,0,0,0}; forward_substitution(a,b,y); Console.WriteLine("Y"); disp_vector(y); // Ux=y から x を求める (後退代入) double[] x = {0,0,0,0}; backward_substitution(a,y,x); Console.WriteLine("X"); disp_vector(x); } // 前進消去 private static void forward_elimination(double[,] a, double[] b) { for (int pivot = 0; pivot < N - 1; pivot++) { for (int row = pivot + 1; row < N; row++) { double s = a[row,pivot] / a[pivot,pivot]; for (int col = pivot; col < N; col++) a[row,col] -= a[pivot,col] * s; // これが 上三角行列 a[row,pivot] = s; // これが 下三角行列 // b[row] -= b[pivot] * s; // この値は変更しない } } } // 前進代入 private static void forward_substitution(double[,] a, double[] b, double[] y) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row,col] * y[col]; y[row] = b[row]; } } // 後退代入 private static void backward_substitution(double[,] a, double[] y, double[] x) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[row,col] * x[col]; x[row] = y[row] / a[row,row]; } } // ピボット選択 private static void pivoting(double[,] a, double[] b) { for(int pivot = 0; pivot < N; pivot++) { // 各列で 一番値が大きい行を 探す int max_row = pivot; double max_val = 0; for (int 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) { double tmp; for (int 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次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { int i = 0; foreach (double col in matrix) { Console.Write(string.Format("{0,14:F10}\t", col)); if (++i >= N) { i = 0; Console.WriteLine(); } } } }
using System; public class CS1006 { private const int N = 4; // コレスキー法 public static void Main() { double[,] a = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}}; double[] b = {34,68,96,125}; Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 前進消去 forward_elimination(a,b); Console.WriteLine("LL^T"); disp_matrix(a); // Ly=b から y を求める (前進代入) double[] y = {0,0,0,0}; forward_substitution(a,b,y); Console.WriteLine("Y"); disp_vector(y); // L^Tx=y から x を求める (後退代入) double[] x = {0,0,0,0}; backward_substitution(a,y,x); Console.WriteLine("X"); disp_vector(x); } // 前進消去 private static void forward_elimination(double[,] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) { double s = 0; for (int col = 0; col < pivot; col++) s += a[pivot,col] * a[pivot,col]; // ここで根号の中が負の値になると計算できない! a[pivot,pivot] = Math.Sqrt(a[pivot,pivot] - s); for (int row = pivot + 1; row < N; row++) { s = 0; for (int 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]; } } } // 前進代入 private static void forward_substitution(double[,] a, double[] b, double[] y) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row,col] * y[col]; y[row] = b[row] / a[row,row]; } } // 後退代入 private static void backward_substitution(double[,] a, double[] y, double[] x) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[col,row] * x[col]; x[row] = y[row] / a[row,row]; } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { int i = 0; foreach (double col in matrix) { Console.Write(string.Format("{0,14:F10}\t", col)); if (++i >= N) { i = 0; Console.WriteLine(); } } } }
using System; public class CS1007 { private const int N = 4; // 修正コレスキー法 public static void Main() { double[,] a = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}}; double[] b = {34,68,96,125}; Console.WriteLine("A"); disp_matrix(a); Console.WriteLine("B"); disp_vector(b); Console.WriteLine(); // 前進消去 forward_elimination(a,b); Console.WriteLine("LDL^T"); disp_matrix(a); // Ly=b から y を求める (前進代入) double[] y = {0,0,0,0}; forward_substitution(a,b,y); Console.WriteLine("Y"); disp_vector(y); // DL^Tx=y から x を求める (後退代入) double[] x = {0,0,0,0}; backward_substitution(a,y,x); Console.WriteLine("X"); disp_vector(x); } // 前進消去 private static void forward_elimination(double[,] a, double[] b) { for (int pivot = 0; pivot < N; pivot++) { double s; // pivot < k の場合 for (int col = 0; col < pivot; col++) { s = a[pivot,col]; for (int 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 (int k = 0; k < pivot; k++) s -= a[pivot,k] * a[pivot,k] * a[k,k]; a[pivot,pivot] = s; } } // 前進代入 private static void forward_substitution(double[,] a, double[] b, double[] y) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row,col] * y[col]; y[row] = b[row]; } } // 後退代入 private static void backward_substitution(double[,] a, double[] y, double[] x) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[col,row] * a[row,row] * x[col]; x[row] = y[row] / a[row,row]; } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { int i = 0; foreach (double col in matrix) { Console.Write(string.Format("{0,14:F10}\t", col)); if (++i >= N) { i = 0; Console.WriteLine(); } } } }
using System; public class CS1101 { private const int N = 4; // ベキ乗法で最大固有値を求める public static void Main() { double[,] 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}}; double[] x = {1.0 ,0.0 ,0.0 ,0.0}; // ベキ乗法 double lambda = power(a, x); Console.WriteLine(); Console.WriteLine("eigenvalue"); Console.WriteLine(string.Format("{0,14:F10}", lambda)); Console.WriteLine("eigenvector"); disp_vector(x); } // ベキ乗法 private static double power(double[,] a, double[] x0) { double lambda = 0.0; // 正規化 (ベクトル x0 の長さを1にする) normarize(x0); double e0 = 0.0; for (int i = 0; i < N; i++) e0 += x0[i]; for (int k = 1; k <= 200; k++) { // 1次元配列を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_vector(x0); // 行列の積 x1 = A × x0 double[] x1 = new double[N]; for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) x1[i] += a[i, j] * x0[j]; // 内積 double p0 = 0.0; double p1 = 0.0; for (int i = 0; i < N; i++) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p0 / p1; // 正規化 (ベクトル x1 の長さを1にする) normarize(x1); // 収束判定 double e1 = 0.0; for (int i = 0; i < N; i++) e1 += x1[i]; if (Math.Abs(e0 - e1) < 0.00000000001) break; for (int i = 0; i < N; i++) x0[i] = x1[i]; e0 = e1; } return lambda; } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 正規化 (ベクトルの長さを1にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.Sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
using System; public class CS1102 { private const int N = 4; // 逆ベキ乗法で最小固有値を求める public static void Main() { double[,] 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}}; double[] x = {1.0 ,0.0 ,0.0 ,0.0}; // LU分解 forward_elimination(a); // 逆ベキ乗法 double lambda = inverse(a, x); Console.WriteLine(); Console.WriteLine("eigenvalue"); Console.WriteLine(string.Format("{0,14:F10}", lambda)); Console.WriteLine("eigenvector"); disp_vector(x); } // 逆ベキ乗法 private static double inverse(double[,] a, double[] x0) { double lambda = 0.0; // 正規化 (ベクトル x0 の長さを1にする) normarize(x0); double e0 = 0.0; for (int i = 0; i < N; i++) e0 += x0[i]; for (int k = 1; k <= 200; k++) { // 1次元配列を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_vector(x0); // Ly = b から y を求める (前進代入) double[] b = new double[N]; double[] y = new double[N]; for (int i = 0; i < N; i++) b[i] = x0[i]; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) double[] x1 = new double[N]; backward_substitution(a,x1,y); // 内積 double p0 = 0.0; double p1 = 0.0; for (int i = 0; i < N; i++) { p0 += x1[i] * x1[i]; p1 += x1[i] * x0[i]; } // 固有値 lambda = p1 / p0; // 正規化 (ベクトル x1 の長さを1にする) normarize(x1); // 収束判定 double e1 = 0.0; for (int i = 0; i < N; i++) e1 += x1[i]; if (Math.Abs(e0 - e1) < 0.00000000001) break; for (int i = 0; i < N; i++) x0[i] = x1[i]; e0 = e1; } return lambda; } // LU分解 private static void forward_elimination(double[,] a) { for (int pivot = 0; pivot < N - 1; pivot++) { for (int row = pivot + 1; row < N; row++) { double s = a[row,pivot] / a[pivot,pivot]; for (int col = pivot; col < N; col++) a[row,col] -= a[pivot,col] * s; // これが 上三角行列 a[row,pivot] = s; // これが 下三角行列 } } } // 前進代入 private static void forward_substitution(double[,] a, double[] y, double[] b) { for (int row = 0; row < N; row++) { for (int col = 0; col < row; col++) b[row] -= a[row,col] * y[col]; y[row] = b[row]; } } // 後退代入 private static void backward_substitution(double[,] a, double[] x, double[] y) { for (int row = N - 1; row >= 0; row--) { for (int col = N - 1; col > row; col--) y[row] -= a[row,col] * x[col]; x[row] = y[row] / a[row,row]; } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 正規化 (ベクトルの長さを1にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.Sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
using System; public class CS1103 { private const int N = 4; // LR分解で固有値を求める public static void Main() { double[,] 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}}; double[,] l = new double[N,N]; double[,] u = new double[N,N]; for (int k = 1; k <= 200; k++) { // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_eigenvalue(a); // 収束判定 double e = 0.0; for (int i = 1; i < N; i++) for (int j = 0; j < i; j++) e += Math.Abs(a[i, j]); if (e < 0.00000000001) break; } Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_eigenvalue(a); } // LU分解 private static void decomp(double[,] a, double[,] l, double[,] u) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { l[i, j] = 0.0; u[i, j] = 0.0; } } l[0, 0] = 1.0; for (int j = 0; j < N; j++) u[0, j] = a[0, j]; for (int 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 (int i = 1; i < N; i++) { l[i, i] = 1.0; double t = a[i, i]; for (int k = 0; k <= i; k++) t -= l[i, k] * u[k, i]; u[i, i] = t; for (int j = i + 1; j < N; j++) { u[j, i] = 0.0; l[i, j] = 0.0; t = a[j, i]; for (int k = 0; k <= i; k++) t -= l[j, k] * u[k, i]; l[j, i] = t / u[i, i]; t = a[i, j]; for (int k = 0; k <= i; k++) t -= l[i, k] * u[k, j]; u[i, j] = t; } } } // 行列の積 private static void multiply(double[,] a, double[,] b, double[,] c) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double s = 0.0; for (int k = 0; k < N; k++) s += a[i, k] * b[k, j]; c[i,j] = s; } } } // 対角要素を表示 private static void disp_eigenvalue(double[,] a) { for (int i = 0; i < N; i++) Console.Write(string.Format("{0,14:F10}\t", a[i,i])); Console.WriteLine(); } }
using System; public class CS1104 { private const int N = 4; // QR分解で固有値を求める public static void Main() { double[,] 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}}; double[,] q = new double[N,N]; double[,] r = new double[N,N]; for (int k = 1; k <= 100; k++) { // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, a); // 対角要素を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_eigenvalue(a); // 収束判定 double e = 0.0; for (int i = 0; i < N; i++) for (int j = 0; j < i; j++) e += Math.Abs(a[i, j]); if (e < 0.00000000001) break; } Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_eigenvalue(a); } // QR分解 private static void decomp(double[,] a, double[,] q, double[,] r) { double[] x = new double[N]; for (int k = 0; k < N; k++) { for (int i = 0; i < N; i++) x[i] = a[i, k]; for (int j = 0; j < k; j++) { double t = 0.0; for (int i = 0; i < N; i++) t += a[i, k] * q[i, j]; r[j, k] = t; r[k, j] = 0.0; for (int i = 0; i < N; i++) x[i] -= t * q[i, j]; } double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; r[k,k] = Math.Sqrt(s); for (int i = 0; i < N; i++) q[i, k] = x[i] / r[k, k]; } } // 行列の積 private static void multiply(double[,] a, double[,] b, double[,] c) { for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double s = 0.0; for (int k = 0; k < N; k++) s += a[i, k] * b[k, j]; c[i,j] = s; } } } // 対角要素を表示 private static void disp_eigenvalue(double[,] a) { for (int i = 0; i < N; i++) Console.Write(string.Format("{0,14:F10}\t", a[i,i])); Console.WriteLine(); } }
using System; public class CS1105 { private const int N = 4; // ヤコビ法で固有値を求める public static void Main() { double[,] 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}}; double[,] 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); Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_eigenvalue(a); Console.WriteLine(); Console.WriteLine("eigenvector"); disp_eigenvector(v); } // ヤコビ法 private static void jacobi(double[,] a, double[,] v) { for (int k = 1; k <= 100; k++) { // 最大値を探す int p = 0; int q = 0; double max_val = 0.0; for (int i = 0; i < N; i++) { for (int 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; } } } // θ を求める double 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 double c = Math.Cos(t); double s = Math.Sin(t); // U^t × A double t1 = 0.0; double t2 = 0.0; for (int 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 (int 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; } // 対角要素を表示 Console.Write(string.Format("{0,3:D}\t", k)); disp_eigenvalue(a); // 収束判定 if (max_val < 0.00000000001) break; } } // 対角要素を表示 private static void disp_eigenvalue(double[,] a) { for (int i = 0; i < N; i++) Console.Write(string.Format("{0,14:F10}\t", a[i,i])); Console.WriteLine(); } // 固有ベクトルを表示 private static void disp_eigenvector(double[,] matrix) { for (int i = 0; i < N; i++) { double[] row = new double[N]; for (int j = 0; j < N; j++) row[j] = matrix[i, j] ; normarize(row); disp_vector(row); } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 正規化 (ベクトルの長さを1にする) private static void normarize(double[] x) { double s = 0.0; for (int i = 0; i < N; i++) s += x[i] * x[i]; s = Math.Sqrt(s); for (int i = 0; i < N; i++) x[i] /= s; } }
using System; public class CS1106 { private const int N = 4; // ハウスホルダー変換とQR分解で固有値を求める public static void Main() { double[,] 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}}; double[] d = new double[N]; double[] e = new double[N]; // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); Console.WriteLine(); Console.WriteLine("eigenvalue"); disp_vector(d); Console.WriteLine(); Console.WriteLine("eigenvector"); disp_matrix(a); } // ハウスホルダー変換 private static void tridiagonalize(double[,] a, double[] d, double[] e) { double[] w = new double[N]; double[] v = new double[N]; for (int k = 0; k < N - 2; k++) { d[k] = a[k, k]; double t = 0.0; for (int i = k + 1; i < N; i++) { w[i] = a[i, k]; t += w[i] * w[i]; } double 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 (int i = k + 1; i < N; i++) w[i] *= s; for (int i = k + 1; i < N; i++) { s = 0.0; for (int 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 (int i = k + 1; i < N; i++) s += w[i] * v[i]; s /= 2.0; for (int i = k + 1; i < N; i++) v[i] -= s * w[i]; for (int i = k + 1; i < N; i++) for (int j = k + 1; j <= i; j++) a[i, j] -= w[i] * v[j] + w[j] * v[i]; // 次の行は固有ベクトルを求めないなら不要 for (int 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 (int k = N - 1; k >= 0; k--) { w = new double[] {0.0 ,0.0 ,0.0 ,0.0}; if (k < N - 2) { for (int i = k + 1; i < N; i++) w[i] = a[i, k]; for (int i = k + 1; i < N; i++) { double s = 0.0; for (int j = k + 1; j < N; j++) s += a[i, j] * w[j]; v[i] = s; } for (int i = k + 1; i < N; i++) for (int j = k + 1; j < N; j++) a[i, j] -= v[i] * w[j]; } for (int i = 0; i < N; i++) a[i, k] = 0.0; a[k, k] = 1.0; } } // QR分解 private static void decomp(double[,] a, double[] d, double[] e) { e[0] = 1.0; int h = N - 1; while (Math.Abs(e[h]) < 0.00000000001) h--; while (h > 0) { e[0] = 0.0; int l = h - 1; while (Math.Abs(e[l]) >= 0.00000000001) l--; for (int j = 1; j <= 100; j++) { double w = (d[h - 1] - d[h]) / 2.0; double s = Math.Sqrt(w * w + e[h] * e[h]); if (w < 0.0) s = -s; double x = d[l] - d[h] + e[h] * e[h] / (w + s); double y = e[l + 1]; double z = 0.0; double t = 0.0; double u = 0.0; for (int 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 (int 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]; } } Console.Write(string.Format("{0,3:D}\t", j)); disp_vector(d); // 収束判定 if (Math.Abs(e[h]) < 0.00000000001) break; } e[0] = 1.0; while (Math.Abs(e[h]) < 0.00000000001) h--; } // 次の行は固有ベクトルを求めないなら不要 for (int k = 0; k < N - 1; k++) { int l = k; for (int i = k + 1; i < N; i++) if (d[i] > d[l]) l = i; double t = d[k]; d[k] = d[l]; d[l] = t; for (int i = 0; i < N; i++) { t = a[k, i]; a[k, i] = a[l, i]; a[l, i] = t; } } } // 1次元配列を表示 private static void disp_vector(double[] row) { foreach (double col in row) Console.Write(string.Format("{0,14:F10}\t", col)); Console.WriteLine(); } // 2次元配列を表示 private static void disp_matrix(double[,] matrix) { for (int row = 0; row < N; row++) { for (int col = 0; col < N; col++) Console.Write(string.Format("{0,14:F10}\t", matrix[row,col])); Console.WriteLine(); } } }