program Pas0101(arg); uses SysUtils, Math; begin writeln(3 + 5); writeln(3 - 5); writeln(3 * 5); writeln(power(3, 5)); writeln(5 / 3); writeln(5 div 3); writeln(5 mod 3); write(IntToStr(3 * 5) + #13#10); writeln(format('%3d', [3 * 5])); writeln(format('%23.20f', [5 / 3])); end.
program Pas0102(arg); uses SysUtils; var i:Integer = 3 * 5; begin writeln(format('3 * 5 = %d', [i])); end.
program Pas0103(arg); uses SysUtils; var i:Integer; begin for i := 1 to 9 do begin write(format('%d, ', [i])); end; writeln(); end.
program Pas0104(arg); uses SysUtils; var i:Integer; begin for i := 1 to 9 do begin if (i mod 3 = 0) then write(format('%d, ', [i])); end; writeln(); end.
program Pas0105(arg); var i:Integer; sum:Integer = 0; begin for i := 1 to 99 do begin if (i mod 3 = 0) then sum := sum + i; end; writeln(sum); end.
program Pas0301(arg); {$MODE delphi} // 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 function sn(a:Integer; lim:Integer):Integer; var n, l:Integer; begin n := lim div a; // 項数:n = 上限:lim / 公差:a l := n * a; // 末項:l = 項数:n * 公差:a result := (a + l) * n div 2; // 総和:sn = (初項:a + 末項:l) * 項数:n / 2 end; begin // 3 の倍数の合計 writeln(sn(3, 999)); end.
program Pas0302(arg); var n:Integer; begin // 10000 までの 自然数の和 // 項目数 n = 10000 n := 10000; writeln( n * (n + 1) div 2 ); end.
program Pas0303(arg); var n:Integer; begin // 10000 までの 偶数の和 // 項目数 n = 5000 n := 10000 div 2; writeln( n * (n + 1) ); end.
program Pas0304(arg); uses SysUtils, Math; var n:Integer; begin // 10000 までの 奇数の和 // 項目数 n = 5000 n := 10000 div 2; writeln( format('%g', [power(n, 2)]) ); // power は Extended 型 end.
program Pas0305(arg); var n:Integer; begin // 1000 までの 自然数の2乗の和 n := 1000; writeln( n * (n + 1) * (2 * n + 1) div 6 ); end.
program Pas0306(arg); uses SysUtils, Math; var n:Integer; begin // 100 までの 自然数の3乗の和 n := 100; writeln( format('%g', [(power(n, 2) * power((n + 1), 2) / 4)]) ); end.
program Pas0303(arg); uses SysUtils, Math; var n:Integer; a:Integer; r:Integer; begin // 初項 2, 公比 3, 項数 10 の等比数列の和 n := 10; a := 2; r := 3; writeln( format('%g', [(a * (power(r, n) - 1)) / (r - 1)]) ); end.
program Pas0401(arg); const a:Integer = 5; // 初項 5 d:Integer = 3; // 公差 3 n:Integer = 10; // 項数 10 var p:Int64 = 1; // 積 m:Integer; i:Integer; begin for i := 1 to n do begin m := a + (d * (i - 1)); p := p * m; end; writeln(p); end.
program Pas0402(arg); {$MODE delphi} function product(m:Integer; d:Integer; n:Integer):Int64; begin if n = 0 then result := 1 else result := m * product(m + d, d, n - 1); end; begin // 初項 5, 公差 3, 項数 10 の数列の積を表示する writeln(product(5, 3, 10)); end.
program Pas0403(arg); // 階乗を求める関数 function Fact(n: Integer): Longint; begin if n <= 1 then Fact := 1 else Fact := n * Fact(n - 1); end; begin // 10の階乗 writeln(Fact(10)); writeln(10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1); end.
program Pas0404(arg); // 下降階乗冪 function FallingFact(x: Integer; n: Integer): Longint; begin if n <= 1 then FallingFact := x else FallingFact := x * FallingFact(x - 1, n - 1); end; begin // 10 から 6 までの 総乗 writeln(FallingFact(10, 5)); writeln(10 * 9 * 8 * 7 * 6); end.
program Pas0405(arg); // 上昇階乗冪 function RisingFact(x: Integer; n: Integer): Longint; begin if n <= 1 then RisingFact := x else RisingFact := x * RisingFact(x + 1, n - 1); end; begin // 10 から 14 までの 総乗 writeln(RisingFact(10, 5)); writeln(10 * 11 * 12 * 13 * 14); end.
program Pas0406(arg); uses SysUtils; // 階乗 function Fact(n: Integer): Longint; begin if n <= 1 then Fact := 1 else Fact := n * Fact(n - 1); end; // 下降階乗冪 function FallingFact(x: Integer; n: Integer): Longint; begin if n <= 1 then FallingFact := x else FallingFact := x * FallingFact(x - 1, n - 1); end; var n: Integer; r: Integer; begin // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) n := 10; r := 5; writeln(format('%g', [Fact(n) / Fact(n - r)])); writeln(FallingFact(n, r)); end.
program Pas0407(arg); uses SysUtils, Math; var n: Integer; r: Integer; begin // 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) n := 10; r := 5; writeln(format('%g', [power(n, r)])); end.
program Pas040101(arg); // 組合せ function Comb(n: Integer; r: Integer): Longint; begin if (r = 0) or (r = n) then Comb := 1 else if r = 1 then Comb := n else Comb := Comb(n - 1, r - 1) + Comb(n - 1, r); end; var n: Integer; r: Integer; begin // 組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) n := 10; r := 5; writeln(Comb(n, r)); end.
program Pas0409(arg); // 組合せ function Comb(n: Integer; r: Integer): Longint; begin if (r = 0) or (r = n) then Comb := 1 else if r = 1 then Comb := n else Comb := Comb(n - 1, r - 1) + Comb(n - 1, r); end; var n: Integer; r: Integer; begin // 重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) n := 10; r := 5; writeln(Comb(n + r - 1, r)); end.
program Pas0501(arg); {$MODE delphi} uses SysUtils, Math; // 自作の正弦関数 function mySin(x:Double; n:Integer; nega:Boolean; numerator:Double; denominator:Double; y:Double):Double; var m: Integer; a: Double; begin m := 2 * n; denominator := denominator * (m + 1) * m; numerator := numerator * x * x; a := numerator / denominator; // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) then result := y else begin if not nega then a := -a; result := y + mySin(x, n + 1, not nega, numerator, denominator, a); end; end; var i: Integer; degree: Integer; radian: Double; d1: Double; d2: Double; begin for i := 0 to 24 do begin degree := i * 15; if (degree mod 30 = 0) or (degree mod 45 = 0) then begin radian := DegToRad(degree); // 自作の正弦関数 d1 := mySin(radian, 1, false, radian, 1.0, radian); // 標準の正弦関数 d2 := Sin(radian); // 標準関数との差異 writeln(format('%3d : %13.10f - %13.10f = %13.10f', [degree, d1, d2, d1 - d2])); end; end; end.
program Pas0502(arg); {$MODE delphi} uses SysUtils, Math; // 自作の余弦関数 function myCos(x:Double; n:Integer; nega:Boolean; numerator:Double; denominator:Double; y:Double):Double; var m: Integer; a: Double; begin m := 2 * n; denominator := denominator * m * (m - 1); numerator := numerator * x * x; a := numerator / denominator; // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) then result := y else begin if not nega then a := -a; result := y + myCos(x, n + 1, not nega, numerator, denominator, a); end; end; var i: Integer; degree: Integer; radian: Double; d1: Double; d2: Double; begin for i := 0 to 24 do begin degree := i * 15; if (degree mod 30 = 0) or (degree mod 45 = 0) then begin radian := DegToRad(degree); // 自作の余弦関数 d1 := myCos(radian, 1, false, 1.0, 1.0, 1.0); // 標準の余弦関数 d2 := Cos(radian); // 標準関数との差異 writeln(format('%3d : %13.10f - %13.10f = %13.10f', [degree, d1, d2, d1 - d2])); end; end; end.
program Pas0503(arg); {$MODE delphi} uses SysUtils, Math; // 自作の正接関数 function myTan(x:Double; x2:Double; n:Integer; t:Double):Double; begin t := x2 / (n - t); n := n - 2; if (n <= 1) then result := x / (1 - t) else result := myTan(x, x2, n, t); end; var i: Integer; degree: Integer; radian: Double; x2: Double; d1: Double; d2: Double; begin for i := 0 to 12 do begin if (i * 15 mod 180 <> 0) then begin degree := i * 15 - 90; radian := DegToRad(degree); // 自作の正接関数 x2 := radian * radian; d1 := myTan(radian, x2, 15, 0.0); // 15:必要な精度が得られる十分大きな奇数 // 標準の正接関数 d2 := Tan(radian); // 標準関数との差異 writeln(format('%3d : %13.10f - %13.10f = %13.10f', [degree, d1, d2, d1 - d2])); end; end; end.
program Pas0504(arg); {$MODE delphi} uses SysUtils, Math; // 自作の指数関数 function myExp(x:Double; n:Integer; numerator:Double; denominator:Double; y:Double):Double; var a: Double; begin denominator := denominator * n; numerator := numerator * x; a := numerator / denominator; // 十分な精度になったら処理を抜ける if (Abs(a) <= 0.00000000001) then result := y else result := y + myExp(x, n + 1, numerator, denominator, a); end; var i: Integer; x: Double; d1: Double; d2: Double; begin for i := 0 to 20 do begin x := (i - 10) / 4.0; // 標準の指数関数 d1 := Exp(x); // 自作の指数関数 d2 := myExp(x, 1, 1.0, 1.0, 1.0); // 標準関数との差異 writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2])); end; end.
program Pas0505(arg); {$MODE delphi} uses SysUtils, Math; // 自作の指数関数 function myExp(x:Double; x2:Double; n:Integer; t:Double):Double; begin t := x2 / (n + t); n := n - 4; if (n < 6) then result := 1 + ((2 * x) / (2 - x + t)) else result := myExp(x, x2, n, t); end; var i: Integer; x: Double; x2: Double; d1: Double; d2: Double; begin for i := 0 to 20 do begin x := (i - 10) / 4.0; // 標準の指数関数 d1 := Exp(x); // 自作の指数関数 x2 := x * x; d2 := myExp(x, x2, 30, 0.0); // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる // 標準関数との差異 writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2])); end; end.
program Pas0506(arg); {$MODE delphi} uses SysUtils, Math; // 自作の対数関数 function myLog(x2:Double; numerator:Double; denominator:Double; y:Double):Double; var a :Double; begin denominator := denominator + 2; numerator := numerator * x2 * x2; a := numerator / denominator; // 十分な精度になったら処理を抜ける if (Abs(a) <= 0.00000000001) then result := y else result := y + myLog(x2, numerator, denominator, a); end; var i: Integer; x: Double; x2: Double; d1: Double; d2: Double; begin for i := 1 to 20 do begin x := i / 5.0; // 標準の対数関数 d1 := Ln(x); // 自作の対数関数 x2 := (x - 1) / (x + 1); d2 := 2 * myLog(x2, x2, 1.0, x2); // 標準関数との差異 writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2])); end; end.
program Pas0507(arg); {$MODE delphi} uses SysUtils, Math; // 自作の対数関数 function myLog(x:Double; n:Integer; t:Double):Double; var n2 :Integer; x2 :Double; begin n2 := n; x2 := x; if (n > 3) then begin if (n mod 2 = 0) then n2 := 2; x2 := x * (n div 2); end; t := x2 / (n2 + t); if (n <= 2) then result := x / (1 + t) else result := myLog(x, n - 1, t); end; var i: Integer; x: Double; d1: Double; d2: Double; begin for i := 1 to 20 do begin x := i / 5.0; // 標準の対数関数 d1 := Ln(x); // 自作の対数関数 d2 := myLog(x - 1, 27, 0.0); // 27:必要な精度が得られる十分大きな奇数 // 標準関数との差異 writeln(format('%5.2f : %13.10f - %13.10f = %13.10f', [x, d1, d2, d1 - d2])); end; end.
program Pas0508(arg); {$MODE delphi} uses SysUtils, Math; // 自作の双曲線正弦関数 function mySinh(x:Double; n:Integer; numerator:Double; denominator:Double; y:Double):Double; var m: Integer; a: Double; begin m := 2 * n; denominator := denominator * (m + 1) * m; numerator := numerator * x * x; a := numerator / denominator; // 十分な精度になったら処理を抜ける if (Abs(a) <= 0.00000000001) then result := y else result := y + mySinh(x, n + 1, numerator, denominator, a); end; var i: Integer; x: Integer; d1: Double; d2: Double; begin for i := 0 to 20 do begin x := i - 10; // 自作の双曲線正弦関数 d1 := mySinh(x, 1, x, 1.0, x); // 標準の双曲線正弦関数 d2 := Sinh(x); // 標準関数との差異 writeln(format('%3d : %17.10f - %17.10f = %13.10f', [x, d1, d2, d1 - d2])); end; end.
program Pas0509(arg); {$MODE delphi} uses SysUtils, Math; // 自作の双曲線余弦関数 function myCosh(x:Double; n:Integer; numerator:Double; denominator:Double; y:Double):Double; var m: Integer; a: Double; begin m := 2 * n; denominator := denominator * m * (m - 1); numerator := numerator * x * x; a := numerator / denominator; // 十分な精度になったら処理を抜ける if (Abs(a) <= 0.00000000001) then result := y else result := y + myCosh(x, n + 1, numerator, denominator, a); end; var i: Integer; x: Integer; d1: Double; d2: Double; begin for i := 0 to 20 do begin x := i - 10; // 自作の双曲線余弦関数 d1 := myCosh(x, 1, 1.0, 1.0, 1.0); // 標準の双曲線余弦関数 d2 := Cosh(x); // 標準関数との差異 writeln(format('%3d : %17.10f - %17.10f = %13.10f', [x, d1, d2, d1 - d2])); end; end.
program Pas0601(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := 4 / (1 + x * x); end; const a:Double = 0; b:Double = 1; var n, i, j:Integer; h, s, x:Double; begin // 台形則で積分 n := 2; for j := 1 to 10 do begin h := (b - a) / n; s := 0; x := a; for i := 1 to n - 1 do begin x := x + h; s := s + f(x); end; s := h * ((f(a) + f(b)) / 2 + s); n := n * 2; // 結果を π と比較 writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI])); end; end.
program Pas0602(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := 4 / (1 + x * x); end; const a:Double = 0; b:Double = 1; var n, i, j:Integer; h, s, x:Double; begin // 中点則で積分 n := 2; for j := 1 to 10 do begin h := (b - a) / n; s := 0; x := a + (h / 2); for i := 1 to n do begin s := s + f(x); x := x + h; end; s := h * s; n := n * 2; // 結果を π と比較 writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI])); end; end.
program Pas0603(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := 4 / (1 + x * x); end; const a:Double = 0; b:Double = 1; var n, i, j:Integer; h, s, x:Double; s2, s4 :Double; begin // Simpson則で積分 n := 2; for j := 1 to 5 do begin h := (b - a) / n; s2 := 0; s4 := 0; x := a + h; for i := 1 to (n div 2) do begin s4 := s4 + f(x); x := x + h; s2 := s2 + f(x); x := x + h; end; s2 := (s2 - f(b)) * 2 + f(a) + f(b); s4 := s4 * 4; s := (s2 + s4) * h / 3; n := n * 2; // 結果を π と比較 writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI])); end; end.
program Pas0604(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := 4 / (1 + x * x); end; const a:Double = 0; b:Double = 1; var n, i, j:Integer; h, s, x:Double; t:array [1..6, 1..6] of Double; begin // 台形則で積分 n := 2; for i := 1 to 6 do begin h := (b - a) / n; s := 0; x := a; for j := 1 to n - 1 do begin x := x + h; s := s + f(x); end; // 結果を保存 t[i,1] := h * ((f(a) + f(b)) / 2 + s); n := n * 2; end; // Richardsonの補外法 n := 4; for j := 2 to 6 do begin for i := j to 6 do begin t[i,j] := t[i, j - 1] + (t[i, j - 1] - t[i - 1, j - 1]) / (n - 1); if i = j then begin // 結果を π と比較 writeln(format('%2d : %13.10f, %13.10f', [j, t[i,j], t[i,j] - PI])); end; end; n := n * 4; end; end.
program Pas0701(arg); {$MODE delphi} uses SysUtils, Math; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // Lagrange (ラグランジュ) 補間 function lagrange(d:Double; x:array of Double; y:array of Double):Double; var sum, prod :Double; i, j :Integer; begin sum := 0; for i := Low(x) to High(x) do begin prod := y[i]; for j := Low(x) to High(x) do begin if j <> i then prod := prod * (d - x[j]) / (x[i] - x[j]); end; sum := sum + prod; end; result := sum; end; const // データ点の数 - 1 N = 6; var i :Integer; x :array [0..N] of Double; y :array [0..N] of Double; d, d1, d2 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := Low(x) to High(x) do begin d := (i - 1) * 1.5 - 4.5; x[i] := d; y[i] := f(d); end; // 0.5刻みで 与えられていない値を補間 for i := 0 to 18 do begin d := i * 0.5 - 4.5; d1 := f(d); d2 := lagrange(d, x, y); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d, d1, d2, d1 - d2])); end; end.
program Pas0702(arg); {$MODE delphi} uses SysUtils, Math; const // データ点の数 - 1 N = 6; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // Neville (ネヴィル) 補間 function neville(d:Double; x:array of Double; y:array of Double):Double; var w :array [0..N, 0..N] of Double; i, j :Integer; begin for i := Low(y) to High(y) do w[0,i] := y[i]; for j := 1 to High(x) do begin for i := Low(x) to (High(x) - j) do 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]); end; result := w[N,0]; end; var i :Integer; x :array [0..N] of Double; y :array [0..N] of Double; d, d1, d2 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := Low(x) to High(x) do begin d := i * 1.5 - 4.5; x[i] := d; y[i] := f(d); end; // 0.5刻みで 与えられていない値を補間 for i := 0 to 18 do begin d := i * 0.5 - 4.5; d1 := f(d); d2 := neville(d, x, y); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d, d1, d2, d1 - d2])); end; end.
program Pas0703(arg); {$MODE delphi} uses SysUtils, Math; const // データ点の数 - 1 N = 6; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // Newton (ニュートン) 補間 function newton(d:Double; x:array of Double; a:array of Double):Double; var sum, prod :Double; i, j :Integer; begin sum := a[0]; for i := 1 to High(x) do begin prod := a[i]; for j := Low(x) to i-1 do prod := prod * (d - x[j]); sum := sum + prod; end; result := sum; end; var i, j :Integer; x :array [0..N] of Double; y :array [0..N] of Double; a :array [0..N] of Double; d :array [0..N, 0..N] of Double; d1, d2, d3 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := Low(x) to High(x) do begin d1 := i * 1.5 - 4.5; x[i] := d1; y[i] := f(d1); end; // 差分商の表を作る for j := Low(y) to High(y) do d[0,j] := y[j]; for i := 1 to High(x) do begin for j := Low(x) to High(x)-i do d[i,j] := (d[i-1,j+1] - d[i-1,j]) / (x[j+i] - x[j]); end; // n階差分商 for j := Low(a) to High(a) do a[j] := d[j,0]; // 0.5刻みで 与えられていない値を補間 for i := 0 to 18 do begin d1 := i * 0.5 - 4.5; d2 := f(d1); d3 := newton(d1, x, a); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3])); end; end.
program Pas0704(arg); {$MODE delphi} uses SysUtils, Math; const // データ点の数 - 1 N = 6; Nx2 = 13; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // 導関数 function fd(x:Double):Double; begin result := 1 - power(x,2) / 2 + power(x,4) / (4 * 3 * 2); end; // Hermite (エルミート) 補間 function hermite(d:Double; z:array of Double; a:array of Double):Double; var sum, prod :Double; i, j :Integer; begin sum := a[0]; for i := 1 to High(a) do begin prod := a[i]; for j := Low(z) to i-1 do prod := prod * (d - z[j]); sum := sum + prod; end; result := sum; end; var i, j :Integer; x :array [0..N] of Double; y :array [0..N] of Double; yd :array [0..N] of Double; z :array [0..Nx2] of Double; a :array [0..Nx2] of Double; d :array [0..Nx2, 0..Nx2] of Double; d1, d2, d3 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := Low(x) to High(x) do begin d1 := i * 1.5 - 4.5; x[i] := d1; y[i] := f(d1); yd[i] := fd(d1); end; // 差分商の表を作る for i := Low(z) to High(z) do begin j := i div 2; z[i] := x[j]; d[0,i] := y[j]; end; for i := 1 to High(z) do begin for j := Low(z) to High(z)-i do begin if (i = 1) and (j mod 2 = 0) then d[i,j] := yd[j div 2] else d[i,j] := (d[i-1,j+1] - d[i-1,j]) / (z[j+i] - z[j]); end; end; // n階差分商 for j := Low(a) to High(a) do a[j] := d[j,0]; // 0.5刻みで 与えられていない値を補間 for i := 0 to 18 do begin d1 := i * 0.5 - 4.5; d2 := f(d1); d3 := hermite(d1, z, a); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3])); end; end.
program Pas0705(arg); {$MODE delphi} uses SysUtils, Math; const // データ点の数 - 1 N = 6; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // Spline (スプライン) 補間 function spline(d:Double; x:array of Double; y:array of Double; z:array of Double):Double; var d1, d2, d3 :Double; i, k :Integer; begin // 補間関数値がどの区間にあるか k := -1; for i := 1 to High(x) do begin if d <= x[i] then begin k := i - 1; break; end; end; if k < 0 then k := i - 1; d1 := x[k+1] - d; d2 := d - x[k]; d3 := x[k+1] - x[k]; result := (z[k] * power(d1,3) + z[k+1] * power(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; end; var i, j :Integer; x :array [0..N] of Double; y :array [0..N] of Double; a :array [0..N] of Double; b :array [0..N] of Double; c :array [0..N] of Double; d :array [0..N] of Double; g :array [0..N] of Double; s :array [0..N] of Double; z :array [0..N] of Double; d1, d2, d3 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := Low(x) to High(x) do begin d1 := i * 1.5 - 4.5; x[i] := d1; y[i] := f(d1); end; // 3項方程式の係数の表を作る for i := 1 to High(x)-1 do begin 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])); end; // 3項方程式を解く (ト−マス法) g[1] := b[1]; s[1] := d[1]; for i := 2 to High(a)-1 do begin g[i] := b[i] - a[i] * c[i-1] / g[i-1]; s[i] := d[i] - a[i] * s[i-1] / g[i-1]; end; z[0] := 0; z[N] := 0; z[N-1] := s[N-1] / g[N-1]; for i := N-2 downto 1 do begin z[i] := (s[i] - c[i] * z[i+1]) / g[i]; end; // 0.5刻みで 与えられていない値を補間 for i := 0 to 18 do begin d1 := i * 0.5 - 4.5; d2 := f(d1); d3 := spline(d1, x, y, z); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3])); end; end.
program Pas0801(arg); {$MODE delphi} uses SysUtils, Math; const // 重力加速度 g = -9.8; // 空気抵抗係数 k = -0.01; // 時間間隔(秒) h = 0.01; // 角度 degree = 45; // 空気抵抗による水平方向の減速分 function fx(vx:Double; vy:Double):Double; begin result := k * Sqrt(vx * vx + vy * vy) * vx; end; // 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Double; vy:Double):Double; begin result := g + (k * Sqrt(vx * vx + vy * vy) * vy); end; var // 角度 radian:Double; // 初速 v:Double; // 水平方向の速度 vx:array [0..1] of Double; // 鉛直方向の速度 vy:array [0..1] of Double; // 経過秒数 t:Double = 0.0; // 位置 x:Double = 0.0; y:Double = 0.0; i:Integer; begin // 角度 radian := degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 v := 250 * 1000 div 3600; // 水平方向の速度 vx[0] := v * Cos(radian); // 鉛直方向の速度 vy[0] := v * Sin(radian); // Euler法 i := 1; while y >= 0.0 do begin // 経過秒数 t := i * h; // 位置 x := x + h * vx[0]; y := y + h * vy[0]; writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%8.5f', [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]; inc(i); end; end.
program Pas0802(arg); {$MODE delphi} uses SysUtils, Math; const // 重力加速度 g = -9.8; // 空気抵抗係数 k = -0.01; // 時間間隔(秒) h = 0.01; // 角度 degree = 45; // 空気抵抗による水平方向の減速分 function fx(vx:Double; vy:Double):Double; begin result := k * Sqrt(vx * vx + vy * vy) * vx; end; // 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Double; vy:Double):Double; begin result := g + (k * Sqrt(vx * vx + vy * vy) * vy); end; var // 角度 radian:Double; // 初速 v:Double; // 水平方向の速度 vx:array [0..2] of Double; // 鉛直方向の速度 vy:array [0..2] of Double; // 経過秒数 t:Double = 0.0; // 位置 x:array [0..2] of Double; y:array [0..2] of Double; i:Integer; begin // 角度 radian := degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 v := 250 * 1000 div 3600; // 水平方向の速度 vx[0] := v * Cos(radian); // 鉛直方向の速度 vy[0] := v * Sin(radian); // 位置 x[0] := 0.0; y[0] := 0.0; // Heun法 i := 1; while y[0] >= 0.0 do begin // 経過秒数 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]; writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%8.5f', [t, vx[0], vy[0], x[0], y[0]])); inc(i); end; end.
program Pas0803(arg); {$MODE delphi} uses SysUtils, Math; const // 重力加速度 g = -9.8; // 空気抵抗係数 k = -0.01; // 時間間隔(秒) h = 0.01; // 角度 degree = 45; // 空気抵抗による水平方向の減速分 function fx(vx:Double; vy:Double):Double; begin result := k * Sqrt(vx * vx + vy * vy) * vx; end; // 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Double; vy:Double):Double; begin result := g + (k * Sqrt(vx * vx + vy * vy) * vy); end; var // 角度 radian:Double; // 初速 v:Double; // 水平方向の速度 vx:array [0..1] of Double; wx:Double; // 鉛直方向の速度 vy:array [0..1] of Double; wy:Double; // 経過秒数 t:Double = 0.0; // 位置 x:array [0..1] of Double; y:array [0..1] of Double; i:Integer; begin // 角度 radian := degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 v := 250 * 1000 div 3600; // 水平方向の速度 vx[0] := v * Cos(radian); // 鉛直方向の速度 vy[0] := v * Sin(radian); // 位置 x[0] := 0.0; y[0] := 0.0; // 中点法 i := 1; while y[0] >= 0.0 do begin // 経過秒数 t := i * h; // 位置・速度 vx[1] := h * fx(vx[0], vy[0]); vy[1] := h * fy(vx[0], vy[0]); wx := vx[0] + vx[1] / 2; 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; writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]])); inc(i); end; end.
program Pas0804(arg); {$MODE delphi} uses SysUtils, Math; const // 重力加速度 g = -9.8; // 空気抵抗係数 k = -0.01; // 時間間隔(秒) h = 0.01; // 角度 degree = 45; // 空気抵抗による水平方向の減速分 function fx(vx:Double; vy:Double):Double; begin result := k * Sqrt(vx * vx + vy * vy) * vx; end; // 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Double; vy:Double):Double; begin result := g + (k * Sqrt(vx * vx + vy * vy) * vy); end; var // 角度 radian:Double; // 初速 v:Double; // 水平方向の速度 vx:array [0..4] of Double; wx:Double; // 鉛直方向の速度 vy:array [0..4] of Double; wy:Double; // 経過秒数 t:Double = 0.0; // 位置 x:array [0..4] of Double; y:array [0..4] of Double; i:Integer; begin // 角度 radian := degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 v := 250 * 1000 div 3600; // 水平方向の速度 vx[0] := v * Cos(radian); // 鉛直方向の速度 vy[0] := v * Sin(radian); // 位置 x[0] := 0.0; y[0] := 0.0; // Runge-Kutta法 i := 1; while y[0] >= 0.0 do begin // 経過秒数 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]); wx := vx[0] + vx[1] / 2; 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[0] + ( x[1] + x[2] * 2 + x[3] * 2 + x[4]) / 6; y[0] := y[0] + ( y[1] + y[2] * 2 + y[3] * 2 + y[4]) / 6; vx[0] := vx[0] + (vx[1] + vx[2] * 2 + vx[3] * 2 + vx[4]) / 6; vy[0] := vy[0] + (vy[1] + vy[2] * 2 + vy[3] * 2 + vy[4]) / 6; writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]])); inc(i); end; end.
program Pas0805(arg); {$MODE delphi} uses SysUtils, Math; const // 重力加速度 g = -9.8; // 空気抵抗係数 k = -0.01; // 時間間隔(秒) h = 0.01; // 角度 degree = 45; // 空気抵抗による水平方向の減速分 function fx(vx:Double; vy:Double):Double; begin result := k * Sqrt(vx * vx + vy * vy) * vx; end; // 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Double; vy:Double):Double; begin result := g + (k * Sqrt(vx * vx + vy * vy) * vy); end; var // 角度 radian:Double; // 初速 v:Double; // 水平方向の速度 vx:array [0..4] of Double; wx:Double; // 鉛直方向の速度 vy:array [0..4] of Double; wy:Double; // 経過秒数 t:Double = 0.0; // 位置 x:array [0..4] of Double; y:array [0..4] of Double; i:Integer; begin // 角度 radian := degree * PI / 180.0; // 初速 250 km/h -> 秒速に変換 v := 250 * 1000 div 3600; // 水平方向の速度 vx[0] := v * Cos(radian); // 鉛直方向の速度 vy[0] := v * Sin(radian); // 位置 x[0] := 0.0; y[0] := 0.0; // Runge-Kutta-Gill法 i := 1; while y[0] >= 0.0 do begin // 経過秒数 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]); wx := vx[0] + vx[1] / 2; 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] * ((Sqrt(2.0) - 1) / 2) + vx[2] * (1 - 1 / Sqrt(2.0)); wy := vy[0] + vy[1] * ((Sqrt(2.0) - 1) / 2) + vy[2] * (1 - 1 / 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] / Sqrt(2.0) + vx[3] * (1 + 1 / Sqrt(2.0)); wy := vy[0] - vy[2] / Sqrt(2.0) + vy[3] * (1 + 1 / 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[0] + ( x[1] + x[2] * (2 - Sqrt(2.0)) + x[3] * (2 + Sqrt(2.0)) + x[4]) / 6; y[0] := y[0] + ( y[1] + y[2] * (2 - Sqrt(2.0)) + y[3] * (2 + Sqrt(2.0)) + y[4]) / 6; vx[0] := vx[0] + (vx[1] + vx[2] * (2 - Sqrt(2.0)) + vx[3] * (2 + Sqrt(2.0)) + vx[4]) / 6; vy[0] := vy[0] + (vy[1] + vy[2] * (2 - Sqrt(2.0)) + vy[3] * (2 + Sqrt(2.0)) + vy[4]) / 6; writeln(format('%4.2f'#9'%8.5f'#9'%9.5f'#9'%9.5f'#9'%9.5f', [t, vx[0], vy[0], x[0], y[0]])); inc(i); end; end.
program Pas0901(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := x * x - 2.0; end; function bisection(a:Double; b:Double):Double; var c: Double; fc: Double; begin while true do begin { 区間 (a, b) の中点 c = (a + b) / 2 } c := (a + b) / 2; writeln(format('%12.10f'#9'%13.10f', [c, c - Sqrt(2)])); fc := f(c); if Abs(fc) < 0.0000000001 then break; if fc < 0 then begin { f(c) < 0 であれば, 解は区間 (c, b) の中に存在 } a := c; end else begin { f(c) > 0 であれば, 解は区間 (a, c) の中に存在 } b := c; end; end; result := c; end; var a: Double = 1.0; b: Double = 2.0; begin writeln(format('%12.10f', [bisection(a, b)])); end.
program Pas0902(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := x * x - 2.0; end; function falseposition(a:Double; b:Double):Double; var c: Double; fc: Double; begin while true do begin { 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 } c := (a * f(b) - b * f(a)) / (f(b) - f(a)); writeln(format('%12.10f'#9'%13.10f', [c, c - Sqrt(2)])); fc := f(c); if Abs(fc) < 0.0000000001 then break; if fc < 0 then begin { f(c) < 0 であれば, 解は区間 (c, b) の中に存在 } a := c; end else begin { f(c) > 0 であれば, 解は区間 (a, c) の中に存在 } b := c; end; end; result := c; end; var a: Double = 1.0; b: Double = 2.0; begin writeln(format('%12.10f', [falseposition(a, b)])); end.
program Pas0903(arg); {$MODE delphi} uses SysUtils, Math; function g(x:Double):Double; begin result := (x / 2) + (1 / x); end; function iterative(x0:Double):Double; var x1: Double; begin while true do begin x1 := g(x0); writeln(format('%12.10f'#9'%13.10f', [x1, x1 - Sqrt(2)])); if Abs(x1 - x0) < 0.0000000001 then break; x0 := x1; end; result := x1; end; var x: Double = 1.0; begin writeln(format('%12.10f', [iterative(x)])); end.
program Pas0904(arg); {$MODE delphi} uses SysUtils, Math; function f0(x:Double):Double; begin result := x * x - 2; end; function f1(x:Double):Double; begin result := 2 * x; end; function newton(x0:Double):Double; var x1: Double; begin while true do begin x1 := x0 - (f0(x0) / f1(x0)); writeln(format('%12.10f'#9'%13.10f', [x1, x1 - Sqrt(2)])); if Abs(x1 - x0) < 0.0000000001 then break; x0 := x1; end; result := x1; end; var x: Double = 2.0; begin writeln(format('%12.10f', [newton(x)])); end.
program Pas0905(arg); {$MODE delphi} uses SysUtils, Math; function f0(x:Double):Double; begin result := x * x - 2; end; function f1(x:Double):Double; begin result := 2 * x; end; function f2(x:Double):Double; begin result := 2; end; function bailey(x0:Double):Double; var x1: Double; begin while true do begin x1 := x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))); writeln(format('%12.10f'#9'%13.10f', [x1, x1 - Sqrt(2)])); if Abs(x1 - x0) < 0.0000000001 then break; x0 := x1; end; result := x1; end; var x: Double = 2.0; begin writeln(format('%12.10f', [bailey(x)])); end.
program Pas0906(arg); {$MODE delphi} uses SysUtils, Math; function f(x:Double):Double; begin result := x * x - 2; end; function secant(x0:Double; x1:Double):Double; var x2: Double; begin while true do begin x2 := x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)); writeln(format('%12.10f'#9'%13.10f', [x1, x1 - Sqrt(2)])); if Abs(x1 - x0) < 0.0000000001 then break; x0 := x1; x1 := x2; end; result := x2; end; var x0: Double = 1.0; x1: Double = 2.0; begin writeln(format('%12.10f', [secant(x0, x1)])); end.
program Pas1001(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // ヤコビの反復法 procedure jacobi(a:TwoDimArray; b:array of Double; var x0:array of Double); var x1:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); finish:Boolean; i, j:Integer; begin while (true) do begin finish := true; for i := Low(x0) to High(x0) do begin x1[i] := 0.0; for j := Low(x0) to High(x0) do if j <> i then x1[i] := x1[i] + a[i,j] * x0[j]; x1[i] := (b[i] - x1[i]) / a[i,i]; if (Abs(x1[i] - x0[i]) > 0.0000000001) then finish := false; end; for i := Low(x0) to High(x0) do x0[i] := x1[i]; if finish then exit; disp_vector(x0); end; end; var a:TwoDimArray = (( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0), (-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0)); b:array [0..N] of Double = (20.0, 16.0, 8.0, 17.0); c:array [0..N] of Double = ( 0.0, 0.0, 0.0, 0.0); begin // ヤコビの反復法 jacobi(a,b,c); writeln('X'); disp_vector(c); end.
program Pas1002(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // ガウス・ザイデル法 procedure gauss(a:TwoDimArray; b:array of Double; var x0:array of Double); var x1 :Double; finish:Boolean; i, j:Integer; begin while (true) do begin finish := true; for i := Low(x0) to High(x0) do begin x1 := 0.0; for j := Low(x0) to High(x0) do if j <> i then x1 := x1 + a[i,j] * x0[j]; x1 := (b[i] - x1) / a[i,i]; if (Abs(x1 - x0[i]) > 0.0000000001) then finish := false; x0[i] := x1; end; if finish then exit; disp_vector(x0); end; end; var a :TwoDimArray = (( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0), (-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0)); b :array [0..N] of Double = (20.0, 16.0, 8.0, 17.0); c :array [0..N] of Double = ( 0.0, 0.0, 0.0, 0.0); begin // ガウス・ザイデル法 gauss(a,b,c); writeln('X'); disp_vector(c); end.
program Pas1003(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var i,j:Integer; begin for i := Low(matrix) to High(matrix) do begin for j := Low(matrix) to High(matrix) do write(format('%14.10f'#9, [matrix[i,j]])); writeln(); end; end; // 前進消去 procedure forward_elimination(var a:TwoDimArray; var b:array of Double); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) - 1 do begin for row := (pivot + 1) to High(a) do begin s := a[row,pivot] / a[pivot,pivot]; for col := pivot to High(a) do a[row,col] := a[row,col] - a[pivot,col] * s; b[row] := b[row] - b[pivot] * s; end; end; end; // 後退代入 procedure backward_substitution(a:TwoDimArray; var b:array of Double); var row, col:Integer; begin for row := High(a) downto Low(a) do begin for col := High(a) downto row + 1 do b[row] := b[row] - a[row,col] * b[col]; b[row] := b[row] / a[row,row]; end; end; // ピボット選択 procedure pivoting(var a:TwoDimArray; var b:array of Double); var max_val, tmp:Double; pivot, row, col, max_row:Integer; begin for pivot := Low(a) to High(a) do begin // 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0; for row := pivot to High(a) do begin if Abs(a[row,pivot]) > max_val then begin // 一番値が大きい行 max_val := Abs(a[row,pivot]); max_row := row; end; end; // 一番値が大きい行と入れ替え if max_row <> pivot then begin tmp := 0; for col := Low(a) to High(a) do begin tmp := a[max_row,col]; a[max_row,col] := a[pivot,col]; a[pivot,col] := tmp; end; tmp := b[max_row]; b[max_row] := b[pivot]; b[pivot] := tmp; end; end; end; var a :TwoDimArray = ((-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0), ( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0)); b :array [0..N] of Double = (8.0, 17.0, 20.0, 16.0); begin // ピボット選択 pivoting(a,b); writeln('pivoting'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 前進消去 forward_elimination(a,b); writeln('forward elimination'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 後退代入 backward_substitution(a,b); writeln('X'); disp_vector(b); end.
program Pas1004(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var i,j:Integer; begin for i := Low(matrix) to High(matrix) do begin for j := Low(matrix) to High(matrix) do write(format('%14.10f'#9, [matrix[i,j]])); writeln(); end; end; // 前進消去 procedure forward_elimination(var a:TwoDimArray; var b:array of Double); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) do begin for row := Low(a) to High(a) do begin if row = pivot then continue; s := a[row,pivot] / a[pivot,pivot]; for col := pivot to High(a) do a[row,col] := a[row,col] - a[pivot,col] * s; b[row] := b[row] - b[pivot] * s; end; end; end; // 後退代入 procedure backward_substitution(a:TwoDimArray; var b:array of Double); var pivot:Integer; begin for pivot := Low(a) to High(a) do b[pivot] := b[pivot] / a[pivot,pivot]; end; // ピボット選択 procedure pivoting(var a:TwoDimArray; var b:array of Double); var max_val, tmp:Double; pivot, row, col, max_row:Integer; begin for pivot := Low(a) to High(a) do begin // 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0; for row := pivot to High(a) do begin if Abs(a[row,pivot]) > max_val then begin // 一番値が大きい行 max_val := Abs(a[row,pivot]); max_row := row; end; end; // 一番値が大きい行と入れ替え if max_row <> pivot then begin tmp := 0; for col := Low(a) to High(a) do begin tmp := a[max_row,col]; a[max_row,col] := a[pivot,col]; a[pivot,col] := tmp; end; tmp := b[max_row]; b[max_row] := b[pivot]; b[pivot] := tmp; end; end; end; var a :TwoDimArray = (( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0), (-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0)); b :array [0..N] of Double = (20.0, 16.0, 8.0, 17.0); begin // ピボット選択 pivoting(a,b); writeln('pivoting'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 前進消去 forward_elimination(a,b); writeln('forward elimination'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 後退代入 backward_substitution(a,b); writeln('X'); disp_vector(b); end.
program Pas1005(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var i,j:Integer; begin for i := Low(matrix) to High(matrix) do begin for j := Low(matrix) to High(matrix) do write(format('%14.10f'#9, [matrix[i,j]])); writeln(); end; end; // 前進消去 procedure forward_elimination(var a:TwoDimArray; var b:array of Double); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) - 1 do begin for row := (pivot + 1) to High(a) do begin s := a[row,pivot] / a[pivot,pivot]; for col := pivot to High(a) do a[row,col] := a[row,col] - a[pivot,col] * s; // これが 上三角行列 a[row,pivot] := s; // これが 下三角行列 // b[row] := b[row] - b[pivot] * s; // この値は変更しない end; end; end; // 前進代入 procedure forward_substitution(a:TwoDimArray; b:array of Double; var y:array of Double); var row, col:Integer; begin for row := Low(a) to High(a) do begin for col := Low(a) to row do b[row] := b[row] - a[row,col] * y[col]; y[row] := b[row]; end; end; // 後退代入 procedure backward_substitution(a:TwoDimArray; y:array of Double; var x:array of Double); var row, col:Integer; begin for row := High(a) downto Low(a) do begin for col := High(a) downto row + 1 do y[row] := y[row] - a[row,col] * x[col]; x[row] := y[row] / a[row,row]; end; end; // ピボット選択 procedure pivoting(var a:TwoDimArray; var b:array of Double); var max_val, tmp:Double; pivot, row, col, max_row:Integer; begin for pivot := Low(a) to High(a) do begin // 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0; for row := pivot to High(a) do begin if Abs(a[row,pivot]) > max_val then begin // 一番値が大きい行 max_val := Abs(a[row,pivot]); max_row := row; end; end; // 一番値が大きい行と入れ替え if max_row <> pivot then begin tmp := 0; for col := Low(a) to High(a) do begin tmp := a[max_row,col]; a[max_row,col] := a[pivot,col]; a[pivot,col] := tmp; end; tmp := b[max_row]; b[max_row] := b[pivot]; b[pivot] := tmp; end; end; end; var a :TwoDimArray = ((-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0), ( 9.0, 2.0, 1.0, 1.0), (2.0, 8.0, -2.0, 1.0)); b :array [0..N] of Double = (8.0, 17.0, 20.0, 16.0); x :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); y :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); begin // ピボット選択 pivoting(a,b); writeln('pivoting'); writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 前進消去 forward_elimination(a,b); writeln('LU'); disp_matrix(a); // Ly=b から y を求める (前進代入) forward_substitution(a,b,y); writeln('Y'); disp_vector(y); // Ux=y から x を求める (後退代入) backward_substitution(a,y,x); writeln('X'); disp_vector(x); end.
program Pas1006(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var i,j:Integer; begin for i := Low(matrix) to High(matrix) do begin for j := Low(matrix) to High(matrix) do write(format('%14.10f'#9, [matrix[i,j]])); writeln(); end; end; // 前進消去 procedure forward_elimination(var a:TwoDimArray; var b:array of Double); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) do begin s := 0; for col := Low(a) to pivot - 1 do s := s + a[pivot,col] * a[pivot,col]; // ここで根号の中が負の値になると計算できない! a[pivot,pivot] := Sqrt(a[pivot,pivot] - s); for row := (pivot + 1) to High(a) do begin s := 0; for col := Low(a) to pivot - 1 do s := s + a[row,col] * a[pivot,col]; a[row,pivot] := (a[row,pivot] - s) / a[pivot,pivot]; a[pivot,row] := a[row,pivot]; end; end; end; // 前進代入 procedure forward_substitution(a:TwoDimArray; b:array of Double; var y:array of Double); var row, col:Integer; begin for row := Low(a) to High(a) do begin for col := Low(a) to row do b[row] := b[row] - a[row,col] * y[col]; y[row] := b[row] / a[row,row]; end; end; // 後退代入 procedure backward_substitution(a:TwoDimArray; y:array of Double; var x:array of Double); var row, col:Integer; begin for row := High(a) downto Low(a) do begin for col := High(a) downto row + 1 do y[row] := y[row] - a[row,col] * x[col]; x[row] := y[row] / a[row,row]; end; end; var a :TwoDimArray = ((5.0,2.0,3.0,4.0),(2.0,10.0,6.0,7.0),(3.0,6.0,15.0,9.0),(4.0,7.0,9.0,20.0)); b :array [0..N] of Double = (34.0,68.0,96.0,125.0); x :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); y :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); begin writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 前進消去 forward_elimination(a,b); writeln('LL^T'); disp_matrix(a); // Ly=b から y を求める (前進代入) forward_substitution(a,b,y); writeln('Y'); disp_vector(y); // L^Tx=y から x を求める (後退代入) backward_substitution(a,y,x); writeln('X'); disp_vector(x); end.
program Pas1007(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var i,j:Integer; begin for i := Low(matrix) to High(matrix) do begin for j := Low(matrix) to High(matrix) do write(format('%14.10f'#9, [matrix[i,j]])); writeln(); end; end; // 前進消去 procedure forward_elimination(var a:TwoDimArray; var b:array of Double); var s:Double; pivot, row, col, k:Integer; begin for pivot := Low(a) to High(a) do begin // pivot < k の場合 s := 0; for col := Low(a) to pivot - 1 do begin s := a[pivot,col]; for k := Low(a) to col - 1 do s := s - a[pivot,k] * a[col,k] * a[k,k]; a[pivot,col] := s / a[col,col]; a[col,pivot] := a[pivot,col]; end; // pivot == k の場合 s := a[pivot,pivot]; for k := Low(a) to pivot - 1 do s := s - a[pivot,k] * a[pivot,k] * a[k,k]; a[pivot,pivot] := s; end; end; // 前進代入 procedure forward_substitution(a:TwoDimArray; b:array of Double; var y:array of Double); var row, col:Integer; begin for row := Low(a) to High(a) do begin for col := Low(a) to row do b[row] := b[row] - a[row,col] * y[col]; y[row] := b[row]; end; end; // 後退代入 procedure backward_substitution(a:TwoDimArray; y:array of Double; var x:array of Double); var row, col:Integer; begin for row := High(a) downto Low(a) do begin for col := High(a) downto row + 1 do y[row] := y[row] - a[row,col] * a[row,row] * x[col]; x[row] := y[row] / a[row,row]; end; end; var a :TwoDimArray = ((5.0,2.0,3.0,4.0),(2.0,10.0,6.0,7.0),(3.0,6.0,15.0,9.0),(4.0,7.0,9.0,20.0)); b :array [0..N] of Double = (34.0,68.0,96.0,125.0); x :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); y :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); begin writeln('A'); disp_matrix(a); writeln('B'); disp_vector(b); writeln(); // 前進消去 forward_elimination(a,b); writeln('LDL^T'); disp_matrix(a); // Ly=b から y を求める (前進代入) forward_substitution(a,b,y); writeln('Y'); disp_vector(y); // DL^Tx=y から x を求める (後退代入) backward_substitution(a,y,x); writeln('X'); disp_vector(x); end.
program Pas1101(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 正規化 (ベクトルの長さを1にする) procedure normarize(var x:array of Double); var s:Double = 0.0; i:Integer; begin for i := Low(x) to High(x) do s := s + x[i] * x[i]; s := Sqrt(s); for i := Low(x) to High(x) do x[i] := x[i] / s; end; // ベキ乗法 function power(a:TwoDimArray; var x0:array of Double):Double; var lambda:Double = 0.0; p0, p1, e0, e1:Double; i, j, k:Integer; x1:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); begin // 正規化 (ベクトル x0 の長さを1にする) normarize(x0); e0 := 0.0; for i := Low(x0) to High(x0) do e0 := e0 + x0[i]; for k := 1 to 100 do begin // 1次元配列を表示 write(format('%3d'#9, [k])); disp_vector(x0); // 行列の積 x1 = A × x0 for i := Low(x1) to High(x1) do x1[i] := 0.0; for i := Low(x1) to High(x1) do begin for j := Low(x0) to High(x0) do x1[i] := x1[i] + a[i][j] * x0[j]; end; // 内積 p0 := 0.0; p1 := 0.0; for i := Low(x0) to High(x0) do begin p0 := p0 + x1[i] * x1[i]; p1 := p1 + x1[i] * x0[i]; end; // 固有値 lambda := p0 / p1; // 正規化 (ベクトル x1 の長さを1にする) normarize(x1); // 収束判定 e1 := 0.0; for i := Low(x0) to High(x0) do e1 := e1 + x1[i]; if Abs(e1 - e0) < 0.00000000001 then break; for i := Low(x0) to High(x0) do x0[i] := x1[i]; e0 := e1; end; power := lambda; end; // ベキ乗法で最大固有値を求める var a :TwoDimArray = ((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)); x :array [0..N] of Double = (1.0, 0.0, 0.0, 0.0); lambda:Double; begin // ベキ乗法 lambda := power(a, x); writeln(); writeln('eigenvalue'); writeln(format('%14.10f', [lambda])); writeln('eigenvector'); disp_vector(x); end.
program Pas1102(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 正規化 (ベクトルの長さを1にする) procedure normarize(var x:array of Double); var s:Double = 0.0; i:Integer; begin for i := Low(x) to High(x) do s := s + x[i] * x[i]; s := Sqrt(s); for i := Low(x) to High(x) do x[i] := x[i] / s; end; // LU分解 procedure forward_elimination(var a:TwoDimArray); var s:Double; pivot, row, col:Integer; begin for pivot := Low(a) to High(a) - 1 do begin for row := (pivot + 1) to High(a) do begin s := a[row,pivot] / a[pivot,pivot]; for col := pivot to High(a) do a[row,col] := a[row,col] - a[pivot,col] * s; // これが 上三角行列 a[row,pivot] := s; // これが 下三角行列 end; end; end; // 前進代入 (Ly = b から y を求める) procedure forward_substitution(a:TwoDimArray; var y:array of Double; b:array of Double); var row, col:Integer; begin for row := Low(a) to High(a) do begin for col := Low(a) to row do b[row] := b[row] - a[row,col] * y[col]; y[row] := b[row]; end; end; // 後退代入 (Ux = y から x を求める) procedure backward_substitution(a:TwoDimArray; var x:array of Double; var y:array of Double); var row, col:Integer; begin for row := High(a) downto Low(a) do begin for col := High(a) downto row + 1 do y[row] := y[row] - a[row,col] * x[col]; x[row] := y[row] / a[row,row]; end; end; // 逆ベキ乗法 function inverse(a:TwoDimArray; var x0:array of Double):Double; var lambda:Double = 0.0; p0, p1, e0, e1:Double; i, k:Integer; b:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); y:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); x1:array [0..N] of Double = (0.0, 0.0, 0.0, 0.0); begin // 正規化 (ベクトル x0 の長さを1にする) normarize(x0); e0 := 0.0; for i := Low(x0) to High(x0) do e0 := e0 + x0[i]; for k := 1 to 100 do begin // 1次元配列を表示 write(format('%3d'#9, [k])); disp_vector(x0); // Ly = b から y を求める (前進代入) for i := Low(x0) to High(x0) do begin b[i] := x0[i]; y[i] := 0; end; forward_substitution(a,y,b); // Ux = y から x を求める (後退代入) for i := Low(x1) to High(x1) do x1[i] := 0; backward_substitution(a,x1,y); // 内積 p0 := 0.0; p1 := 0.0; for i := Low(x0) to High(x0) do begin p0 := p0 + x1[i] * x1[i]; p1 := p1 + x1[i] * x0[i]; end; // 固有値 lambda := p1 / p0; // 正規化 (ベクトル x1 の長さを1にする) normarize(x1); // 収束判定 e1 := 0.0; for i := Low(x0) to High(x0) do e1 := e1 + x1[i]; if Abs(e1 - e0) < 0.00000000001 then break; for i := Low(x0) to High(x0) do x0[i] := x1[i]; e0 := e1; end; inverse := lambda; end; // 逆ベキ乗法で最小固有値を求める var a :TwoDimArray = ((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)); x :array [0..N] of Double = (1.0, 0.0, 0.0, 0.0); lambda:Double; begin // LU分解 forward_elimination(a); // 逆ベキ乗法 lambda := inverse(a, x); writeln(); writeln('eigenvalue'); writeln(format('%14.10f', [lambda])); writeln('eigenvector'); disp_vector(x); end.
program Pas1103(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 対角要素を表示 procedure disp_eigenvalue(a:TwoDimArray); var i:Integer; begin for i := 0 to N do write(format('%14.10f'#9, [a[i][i]])); writeln(); end; // 行列の積 procedure multiply(a:TwoDimArray; b:TwoDimArray; var c:TwoDimArray); var i, j, k:Integer; s:Double; begin for i := 0 to N do begin for j := 0 to N do begin s := 0.0; for k := 0 to N do s := s + (a[i][k] * b[k][j]); c[i][j] := s; end; end; end; // LU分解 procedure decomp(a:TwoDimArray; var l:TwoDimArray; var u:TwoDimArray); var i, j, k:Integer; t:Double; begin for i := 0 to N do begin for j := 0 to N do begin l[i][j] := 0.0; u[i][j] := 0.0; end; end; l[0][0] := 1.0; for j := 0 to N do u[0][j] := a[0][j]; for i := 1 to N do begin u[i][0] := 0.0; l[0][i] := 0.0; l[i][0] := a[i][0] / u[0][0]; end; for i := 1 to N do begin l[i][i] := 1.0; t := a[i][i]; for k := 0 to i do t := t - (l[i][k] * u[k][i]); u[i][i] := t; for j := (i + 1) to N do begin u[j][i] := 0.0; l[i][j] := 0.0; t := a[j][i]; for k := 0 to i do t := t - (l[j][k] * u[k][i]); l[j][i] := t / u[i][i]; t := a[i][j]; for k := 0 to i do t := t - (l[i][k] * u[k][j]); u[i][j] := t; end; end; end; // LR分解で固有値を求める var a :TwoDimArray = ((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)); l :TwoDimArray = ((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)); u :TwoDimArray = ((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)); e:Double; i, j, k:Integer; begin for k := 1 to 200 do begin // LU分解 decomp(a, l, u); // 行列の積 multiply(u, l, a); // 対角要素を表示 write(format('%3d'#9, [k])); disp_eigenvalue(a); // 収束判定 e := 0.0; for i := 1 to N do begin for j := 0 to (i - 1) do e := e + Abs(a[i][j]); end; if e < 0.00000000001 then break; end; writeln(); writeln('eigenvalue'); disp_eigenvalue(a); end.
program Pas1104(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 対角要素を表示 procedure disp_eigenvalue(a:TwoDimArray); var i:Integer; begin for i := 0 to N do write(format('%14.10f'#9, [a[i][i]])); writeln(); end; // 行列の積 procedure multiply(a:TwoDimArray; b:TwoDimArray; var c:TwoDimArray); var i, j, k:Integer; s:Double; begin for i := 0 to N do begin for j := 0 to N do begin s := 0.0; for k := 0 to N do s := s + a[i][k] * b[k][j]; c[i][j] := s; end; end; end; // QR分解 procedure decomp(a:TwoDimArray; var q:TwoDimArray; var r:TwoDimArray); var i, j, k:Integer; t, s:Double; x :array [0..N] of Double; begin for k := 0 to N do begin for i := 0 to N do x[i] := a[i][k]; for j := 0 to (k - 1) do begin t := 0.0; for i := 0 to N do t := t + a[i][k] * q[i][j]; r[j][k] := t; r[k][j] := 0.0; for i := 0 to N do x[i] := x[i] - t * q[i][j]; end; s := 0.0; for i := 0 to N do s := s + x[i] * x[i]; r[k][k] := Sqrt(s); for i := 0 to N do q[i][k] := x[i] / r[k][k]; end; end; // QR分解で固有値を求める var a :TwoDimArray = ((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)); q :TwoDimArray = ((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)); r :TwoDimArray = ((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)); e:Double; i, j, k:Integer; begin for k := 1 to 200 do begin // QR分解 decomp(a, q, r); // 行列の積 multiply(r, q, a); // 対角要素を表示 write(format('%3d'#9, [k])); disp_eigenvalue(a); // 収束判定 e := 0.0; for i := 1 to N do begin for j := 0 to (i - 1) do e := e + Abs(a[i][j]); end; if e < 0.00000000001 then break; end; writeln(); writeln('eigenvalue'); disp_eigenvalue(a); end.
program Pas1105(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 対角要素を表示 procedure disp_eigenvalue(a:TwoDimArray); var i:Integer; begin for i := 0 to N do write(format('%14.10f'#9, [a[i][i]])); writeln(); end; // 正規化 (ベクトルの長さを1にする) procedure normarize(var x:array of Double); var s:Double = 0.0; i:Integer; begin for i := Low(x) to High(x) do s := s + x[i] * x[i]; s := Sqrt(s); for i := Low(x) to High(x) do x[i] := x[i] / s; end; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 固有ベクトルを表示 procedure disp_eigenvector(matrix:TwoDimArray); var i, j:Integer; row:array [0..N] of Double; begin for i := 0 to N do begin for j := 0 to N do row[j] := matrix[i][j]; normarize(row); disp_vector(row); end; end; // ヤコビ法 procedure jacobi(var a:TwoDimArray; var v:TwoDimArray); var i, j, k:Integer; p, q:Integer; max_val:Double; c, s:Double; t, t1, t2:Double; begin for k := 1 to 100 do begin // 最大値を探す max_val := 0.0; for i := 0 to N do begin for j := i + 1 to N do begin if (max_val < Abs(a[i][j])) then begin max_val := Abs(a[i][j]); p := i; q := j; end; end; end; // θ を求める t := 0.0; if Abs(a[p][p] - a[q][q]) < 0.00000000001 then begin // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t := PI / 4.0; if (a[p][p] < 0) then t := -t; end else begin // a_{pp} ≠ a_{qq} のとき t := Arctan(2.0 * a[p][q] / (a[p][p] - a[q][q])) / 2.0; end; // θ を使って 行列 U を作成し、A = U^t × A × U c := Cos(t); s := Sin(t); // U^t × A t1 := 0.0; t2 := 0.0; for i := 0 to N do begin 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; end; // A × U for i := 0 to N do begin 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; end; // 対角要素を表示 write(format('%3d'#9, [k])); disp_eigenvalue(a); // 収束判定 if max_val < 0.00000000001 then break; end; end; // ヤコビ法で固有値を求める var a :TwoDimArray = ((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)); v :TwoDimArray = ((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)); begin // ヤコビ法 jacobi(a, v); writeln(); writeln('eigenvalue'); disp_eigenvalue(a); writeln(); writeln('eigenvector'); disp_eigenvector(v); end.
program Pas1106(arg); {$MODE delphi} uses SysUtils, Math; const N = 3; type TwoDimArray = array [0..N,0..N] of Double; // 1次元配列を表示 procedure disp_vector(row:array of Double); var i:Integer; begin for i := Low(row) to High(row) do write(format('%14.10f'#9, [row[i]])); writeln(); end; // 2次元配列を表示 procedure disp_matrix(matrix:TwoDimArray); var row, col:Integer; begin for row := 0 to N do begin for col := 0 to N do write(format('%14.10f'#9, [matrix[row][col]])); writeln(); end; end; // ハウスホルダー変換 procedure tridiagonalize(var a:TwoDimArray; var d:array of Double; var e:array of Double); var w:array [0..N] of Double; v:array [0..N] of Double; i, j, k:Integer; s, t:Double; begin for k := 0 to N do v[k] := 0.0; for k := 0 to (N - 2) do begin for i := 0 to N do w[i] := 0.0; d[k] := a[k][k]; t := 0.0; for i := (k + 1) to N do begin w[i] := a[i][k]; t := t + w[i] * w[i]; end; s := Sqrt(t); if w[k + 1] < 0 then s := -s; if Abs(s) < 0.00000000001 then begin e[k + 1] := 0.0; end else begin e[k + 1] := -s; w[k + 1] := w[k + 1] + s; s := 1.0 / Sqrt(w[k + 1] * s); for i := (k + 1) to N do w[i] := w[i] * s; for i := (k + 1) to N do begin s := 0.0; for j := (k + 1) to N do if j <= i then s := s + a[i][j] * w[j] else s := s + a[j][i] * w[j]; v[i] := s; end; s := 0.0; for i := (k + 1) to N do s := s + w[i] * v[i]; s := s / 2.0; for i := (k + 1) to N do v[i] := v[i] - s * w[i]; for i := (k + 1) to N do for j := (k + 1) to i do a[i][j] := a[i][j] - (w[i] * v[j] + w[j] * v[i]); // 次の行は固有ベクトルを求めないなら不要 for i := (k + 1) to N do a[i][k] := w[i]; end; end; d[N - 1] := a[N - 1][N - 1]; d[N] := a[N][N]; e[0] := 0.0; e[N] := a[N][N - 1]; // 次の行は固有ベクトルを求めないなら不要 for k := N downto 0 do begin for i := 0 to N do w[i] := 0.0; if k < N - 1 then begin for i := (k + 1) to N do w[i] := a[i][k]; for i := (k + 1) to N do begin s := 0.0; for j := (k + 1) to N do s := s + a[i][j] * w[j]; v[i] := s; end; for i := (k + 1) to N do for j := (k + 1) to N do a[i][j] := a[i][j] - v[i] * w[j]; end; for i := 0 to N do a[i][k] := 0.0; a[k][k] := 1.0; end; end; // QR分解 procedure decomp(var a:TwoDimArray; var d:array of Double; var e:array of Double); var h, l:Integer; i, j, k:Integer; w, s :Double; x, y, z :Double; t, u :Double; begin e[0] := 1.0; h := N; while (Abs(e[h]) < 0.00000000001) do dec(h); while (h > 0) do begin e[0] := 0.0; l := h - 1; while (Abs(e[l]) >= 0.00000000001) do dec(l); for j := 1 to 100 do begin w := (d[h - 1] - d[h]) / 2.0; s := Sqrt(w * w + e[h] * e[h]); if w < 0.0 then s := -s; x := d[l] - d[h] + e[h] * e[h] / (w + s); y := e[l + 1]; z := 0.0; t := 0.0; u := 0.0; for k := l to (h - 1) do begin if Abs(x) >= Abs(y) then begin t := -y / x; u := 1.0 / Sqrt(t * t + 1.0); s := t * u; end else begin t := -x / y; s := 1.0 / Sqrt(t * t + 1.0); if t < 0 then s := -s; u := t * s; end; w := d[k] - d[k + 1]; t := (w * s + 2.0 * 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 i := 0 to N do begin x := a[k , i]; y := a[k + 1, i]; a[k , i] := u * x - s * y; a[k + 1, i] := s * x + u * y; end; if k < N - 1 then begin x := e[k + 1]; y := -s * e[k + 2]; z := y; e[k + 2] := u * e[k + 2]; end; end; write(format('%3d'#9, [j])); disp_vector(d); // 収束判定 if (Abs(e[h]) < 0.00000000001) then break; end; e[0] := 1.0; while (Abs(e[h]) < 0.00000000001) do dec(h); end; // 次の行は固有ベクトルを求めないなら不要 for k := 0 to (N - 1) do begin l := k; for i := (k + 1) to N do if d[i] > d[l] then l := i; t := d[k]; d[k] := d[l]; d[l] := t; for i := 0 to N do begin t := a[k][i]; a[k][i] := a[l][i]; a[l][i] := t; end; end; end; // ハウスホルダー変換とQR分解で固有値を求める var a :TwoDimArray = ((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)); d:array [0..N] of Double; e:array [0..N] of Double; begin // ハウスホルダー変換 tridiagonalize(a, d, e); // QR分解 decomp(a, d, e); writeln(); writeln('eigenvalue'); disp_vector(d); writeln(); writeln('eigenvector'); disp_matrix(a); end.