with TEXT_IO, Ada.Integer_Text_IO, Ada.Float_Text_IO; use TEXT_IO, Ada.Integer_Text_IO, Ada.Float_Text_IO; procedure Ada0101 is begin Put(3 + 5); New_Line; Put(3 + 5, Width=> 0); New_Line; Put(3 * 5, Width=> 0); New_Line; Put(3 ** 5, Width=> 0); New_Line; Put(5 / 3, Width=> 0); New_Line; Put(5 mod 3, Width=> 0); New_Line; Put(5 rem 3, Width=> 0); New_Line; Put(Float(5) / Float(3) ); New_Line; Put(Float(5) / Float(3), Exp=>0); New_Line; Put(Float(5) / Float(3), Aft=>10, Exp=>0); New_Line; Put(Float(5) / Float(3), Fore=>5, Aft=>10, Exp=>0); New_Line; end Ada0101;
with TEXT_IO, Ada.Integer_Text_IO; use TEXT_IO, Ada.Integer_Text_IO; procedure Ada0102 is i:Integer := 3 * 5; begin Put("3 + 5 = "); Put(i, Width=> 0); New_Line; end Ada0102;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; procedure Ada0601 is a : Constant Long_Float := 0.0; b : Constant Long_Float := 1.0; h, s, x : Long_Float; n : Integer; function f(x:Long_Float) return Long_Float is begin return 4.0 / (1.0 + x * x); end f; begin -- 台形則で積分 n := 2; for j in 1..10 loop h := (b - a) / Long_Float(n); s := 0.0; x := a; for i in 1..(n - 1) loop x := x + h; s := s + f(x); end loop; s := h * ((f(a) + f(b)) / 2.0 + s); n := n * 2; -- 結果を π と比較 Put(j, Width=> 2); Put(" : "); Put(s , Fore=>3, Aft=>10, Exp=>0); Put(", "); Put(s - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0); New_Line; end loop; end Ada0601;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; procedure Ada0602 is a : Constant Long_Float := 0.0; b : Constant Long_Float := 1.0; h, s, x : Long_Float; n : Integer; function f(x:Long_Float) return Long_Float is begin return 4.0 / (1.0 + x * x); end f; begin -- 中点則で積分 n := 2; for j in 1..10 loop h := (b - a) / Long_Float(n); s := 0.0; x := a + (h / 2.0); for i in 1..n loop s := s + f(x); x := x + h; end loop; s := h * s; n := n * 2; -- 結果を π と比較 Put(j, Width=> 2); Put(" : "); Put(s , Fore=>3, Aft=>10, Exp=>0); Put(", "); Put(s - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0); New_Line; end loop; end Ada0602;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; procedure Ada0603 is a : Constant Long_Float := 0.0; b : Constant Long_Float := 1.0; h, s, x : Long_Float; s2, s4 : Long_Float; n : Integer; function f(x:Long_Float) return Long_Float is begin return 4.0 / (1.0 + x * x); end f; begin -- Simpson則で積分 n := 2; for j in 1..5 loop h := (b - a) / Long_Float(n); s2 := 0.0; s4 := 0.0; x := a + h; for i in 1..(n / 2) loop s4 := s4 + f(x); x := x + h; s2 := s2 + f(x); x := x + h; end loop; s2 := (s2 - f(b)) * 2.0 + f(a) + f(b); s4 := s4 * 4.0; s := (s2 + s4) * h / 3.0; n := n * 2; -- 結果を π と比較 Put(j, Width=> 2); Put(" : "); Put(s , Fore=>3, Aft=>10, Exp=>0); Put(", "); Put(s - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0); New_Line; end loop; end Ada0603;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Integer_Text_IO, Ada.Numerics; procedure Ada0604 is a : Constant Long_Float := 0.0; b : Constant Long_Float := 1.0; h, s, x : Long_Float; n : Integer; t : array (1..6, 1..6) of Long_Float; function f(x:Long_Float) return Long_Float is begin return 4.0 / (1.0 + x * x); end f; begin -- 台形則で積分 n := 2; for i in 1..6 loop h := (b - a) / Long_Float(n); s := 0.0; x := a; for j in 1..(n - 1) loop x := x + h; s := s + f(x); end loop; -- 結果を保存 t(i,1) := h * ((f(a) + f(b)) / 2.0 + s); n := n * 2; end loop; -- Richardsonの補外法 n := 4; for j in 2..6 loop for i in j..6 loop t(i,j) := t(i, j - 1) + (t(i, j - 1) - t(i - 1, j - 1)) / Long_Float(n - 1); if i = j then -- 結果を π と比較 Put(j, Width=> 2); Put(" : "); Put(t(i,j) , Fore=>3, Aft=>10, Exp=>0); Put(", "); Put(t(i,j) - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0); New_Line; end if; end loop; n := n * 4; end loop; end Ada0604;
with TEXT_IO, Ada.Long_Float_Text_IO; use TEXT_IO, Ada.Long_Float_Text_IO; procedure Ada0701 is -- データ点の数 - 1 N : Constant Integer := 6; type Long_Float_Array is array (0..N) of Long_Float; x : Long_Float_Array; y : Long_Float_Array; d, d1, d2 : Long_Float; -- 元の関数 function f(x:Long_Float) return Long_Float is begin return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2); end f; -- Lagrange (ラグランジュ) 補間 function lagrange(d:Long_Float; x:Long_Float_Array; y:Long_Float_Array) return Long_Float is sum, prod :Long_Float; begin sum := Long_Float(0.0); for i in x'Range loop prod := y(i); for j in x'Range loop if j /= i then prod := prod * (d - x(j)) / (x(i) - x(j)); end if; end loop; sum := sum + prod; end loop; return sum; end; begin -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i in x'Range loop d := Long_Float(i) * 1.5 - 4.5; x(i) := d; y(i) := f(d); end loop; -- 0.5刻みで 与えられていない値を補間 for i in 0..18 loop d := Long_Float(i) * 0.5 - 4.5; d1 := f(d); d2 := lagrange(d, x, y); -- 元の関数と比較 Put(d, Fore=>2, Aft=>2, Exp=>0); Put(Ascii.HT); Put(d1, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d2, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d1 - d2, Fore=>3, Aft=>5, Exp=>0); New_Line; end loop; end Ada0701;
with TEXT_IO, Ada.Long_Float_Text_IO; use TEXT_IO, Ada.Long_Float_Text_IO; procedure Ada0702 is -- データ点の数 - 1 N : Constant Integer := 6; type Long_Float_Array is array (0..N) of Long_Float; x : Long_Float_Array; y : Long_Float_Array; d, d1, d2 : Long_Float; -- 元の関数 function f(x:Long_Float) return Long_Float is begin return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2); end f; -- Neville (ネヴィル) 補間 function neville(d:Long_Float; x:Long_Float_Array; y:Long_Float_Array) return Long_Float is w : array (0..N, 0..N) of Long_Float; begin for i in y'Range loop w(0, i) := y(i); end loop; for j in 1 .. x'Last loop for i in x'First .. x'Last-j loop 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 loop; end loop; return w(N,0); end; begin -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i in x'Range loop d := Long_Float(i) * 1.5 - 4.5; x(i) := d; y(i) := f(d); end loop; -- 0.5刻みで 与えられていない値を補間 for i in 0..18 loop d := Long_Float(i) * 0.5 - 4.5; d1 := f(d); d2 := neville(d, x, y); -- 元の関数と比較 Put(d, Fore=>2, Aft=>2, Exp=>0); Put(Ascii.HT); Put(d1, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d2, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d1 - d2, Fore=>3, Aft=>5, Exp=>0); New_Line; end loop; end Ada0702;
with TEXT_IO, Ada.Long_Float_Text_IO; use TEXT_IO, Ada.Long_Float_Text_IO; procedure Ada0703 is -- データ点の数 - 1 N : Constant Integer := 6; type Long_Float_Array is array (0..N) of Long_Float; x : Long_Float_Array; y : Long_Float_Array; a : Long_Float_Array; d : array (0..N, 0..N) of Long_Float; d1, d2, d3 : Long_Float; -- 元の関数 function f(x:Long_Float) return Long_Float is begin return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2); end f; -- Newton (ニュートン) 補間 function newton(d:Long_Float; x:Long_Float_Array; a:Long_Float_Array) return Long_Float is sum, prod :Long_Float; begin sum := a(0); for i in 1 .. x'Last loop prod := a(i); for j in x'First .. i-1 loop prod := prod * (d - x(j)); end loop; sum := sum + prod; end loop; return sum; end; begin -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i in x'Range loop d1 := Long_Float(i) * 1.5 - 4.5; x(i) := d1; y(i) := f(d1); end loop; -- 差分商の表を作る for j in y'Range loop d(0,j) := y(j); end loop; for i in 1 .. x'Last loop for j in x'First .. x'Last-i loop d(i,j) := (d(i-1,j+1) - d(i-1,j)) / (x(j+i) - x(j)); end loop; end loop; -- n階差分商 for j in a'Range loop a(j) := d(j,0); end loop; -- 0.5刻みで 与えられていない値を補間 for i in 0..18 loop d1 := Long_Float(i) * 0.5 - 4.5; d2 := f(d1); d3 := newton(d1, x, a); -- 元の関数と比較 Put(d1, Fore=>2, Aft=>2, Exp=>0); Put(Ascii.HT); Put(d2, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d3, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d2 - d3, Fore=>3, Aft=>5, Exp=>0); New_Line; end loop; end Ada0703;
with TEXT_IO, Ada.Long_Float_Text_IO; use TEXT_IO, Ada.Long_Float_Text_IO; procedure Ada0704 is -- データ点の数 - 1 N : Constant Integer := 6; Nx2 : Constant Integer := 13; type Long_Float_ArrayN is array (0..N) of Long_Float; type Long_Float_ArrayNx2 is array (0..Nx2) of Long_Float; x : Long_Float_ArrayN; y : Long_Float_ArrayN; yd : Long_Float_ArrayN; z : Long_Float_ArrayNx2; a : Long_Float_ArrayNx2; d : array (0..Nx2, 0..Nx2) of Long_Float; d1, d2, d3 : Long_Float; k : Integer; -- 元の関数 function f(x:Long_Float) return Long_Float is begin return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2); end f; -- 導関数 function fd(x:Long_Float) return Long_Float is begin return 1.0 - Long_Float(x ** 2) / Long_Float(2) + Long_Float(x ** 4) / Long_Float(4 * 3 * 2); end fd; -- Hermite (エルミート) 補間 function hermite(d:Long_Float; z:Long_Float_ArrayNx2; a:Long_Float_ArrayNx2) return Long_Float is sum, prod :Long_Float; begin sum := a(0); for i in 1 .. a'Last loop prod := a(i); for j in z'First .. i-1 loop prod := prod * (d - z(j)); end loop; sum := sum + prod; end loop; return sum; end; begin -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i in x'Range loop d1 := Long_Float(i) * 1.5 - 4.5; x(i) := d1; y(i) := f(d1); yd(i) := fd(d1); end loop; -- 差分商の表を作る for i in z'Range loop k := i / 2; z(i) := x(k); d(0,i) := y(k); end loop; for i in 1 .. z'Last loop for j in z'First .. z'Last-i loop if (i = 1) and (j mod 2 = 0) then d(i,j) := yd(j / 2); else d(i,j) := (d(i-1,j+1) - d(i-1,j)) / (z(j+i) - z(j)); end if; end loop; end loop; -- n階差分商 for j in a'Range loop a(j) := d(j,0); end loop; -- 0.5刻みで 与えられていない値を補間 for i in 0..18 loop d1 := Long_Float(i) * 0.5 - 4.5; d2 := f(d1); d3 := hermite(d1, z, a); -- 元の関数と比較 Put(d1, Fore=>2, Aft=>2, Exp=>0); Put(Ascii.HT); Put(d2, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d3, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d2 - d3, Fore=>3, Aft=>5, Exp=>0); New_Line; end loop; end Ada0704;
with TEXT_IO, Ada.Long_Float_Text_IO; use TEXT_IO, Ada.Long_Float_Text_IO; procedure Ada0705 is -- データ点の数 - 1 N : Constant Integer := 6; type Long_Float_Array is array (0..N) of Long_Float; x : Long_Float_Array; y : Long_Float_Array; a : Long_Float_Array; b : Long_Float_Array; c : Long_Float_Array; d : Long_Float_Array; g : Long_Float_Array; s : Long_Float_Array; z : Long_Float_Array; d1, d2, d3 : Long_Float; k : Integer; -- 元の関数 function f(x:Long_Float) return Long_Float is begin return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2); end f; -- Spline (スプライン) 補間 function spline(d:Long_Float; x:Long_Float_Array; y:Long_Float_Array; z:Long_Float_Array) return Long_Float is d1, d2, d3 : Long_Float; k : Integer; begin -- 補間関数値がどの区間にあるか k := -1; for i in 1 .. x'Last loop if d <= x(i) then k := i - 1; exit; end if; end loop; if k < 0 then k := N; end if; d1 := x(k+1) - d; d2 := d - x(k); d3 := x(k+1) - x(k); return (z(k) * (d1 ** 3) + z(k+1) * (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; begin -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i in x'Range loop d1 := Long_Float(i) * 1.5 - 4.5; x(i) := d1; y(i) := f(d1); end loop; -- 3項方程式の係数の表を作る for i in 1..x'Last-1 loop 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 loop; -- 3項方程式を解く (ト−マス法) g(1) := b(1); s(1) := d(1); for i in 2..a'Last-1 loop 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 loop; z(0) := 0.0; z(N) := 0.0; z(N-1) := s(N-1) / g(N-1); k := N - 2; while k > 0 loop z(k) := (s(k) - c(k) * z(k+1)) / g(k); k := k - 1; end loop; -- 0.5刻みで 与えられていない値を補間 for i in 0..18 loop d1 := Long_Float(i) * 0.5 - 4.5; d2 := f(d1); d3 := spline(d1, x, y, z); -- 元の関数と比較 Put(d1, Fore=>2, Aft=>2, Exp=>0); Put(Ascii.HT); Put(d2, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d3, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d2 - d3, Fore=>3, Aft=>5, Exp=>0); New_Line; end loop; end Ada0705;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; procedure Ada0801 is -- 重力加速度 g : Constant Long_Float := -9.8; -- 空気抵抗係数 k : Constant Long_Float := -0.01; -- 時間間隔(秒) h : Constant Long_Float := 0.01; -- 角度 degree : Constant Long_Float := 45.0; -- 空気抵抗による水平方向の減速分 function fx(vx:Long_Float; vy:Long_Float) return Long_Float is begin return k * Sqrt(vx * vx + vy * vy) * vx; end fx; -- 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Long_Float; vy:Long_Float) return Long_Float is begin return g + (k * Sqrt(vx * vx + vy * vy) * vy); end fy; -- 角度 radian:Long_Float; -- 初速 v:Long_Float; -- 水平方向の速度 vx:array (0..1) of Long_Float; -- 鉛直方向の速度 vy:array (0..1) of Long_Float; -- 経過秒数 t:Long_Float := 0.0; -- 位置 x:Long_Float := 0.0; y:Long_Float := 0.0; i:Integer; begin -- 角度 radian := degree * Pi / 180.0; -- 初速 250 km/h -> 秒速に変換 v := Long_Float(250 * 1000 / 3600); -- 水平方向の速度 vx(0) := v * Cos(radian); -- 鉛直方向の速度 vy(0) := v * Sin(radian); -- Euler法 i := 1; while y >= 0.0 loop -- 経過秒数 t := Long_Float(i) * h; -- 位置 x := x + h * vx(0); y := y + h * vy(0); Put(t, Fore=>1, Aft=>2, Exp=>0); Put(Ascii.HT); Put(vx(1), Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(vy(1), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(x, Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(y, Fore=>3, Aft=>5, Exp=>0); New_Line; -- 速度 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); i := i + 1; end loop; end Ada0801;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; procedure Ada0802 is -- 重力加速度 g : Constant Long_Float := -9.8; -- 空気抵抗係数 k : Constant Long_Float := -0.01; -- 時間間隔(秒) h : Constant Long_Float := 0.01; -- 角度 degree : Constant Long_Float := 45.0; -- 空気抵抗による水平方向の減速分 function fx(vx:Long_Float; vy:Long_Float) return Long_Float is begin return k * Sqrt(vx * vx + vy * vy) * vx; end fx; -- 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Long_Float; vy:Long_Float) return Long_Float is begin return g + (k * Sqrt(vx * vx + vy * vy) * vy); end fy; -- 角度 radian:Long_Float; -- 初速 v:Long_Float; -- 水平方向の速度 vx:array (0..2) of Long_Float; -- 鉛直方向の速度 vy:array (0..2) of Long_Float; -- 経過秒数 t:Long_Float := 0.0; -- 位置 x:array (0..2) of Long_Float; y:array (0..2) of Long_Float; i:Integer; begin -- 角度 radian := degree * Pi / 180.0; -- 初速 250 km/h -> 秒速に変換 v := Long_Float(250 * 1000 / 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 loop -- 経過秒数 t := Long_Float(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.0; y(2) := y(0) + h * ( vy(0) + vy(1) ) / 2.0; vx(2) := vx(0) + h * (fx(vx(0), vy(0)) + fx(vx(1), vy(1))) / 2.0; vy(2) := vy(0) + h * (fy(vx(0), vy(0)) + fy(vx(1), vy(1))) / 2.0; x(0) := x(2); y(0) := y(2); vx(0) := vx(2); vy(0) := vy(2); Put(t, Fore=>1, Aft=>2, Exp=>0); Put(Ascii.HT); Put(vx(0), Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(vy(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(x(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(y(0), Fore=>4, Aft=>5, Exp=>0); New_Line; i := i + 1; end loop; end Ada0802;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; procedure Ada0803 is -- 重力加速度 g : Constant Long_Float := -9.8; -- 空気抵抗係数 k : Constant Long_Float := -0.01; -- 時間間隔(秒) h : Constant Long_Float := 0.01; -- 角度 degree : Constant Long_Float := 45.0; -- 空気抵抗による水平方向の減速分 function fx(vx:Long_Float; vy:Long_Float) return Long_Float is begin return k * Sqrt(vx * vx + vy * vy) * vx; end fx; -- 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Long_Float; vy:Long_Float) return Long_Float is begin return g + (k * Sqrt(vx * vx + vy * vy) * vy); end fy; -- 角度 radian:Long_Float; -- 初速 v:Long_Float; -- 水平方向の速度 vx:array (0..1) of Long_Float; wx:Long_Float; -- 鉛直方向の速度 vy:array (0..1) of Long_Float; wy:Long_Float; -- 経過秒数 t:Long_Float := 0.0; -- 位置 x:array (0..1) of Long_Float; y:array (0..1) of Long_Float; i:Integer; begin -- 角度 radian := degree * Pi / 180.0; -- 初速 250 km/h -> 秒速に変換 v := Long_Float(250 * 1000 / 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 loop -- 経過秒数 t := Long_Float(i) * h; -- 位置・速度 vx(1) := h * fx(vx(0), vy(0)); vy(1) := h * fy(vx(0), vy(0)); wx := vx(0) + vx(1) / 2.0; wy := vy(0) + vy(1) / 2.0; vx(0) := vx(0) + h * fx(wx, wy); vy(0) := vy(0) + h * fy(wx, wy); x(0) := x(0) + h * wx; y(0) := y(0) + h * wy; Put(t, Fore=>1, Aft=>2, Exp=>0); Put(Ascii.HT); Put(vx(0), Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(vy(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(x(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(y(0), Fore=>4, Aft=>5, Exp=>0); New_Line; i := i + 1; end loop; end Ada0803;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; procedure Ada0804 is -- 重力加速度 g : Constant Long_Float := -9.8; -- 空気抵抗係数 k : Constant Long_Float := -0.01; -- 時間間隔(秒) h : Constant Long_Float := 0.01; -- 角度 degree : Constant Long_Float := 45.0; -- 空気抵抗による水平方向の減速分 function fx(vx:Long_Float; vy:Long_Float) return Long_Float is begin return k * Sqrt(vx * vx + vy * vy) * vx; end fx; -- 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Long_Float; vy:Long_Float) return Long_Float is begin return g + (k * Sqrt(vx * vx + vy * vy) * vy); end fy; -- 角度 radian:Long_Float; -- 初速 v:Long_Float; -- 水平方向の速度 vx:array (0..4) of Long_Float; wx:Long_Float; -- 鉛直方向の速度 vy:array (0..4) of Long_Float; wy:Long_Float; -- 経過秒数 t:Long_Float := 0.0; -- 位置 x:array (0..4) of Long_Float; y:array (0..4) of Long_Float; i:Integer; begin -- 角度 radian := degree * Pi / 180.0; -- 初速 250 km/h -> 秒速に変換 v := Long_Float(250 * 1000 / 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 loop -- 経過秒数 t := Long_Float(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.0; wy := vy(0) + vy(1) / 2.0; 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.0; wy := vy(0) + vy(2) / 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(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.0 + x(3) * 2.0 + x(4)) / 6.0; y(0) := y(0) + ( y(1) + y(2) * 2.0 + y(3) * 2.0 + y(4)) / 6.0; vx(0) := vx(0) + (vx(1) + vx(2) * 2.0 + vx(3) * 2.0 + vx(4)) / 6.0; vy(0) := vy(0) + (vy(1) + vy(2) * 2.0 + vy(3) * 2.0 + vy(4)) / 6.0; Put(t, Fore=>1, Aft=>2, Exp=>0); Put(Ascii.HT); Put(vx(0), Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(vy(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(x(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(y(0), Fore=>4, Aft=>5, Exp=>0); New_Line; i := i + 1; end loop; end Ada0804;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions; procedure Ada0805 is -- 重力加速度 g : Constant Long_Float := -9.8; -- 空気抵抗係数 k : Constant Long_Float := -0.01; -- 時間間隔(秒) h : Constant Long_Float := 0.01; -- 角度 degree : Constant Long_Float := 45.0; -- 空気抵抗による水平方向の減速分 function fx(vx:Long_Float; vy:Long_Float) return Long_Float is begin return k * Sqrt(vx * vx + vy * vy) * vx; end fx; -- 重力と空気抵抗による鉛直方向の減速分 function fy(vx:Long_Float; vy:Long_Float) return Long_Float is begin return g + (k * Sqrt(vx * vx + vy * vy) * vy); end fy; -- 角度 radian:Long_Float; -- 初速 v:Long_Float; -- 水平方向の速度 vx:array (0..4) of Long_Float; wx:Long_Float; -- 鉛直方向の速度 vy:array (0..4) of Long_Float; wy:Long_Float; -- 経過秒数 t:Long_Float := 0.0; -- 位置 x:array (0..4) of Long_Float; y:array (0..4) of Long_Float; i:Integer; begin -- 角度 radian := degree * Pi / 180.0; -- 初速 250 km/h -> 秒速に変換 v := Long_Float(250 * 1000 / 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 loop -- 経過秒数 t := Long_Float(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.0; wy := vy(0) + vy(1) / 2.0; 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.0) / 2.0) + vx(2) * (1.0 - 1.0 / Sqrt(2.0)); wy := vy(0) + vy(1) * ((Sqrt(2.0) - 1.0) / 2.0) + vy(2) * (1.0 - 1.0 / 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.0 + 1.0 / Sqrt(2.0)); wy := vy(0) - vy(2) / Sqrt(2.0) + vy(3) * (1.0 + 1.0 / 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.0 - Sqrt(2.0)) + x(3) * (2.0 + Sqrt(2.0)) + x(4)) / 6.0; y(0) := y(0) + ( y(1) + y(2) * (2.0 - Sqrt(2.0)) + y(3) * (2.0 + Sqrt(2.0)) + y(4)) / 6.0; vx(0) := vx(0) + (vx(1) + vx(2) * (2.0 - Sqrt(2.0)) + vx(3) * (2.0 + Sqrt(2.0)) + vx(4)) / 6.0; vy(0) := vy(0) + (vy(1) + vy(2) * (2.0 - Sqrt(2.0)) + vy(3) * (2.0 + Sqrt(2.0)) + vy(4)) / 6.0; Put(t, Fore=>1, Aft=>2, Exp=>0); Put(Ascii.HT); Put(vx(0), Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(vy(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(x(0), Fore=>4, Aft=>5, Exp=>0); Put(Ascii.HT); Put(y(0), Fore=>4, Aft=>5, Exp=>0); New_Line; i := i + 1; end loop; end Ada0805;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada0901 is function f(x:Long_Float) return Long_Float is begin return x * x - 2.0; end f; function bisection(x0:Long_Float; x1:Long_Float) return Long_Float is a: Long_Float; b: Long_Float; c: Long_Float; fc: Long_Float; begin a := x0; b := x1; while true loop -- 区間 (a, b) の中点 c = (a + b) / 2 c := (a + b) / 2.0; Put(c, Fore=>2, Aft=>10, Exp=>0); Put(Ascii.HT); Put(c - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0); New_Line; fc := f(c); if Abs(fc) < 0.0000000001 then exit; end if; if fc < 0.0 then -- f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a := c; else -- f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b := c; end if; end loop; return c; end bisection; a: Long_Float := 1.0; b: Long_Float := 2.0; begin Put(bisection(a, b), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0901;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada0902 is function f(x:Long_Float) return Long_Float is begin return x * x - 2.0; end f; function falseposition(x0:Long_Float; x1:Long_Float) return Long_Float is a: Long_Float; b: Long_Float; c: Long_Float; fc: Long_Float; begin a := x0; b := x1; while true loop -- 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c := (a * f(b) - b * f(a)) / (f(b) - f(a)); Put(c, Fore=>2, Aft=>10, Exp=>0); Put(Ascii.HT); Put(c - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0); New_Line; fc := f(c); if Abs(fc) < 0.0000000001 then exit; end if; if fc < 0.0 then -- f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a := c; else -- f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b := c; end if; end loop; return c; end falseposition; a: Long_Float := 1.0; b: Long_Float := 2.0; begin Put(falseposition(a, b), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0902;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada0903 is function g(x:Long_Float) return Long_Float is begin return (x / 2.0) + (1.0 / x); end g; function iterative(x:Long_Float) return Long_Float is x0: Long_Float; x1: Long_Float; begin x0 := x; while true loop x1 := g(x0); Put(x1, Fore=>2, Aft=>10, Exp=>0); Put(Ascii.HT); Put(x1 - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0); New_Line; if Abs(x1 - x0) < 0.0000000001 then exit; end if; x0 := x1; end loop; return x1; end iterative; x: Long_Float := 1.0; begin Put(iterative(x), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0903;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada0904 is function f0(x:Long_Float) return Long_Float is begin return x * x - 2.0; end f0; function f1(x:Long_Float) return Long_Float is begin return 2.0 * x; end f1; function newton(x:Long_Float) return Long_Float is x0: Long_Float; x1: Long_Float; begin x0 := x; while true loop x1 := x0 - (f0(x0) / f1(x0)); Put(x1, Fore=>2, Aft=>10, Exp=>0); Put(Ascii.HT); Put(x1 - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0); New_Line; if Abs(x1 - x0) < 0.0000000001 then exit; end if; x0 := x1; end loop; return x1; end newton; x: Long_Float := 2.0; begin Put(newton(x), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0904;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada0905 is function f0(x:Long_Float) return Long_Float is begin return x * x - 2.0; end f0; function f1(x:Long_Float) return Long_Float is begin return 2.0 * x; end f1; function f2(x:Long_Float) return Long_Float is begin return 2.0; end f2; function bailey(x:Long_Float) return Long_Float is x0: Long_Float; x1: Long_Float; begin x0 := x; while true loop x1 := x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2.0 * f1(x0))))); Put(x1, Fore=>2, Aft=>10, Exp=>0); Put(Ascii.HT); Put(x1 - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0); New_Line; if Abs(x1 - x0) < 0.0000000001 then exit; end if; x0 := x1; end loop; return x1; end bailey; x: Long_Float := 2.0; begin Put(bailey(x), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0905;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada0906 is function f(x:Long_Float) return Long_Float is begin return x * x - 2.0; end f; function secant(a:Long_Float; b:Long_Float) return Long_Float is x0: Long_Float; x1: Long_Float; x2: Long_Float; begin x0 := a; x1 := b; while true loop x2 := x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)); Put(x1, Fore=>2, Aft=>10, Exp=>0); Put(Ascii.HT); Put(x1 - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0); New_Line; if Abs(x1 - x0) < 0.0000000001 then exit; end if; x0 := x1; x1 := x2; end loop; return x2; end secant; x0: Long_Float := 1.0; x1: Long_Float := 2.0; begin Put(secant(x0, x1), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0906;
with TEXT_IO, Ada.Long_Float_Text_IO; use TEXT_IO, Ada.Long_Float_Text_IO; procedure Ada1001 is N:Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (20.0, 16.0, 8.0, 17.0); c:Long_Float_Array := ( 0.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- ヤコビの反復法 procedure jacobi(a:Long_Float_TwoDimArray; b:Long_Float_Array; x0:in out Long_Float_Array) is x1:Long_Float_Array := ( 0.0, 0.0, 0.0, 0.0); finish:Boolean; begin while true loop finish := true; for i in x0'Range loop x1(i) := 0.0; for j in x0'Range loop if j /= i then x1(i) := x1(i) + a(i,j) * x0(j); end if; end loop; x1(i) := (b(i) - x1(i)) / a(i,i); if (Abs(x1(i) - x0(i)) > 0.0000000001) then finish := false; end if; end loop; for i in x0'Range loop x0(i) := x1(i); end loop; if finish then exit; end if; disp_vector(x0); end loop; end jacobi; begin -- ヤコビの反復法 jacobi(a,b,c); Put_Line("X"); disp_vector(c); end Ada1001;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1002 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (20.0, 16.0, 8.0, 17.0); c:Long_Float_Array := ( 0.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- ガウス・ザイデル法 procedure gauss(a:Long_Float_TwoDimArray; b:Long_Float_Array; x0:in out Long_Float_Array) is x1:Long_Float; finish:Boolean; begin while true loop finish := true; for i in x0'Range loop x1 := 0.0; for j in x0'Range loop if j /= i then x1 := x1 + a(i,j) * x0(j); end if; end loop; x1 := (b(i) - x1) / a(i,i); if (Abs(x1 - x0(i)) > 0.0000000001) then finish := false; end if; x0(i) := x1; end loop; if finish then exit; end if; disp_vector(x0); end loop; end gauss; begin -- ガウス・ザイデル法 gauss(a,b,c); Put_Line("X"); disp_vector(c); end Ada1002;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1003 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (8.0, 17.0, 20.0, 16.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for i in matrix'Range loop for j in matrix'Range loop Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- 前進消去 procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is s:Long_Float; begin for pivot in a'First..(a'Last - 1) loop for row in (pivot + 1)..a'Last loop s := a(row,pivot) / a(pivot,pivot); for col in pivot..a'Last loop a(row,col) := a(row,col) - a(pivot,col) * s; end loop; b(row) := b(row) - b(pivot) * s; end loop; end loop; end forward_elimination; -- 後退代入 procedure backward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array) is begin for row in reverse a'Range loop for col in reverse (row + 1)..a'Last loop b(row) := b(row) - a(row,col) * b(col); end loop; b(row) := b(row) / a(row,row); end loop; end backward_substitution; -- ピボット選択 procedure pivoting(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is max_val, tmp:Long_Float; max_row:Integer; begin for pivot in a'Range loop -- 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0.0; for row in pivot..a'Last loop if Abs(a(row,pivot)) > max_val then -- 一番値が大きい行 max_val := Abs(a(row,pivot)); max_row := row; end if; end loop; -- 一番値が大きい行と入れ替え if max_row /= pivot then tmp := 0.0; for col in a'Range loop tmp := a(max_row,col); a(max_row,col) := a(pivot,col); a(pivot,col) := tmp; end loop; tmp := b(max_row); b(max_row) := b(pivot); b(pivot) := tmp; end if; end loop; end pivoting; begin -- ピボット選択 pivoting(a,b); Put_Line("pivoting"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 前進消去 forward_elimination(a,b); Put_Line("forward elimination"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 後退代入 backward_substitution(a,b); Put_Line("X"); disp_vector(b); end Ada1003;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1004 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (8.0, 17.0, 20.0, 16.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for i in matrix'Range loop for j in matrix'Range loop Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- 前進消去 procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is s:Long_Float; begin for pivot in a'Range loop for row in a'Range loop if row /= pivot then s := a(row,pivot) / a(pivot,pivot); for col in pivot..a'Last loop a(row,col) := a(row,col) - a(pivot,col) * s; end loop; b(row) := b(row) - b(pivot) * s; end if; end loop; end loop; end forward_elimination; -- 後退代入 procedure backward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array) is begin for pivot in a'Range loop b(pivot) := b(pivot) / a(pivot,pivot); end loop; end backward_substitution; -- ピボット選択 procedure pivoting(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is max_val, tmp:Long_Float; max_row:Integer; begin for pivot in a'Range loop -- 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0.0; for row in pivot..a'Last loop if Abs(a(row,pivot)) > max_val then -- 一番値が大きい行 max_val := Abs(a(row,pivot)); max_row := row; end if; end loop; -- 一番値が大きい行と入れ替え if max_row /= pivot then tmp := 0.0; for col in a'Range loop tmp := a(max_row,col); a(max_row,col) := a(pivot,col); a(pivot,col) := tmp; end loop; tmp := b(max_row); b(max_row) := b(pivot); b(pivot) := tmp; end if; end loop; end pivoting; begin -- ピボット選択 pivoting(a,b); Put_Line("pivoting"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 前進消去 forward_elimination(a,b); Put_Line("forward elimination"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 後退代入 backward_substitution(a,b); Put_Line("X"); disp_vector(b); end Ada1004;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1005 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (8.0, 17.0, 20.0, 16.0); y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); x:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for i in matrix'Range loop for j in matrix'Range loop Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- 前進消去 procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is s:Long_Float; begin for pivot in a'First..(a'Last - 1) loop for row in (pivot + 1)..a'Last loop s := a(row,pivot) / a(pivot,pivot); for col in pivot..a'Last loop a(row,col) := a(row,col) - a(pivot,col) * s; -- これが 上三角行列 end loop; a(row,pivot) := s; -- これが 下三角行列 -- b(row) := b(row) - b(pivot) * s; -- この値は変更しない end loop; end loop; end forward_elimination; -- 前進代入 procedure forward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array; y:in out Long_Float_Array) is begin for row in a'Range loop for col in a'First..row loop b(row) := b(row) - a(row,col) * y(col); end loop; y(row) := b(row); end loop; end forward_substitution; -- 後退代入 procedure backward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; x:in out Long_Float_Array) is begin for row in reverse a'Range loop for col in reverse (row + 1)..a'Last loop y(row) := y(row) - a(row,col) * x(col); end loop; x(row) := y(row) / a(row,row); end loop; end backward_substitution; -- ピボット選択 procedure pivoting(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is max_val, tmp:Long_Float; max_row:Integer; begin for pivot in a'Range loop -- 各列で 一番値が大きい行を 探す max_row := pivot; max_val := 0.0; for row in pivot..a'Last loop if Abs(a(row,pivot)) > max_val then -- 一番値が大きい行 max_val := Abs(a(row,pivot)); max_row := row; end if; end loop; -- 一番値が大きい行と入れ替え if max_row /= pivot then tmp := 0.0; for col in a'Range loop tmp := a(max_row,col); a(max_row,col) := a(pivot,col); a(pivot,col) := tmp; end loop; tmp := b(max_row); b(max_row) := b(pivot); b(pivot) := tmp; end if; end loop; end pivoting; begin -- ピボット選択 pivoting(a,b); Put_Line("pivoting"); Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 前進消去 forward_elimination(a,b); Put_Line("LU"); disp_matrix(a); -- Ly=b から y を求める (前進代入) forward_substitution(a,b,y); Put_Line("Y"); disp_vector(y); -- Ux=y から x を求める (後退代入) backward_substitution(a,y,x); Put_Line("X"); disp_vector(x); end Ada1005;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1006 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (34.0,68.0,96.0,125.0); y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); x:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for i in matrix'Range loop for j in matrix'Range loop Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- 前進消去 procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is s:Long_Float; begin for pivot in a'Range loop s := 0.0; for col in a'First..(pivot - 1) loop s := s + a(pivot,col) * a(pivot,col); end loop; -- ここで根号の中が負の値になると計算できない! a(pivot,pivot) := Sqrt(a(pivot,pivot) - s); for row in (pivot + 1)..a'Last loop s := 0.0; for col in a'First..(pivot - 1) loop s := s + a(row,col) * a(pivot,col); end loop; a(row,pivot) := (a(row,pivot) - s) / a(pivot,pivot); a(pivot,row) := a(row,pivot); end loop; end loop; end forward_elimination; -- 前進代入 procedure forward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array; y:in out Long_Float_Array) is begin for row in a'Range loop for col in a'First..row loop b(row) := b(row) - a(row,col) * y(col); end loop; y(row) := b(row) / a(row,row); end loop; end forward_substitution; -- 後退代入 procedure backward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; x:in out Long_Float_Array) is begin for row in reverse a'Range loop for col in reverse (row + 1)..a'Last loop y(row) := y(row) - a(row,col) * x(col); end loop; x(row) := y(row) / a(row,row); end loop; end backward_substitution; begin Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 前進消去 forward_elimination(a,b); Put_Line("LL^T"); disp_matrix(a); -- Ly=b から y を求める (前進代入) forward_substitution(a,b,y); Put_Line("Y"); disp_vector(y); -- L^Tx=y から x を求める (後退代入) backward_substitution(a,y,x); Put_Line("X"); disp_vector(x); end Ada1006;
with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1007 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (34.0,68.0,96.0,125.0); y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); x:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for i in matrix'Range loop for j in matrix'Range loop Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- 前進消去 procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is s:Long_Float; begin for pivot in a'Range loop -- pivot < k の場合 s := 0.0; for col in a'First..(pivot - 1) loop s := a(pivot,col); for k in a'First..(col - 1) loop s := s - a(pivot,k) * a(col,k) * a(k,k); end loop; a(pivot,col) := s / a(col,col); a(col,pivot) := a(pivot,col); end loop; -- pivot == k の場合 s := a(pivot,pivot); for k in a'First..(pivot - 1) loop s := s - a(pivot,k) * a(pivot,k) * a(k,k); end loop; a(pivot,pivot) := s; end loop; end forward_elimination; -- 前進代入 procedure forward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array; y:in out Long_Float_Array) is begin for row in a'Range loop for col in a'First..row loop b(row) := b(row) - a(row,col) * y(col); end loop; y(row) := b(row); end loop; end forward_substitution; -- 後退代入 procedure backward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; x:in out Long_Float_Array) is begin for row in reverse a'Range loop for col in reverse (row + 1)..a'Last loop y(row) := y(row) - a(row,col) * a(row,row) * x(col); end loop; x(row) := y(row) / a(row,row); end loop; end backward_substitution; begin Put_Line("A"); disp_matrix(a); Put_Line("B"); disp_vector(b); Put_Line(""); -- 前進消去 forward_elimination(a,b); Put_Line("LDL^T"); disp_matrix(a); -- Ly=b から y を求める (前進代入) forward_substitution(a,b,y); Put_Line("Y"); disp_vector(y); -- DL^Tx=y から x を求める (後退代入) backward_substitution(a,y,x); Put_Line("X"); disp_vector(x); end Ada1007;
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1101 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (1.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 正規化 (ベクトルの長さを1にする) procedure normarize(x:in out Long_Float_Array) is s:Long_Float; begin s := 0.0; for i in x'Range loop s := s + x(i) * x(i); end loop; s := Sqrt(s); for i in x'Range loop x(i) := x(i) / s; end loop; end normarize; -- ベキ乗法 procedure power(a:Long_Float_TwoDimArray; x0:in out Long_Float_Array; lambda:out Long_Float) is p0, p1, e0, e1:Long_Float; x1:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); begin -- 正規化 (ベクトル x0 の長さを1にする) normarize(x0); e0 := 0.0; for i in x0'Range loop e0 := e0 + x0(i); end loop; for k in 1..100 loop -- 1次元配列を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_vector(x0); -- 行列の積 x1 = A × x0 for i in x1'Range loop x1(i) := 0.0; end loop; for i in x1'Range loop for j in x1'Range loop x1(i) := x1(i) + a(i, j) * x0(j); end loop; end loop; -- 内積 p0 := 0.0; p1 := 0.0; for i in x0'Range loop p0 := p0 + x1(i) * x1(i); p1 := p1 + x1(i) * x0(i); end loop; -- 固有値 lambda := p0 / p1; -- 正規化 (ベクトル x1 の長さを1にする) normarize(x1); -- 収束判定 e1 := 0.0; for i in x0'Range loop e1 := e1 + x1(i); end loop; if Abs(e1 - e0) < 0.00000000001 then exit; end if; for i in x0'Range loop x0(i) := x1(i); end loop; e0 := e1; end loop; end power; -- ベキ乗法で最大固有値を求める lambda:Long_Float := 0.0; begin -- ベキ乗法 power(a, x, lambda); Put_Line(""); Put_Line("eigenvalue"); Put(lambda, Fore=>3, Aft=>10, Exp=>0); New_Line; Put_Line("eigenvector"); disp_vector(x); end Ada1101;
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1102 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_Array := (1.0, 0.0, 0.0, 0.0); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 正規化 (ベクトルの長さを1にする) procedure normarize(x:in out Long_Float_Array) is s:Long_Float; begin s := 0.0; for i in x'Range loop s := s + x(i) * x(i); end loop; s := Sqrt(s); for i in x'Range loop x(i) := x(i) / s; end loop; end normarize; -- LU分解 procedure forward_elimination(a:in out Long_Float_TwoDimArray) is s:Long_Float; begin for pivot in a'First..(a'Last - 1) loop for row in (pivot + 1)..a'Last loop s := a(row,pivot) / a(pivot,pivot); for col in pivot..a'Last loop a(row,col) := a(row,col) - a(pivot,col) * s; -- これが 上三角行列 end loop; a(row,pivot) := s; -- これが 下三角行列 end loop; end loop; end forward_elimination; -- 前進代入 (Ly = b から y を求める) procedure forward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; b:in out Long_Float_Array) is begin for row in a'Range loop for col in a'First..row loop b(row) := b(row) - a(row,col) * y(col); end loop; y(row) := b(row); end loop; end forward_substitution; -- 後退代入 (Ux = y から x を求める) procedure backward_substitution(a:Long_Float_TwoDimArray; x:in out Long_Float_Array; y:in out Long_Float_Array) is begin for row in reverse a'Range loop for col in reverse (row + 1)..a'Last loop y(row) := y(row) - a(row,col) * x(col); end loop; x(row) := y(row) / a(row,row); end loop; end backward_substitution; -- 逆ベキ乗法 procedure inverse(a:Long_Float_TwoDimArray; x0:in out Long_Float_Array; lambda:out Long_Float) is p0, p1, e0, e1:Long_Float; b:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); x1:Long_Float_Array := (0.0, 0.0, 0.0, 0.0); begin -- 正規化 (ベクトル x0 の長さを1にする) normarize(x0); e0 := 0.0; for i in x0'Range loop e0 := e0 + x0(i); end loop; for k in 1..100 loop -- 1次元配列を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_vector(x0); -- Ly = b から y を求める (前進代入) for i in x0'Range loop b(i) := x0(i); y(i) := 0.0; end loop; forward_substitution(a,y,b); -- Ux = y から x を求める (後退代入) for i in x1'Range loop x1(i) := 0.0; end loop; backward_substitution(a,x1,y); -- 内積 p0 := 0.0; p1 := 0.0; for i in x0'Range loop p0 := p0 + x1(i) * x1(i); p1 := p1 + x1(i) * x0(i); end loop; -- 固有値 lambda := p1 / p0; -- 正規化 (ベクトル x1 の長さを1にする) normarize(x1); -- 収束判定 e1 := 0.0; for i in x0'Range loop e1 := e1 + x1(i); end loop; if Abs(e1 - e0) < 0.00000000001 then exit; end if; for i in x0'Range loop x0(i) := x1(i); end loop; e0 := e1; end loop; end inverse; -- 逆ベキ乗法で最小固有値を求める lambda:Long_Float := 0.0; begin -- LU分解 forward_elimination(a); -- 逆ベキ乗法 inverse(a, x, lambda); Put_Line(""); Put_Line("eigenvalue"); Put(lambda, Fore=>3, Aft=>10, Exp=>0); New_Line; Put_Line("eigenvector"); disp_vector(x); end Ada1102;
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1103 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_TwoDimArray; u:Long_Float_TwoDimArray; -- 対角要素を表示 procedure disp_eigenvalue(a:Long_Float_TwoDimArray) is begin for i in 0..N loop Put(a(i, i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; Put_Line(""); end disp_eigenvalue; -- 行列の積 procedure multiply(a:Long_Float_TwoDimArray; b:Long_Float_TwoDimArray; c:in out Long_Float_TwoDimArray) is s:Long_Float; begin for i in 0..N loop for j in 0..N loop s := 0.0; for k in 0..N loop s := s + a(i, k) * b(k, j); end loop; c(i, j) := s; end loop; end loop; end multiply; -- LU分解 procedure decomp(a:Long_Float_TwoDimArray; l:in out Long_Float_TwoDimArray; u:in out Long_Float_TwoDimArray) is t:Long_Float; begin for i in 0..N loop for j in 0..N loop l(i, j) := 0.0; u(i, j) := 0.0; end loop; end loop; l(0, 0) := 1.0; for j in 0..N loop u(0, j) := a(0, j); end loop; for i in 1..N loop u(i, 0) := 0.0; l(0, i) := 0.0; l(i, 0) := a(i, 0) / u(0, 0); end loop; for i in 1..N loop l(i, i) := 1.0; t := a(i, i); for k in 0..i loop t := t - l(i, k) * u(k, i); end loop; u(i, i) := t; for j in (i + 1)..N loop u(j, i) := 0.0; l(i, j) := 0.0; t := a(j, i); for k in 0..i loop t := t - l(j, k) * u(k, i); end loop; l(j, i) := t / u(i, i); t := a(i, j); for k in 0..i loop t := t - l(i, k) * u(k, j); end loop; u(i, j) := t; end loop; end loop; end decomp; -- LR分解で固有値を求める e:Long_Float := 0.0; begin for k in 1..200 loop -- LU分解 decomp(a, l, u); -- 行列の積 multiply(u, l, a); -- 対角要素を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_eigenvalue(a); -- 収束判定 e := 0.0; for i in 1..N loop for j in 0..(i - 1) loop e := e + Abs(a(i, j)); end loop; end loop; if e < 0.00000000001 then exit; end if; end loop; Put_Line(""); Put_Line("eigenvalue"); disp_eigenvalue(a); end Ada1103;
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1104 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_TwoDimArray; r:Long_Float_TwoDimArray; -- 対角要素を表示 procedure disp_eigenvalue(a:Long_Float_TwoDimArray) is begin for i in 0..N loop Put(a(i, i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; Put_Line(""); end disp_eigenvalue; -- 行列の積 procedure multiply(a:Long_Float_TwoDimArray; b:Long_Float_TwoDimArray; c:in out Long_Float_TwoDimArray) is s:Long_Float; begin for i in 0..N loop for j in 0..N loop s := 0.0; for k in 0..N loop s := s + a(i, k) * b(k, j); end loop; c(i, j) := s; end loop; end loop; end multiply; -- QR分解 procedure decomp(a:Long_Float_TwoDimArray; q:in out Long_Float_TwoDimArray; r:in out Long_Float_TwoDimArray) is t, s:Long_Float; x:Long_Float_Array; begin for k in 0..N loop for i in 0..N loop x(i) := a(i, k); end loop; for j in 0..(k - 1) loop t := 0.0; for i in 0..N loop t := t + a(i, k) * q(i, j); end loop; r(j, k) := t; r(k, j) := 0.0; for i in 0..N loop x(i) := x(i) - t * q(i, j); end loop; end loop; s := 0.0; for i in 0..N loop s := s + x(i) * x(i); end loop; r(k, k) := Sqrt(s); for i in 0..N loop q(i, k) := x(i) / r(k, k); end loop; end loop; end decomp; -- QR分解で固有値を求める e:Long_Float := 0.0; begin for k in 1..200 loop -- QR分解 decomp(a, q, r); -- 行列の積 multiply(r, q, a); -- 対角要素を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_eigenvalue(a); -- 収束判定 e := 0.0; for i in 1..N loop for j in 0..(i - 1) loop e := e + Abs(a(i, j)); end loop; end loop; if e < 0.00000000001 then exit; end if; end loop; Put_Line(""); Put_Line("eigenvalue"); disp_eigenvalue(a); end Ada1104;
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1105 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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:Long_Float_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)); -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 正規化 (ベクトルの長さを1にする) procedure normarize(x:in out Long_Float_Array) is s:Long_Float; begin s := 0.0; for i in x'Range loop s := s + x(i) * x(i); end loop; s := Sqrt(s); for i in x'Range loop x(i) := x(i) / s; end loop; end normarize; -- 対角要素を表示 procedure disp_eigenvalue(a:Long_Float_TwoDimArray) is begin for i in 0..N loop Put(a(i, i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_eigenvalue; -- 固有ベクトルを表示 procedure disp_eigenvector(matrix:Long_Float_TwoDimArray) is row:Long_Float_Array; begin for i in 0..N loop for j in 0..N loop row(j) := matrix(i, j); end loop; normarize(row); disp_vector(row); end loop; end disp_eigenvector; -- ヤコビ法 procedure jacobi(a:in out Long_Float_TwoDimArray; v:in out Long_Float_TwoDimArray) is p, q:Integer; max_val:Long_Float; c, s, t, t1, t2:Long_Float; begin for k in 1..100 loop -- 最大値を探す max_val := 0.0; for i in 0..N loop for j in (i + 1)..N loop if (max_val < Abs(a(i, j))) then max_val := Abs(a(i, j)); p := i; q := j; end if; end loop; end loop; -- θ を求める t := 0.0; if Abs(a(p, p) - a(q, q)) < 0.00000000001 then -- a_{pp} = a_{qq} のとき、回転角tをπ/4にする t := Ada.Numerics.Pi / 4.0; if (a(p, p) < 0.0) then t := -t; end if; else -- a_{pp} ≠ a_{qq} のとき t := Arctan(2.0 * a(p, q) / (a(p, p) - a(q, q))) / 2.0; end if; -- θ を使って 行列 U を作成し、A = U^t × A × U c := Cos(t); s := Sin(t); -- U^t × A t1 := 0.0; t2 := 0.0; for i in 0..N loop 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 loop; -- A × U for i in 0..N loop 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 loop; --対角要素を表示 Put(k, Width=> 3); Put(Ascii.HT); disp_eigenvalue(a); -- 収束判定 if max_val < 0.00000000001 then exit; end if; end loop; end jacobi; -- ヤコビ法で固有値を求める begin -- ヤコビ法 jacobi(a, v); Put_Line(""); Put_Line("eigenvalue"); disp_eigenvalue(a); New_Line; Put_Line("eigenvector"); disp_eigenvector(v); end Ada1105;
with TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; use TEXT_IO, Ada.Integer_Text_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions; procedure Ada1106 is N : Constant Integer := 3; type Long_Float_Array is array (0..N) of Long_Float; type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float; a:Long_Float_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, e:Long_Float_Array; -- 1次元配列を表示 procedure disp_vector(row:Long_Float_Array) is begin for i in row'Range loop Put(row(i), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end disp_vector; -- 2次元配列を表示 procedure disp_matrix(matrix:Long_Float_TwoDimArray) is begin for row in 0..N loop for col in 0..N loop Put(matrix(row, col), Fore=>3, Aft=>10, Exp=>0); Put(Ascii.HT); end loop; New_Line; end loop; end disp_matrix; -- ハウスホルダー変換 procedure tridiagonalize(a:in out Long_Float_TwoDimArray; d:in out Long_Float_Array; e:in out Long_Float_Array) is w, v:Long_Float_Array; t, s:Long_Float; begin for k in 0..N loop v(k) := 0.0; end loop; for k in 0..(N - 2) loop for i in 0..N loop w(i) := 0.0; end loop; d(k) := a(k, k); t := 0.0; for i in (k + 1)..N loop w(i) := a(i, k); t := t + w(i) * w(i); end loop; s := Sqrt(t); if w(k + 1) < 0.0 then s := -s; end if; if Abs(s) < 0.00000000001 then e(k + 1) := 0.0; else e(k + 1) := -s; w(k + 1) := w(k + 1) + s; s := 1.0 / Sqrt(w(k + 1) * s); for i in (k + 1)..N loop w(i) := w(i) * s; end loop; for i in (k + 1)..N loop s := 0.0; for j in (k + 1)..N loop if j <= i then s := s + a(i, j) * w(j); else s := s + a(j, i) * w(j); end if; end loop; v(i) := s; end loop; s := 0.0; for i in (k + 1)..N loop s := s + w(i) * v(i); end loop; s := s / 2.0; for i in (k + 1)..N loop v(i) := v(i) - s * w(i); end loop; for i in (k + 1)..N loop for j in (k + 1)..i loop a(i, j) := a(i, j) - (w(i) * v(j) + w(j) * v(i)); end loop; end loop; -- 次の行は固有ベクトルを求めないなら不要 for i in (k + 1)..N loop a(i, k) := w(i); end loop; end if; end loop; 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 in reverse 0..N loop for i in 0..N loop w(i) := 0.0; end loop; if k < N - 1 then for i in (k + 1)..N loop w(i) := a(i, k); end loop; for i in (k + 1)..N loop s := 0.0; for j in (k + 1)..N loop s := s + a(i, j) * w(j); end loop; v(i) := s; end loop; for i in (k + 1)..N loop for j in (k + 1)..N loop a(i, j) := a(i, j) - v(i) * w(j); end loop; end loop; end if; for i in 0..N loop a(i, k) := 0.0; end loop; a(k, k) := 1.0; end loop; end tridiagonalize; -- QR分解 procedure decomp(a:in out Long_Float_TwoDimArray; d:in out Long_Float_Array; e:in out Long_Float_Array) is h, l:Integer; s, t, u, w, x, y, z:Long_Float; begin e(0) := 1.0; h := N; while (Abs(e(h)) < 0.00000000001) loop h := h -1; end loop; while (h > 0) loop e(0) := 0.0; l := h - 1; while (Abs(e(l)) >= 0.00000000001) loop l := l - 1; end loop; for j in 1..100 loop w := (d(h - 1) - d(h)) / 2.0; s := Sqrt(w * w + e(h) * e(h)); if w < 0.0 then s := -s; end if; 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 in l..(h - 1) loop if Abs(x) >= Abs(y) then t := -y / x; u := 1.0 / Sqrt(t * t + 1.0); s := t * u; else t := -x / y; s := 1.0 / Sqrt(t * t + 1.0); if t < 0.0 then s := -s; end if; u := t * s; end if; 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 in 0..N loop 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 loop; if k < N - 1 then x := e(k + 1); y := -s * e(k + 2); z := y; e(k + 2) := u * e(k + 2); end if; end loop; Put(j, Width=> 3); Put(Ascii.HT); disp_vector(d); -- 収束判定 if (Abs(e(h)) < 0.00000000001) then exit; end if; end loop; e(0) := 1.0; while (Abs(e(h)) < 0.00000000001) loop h := h - 1; end loop; end loop; -- 次の行は固有ベクトルを求めないなら不要 for k in 0..(N - 1) loop l := k; for i in (k + 1)..N loop if d(i) > d(l) then l := i; end if; end loop; t := d(k); d(k) := d(l); d(l) := t; for i in 0..N loop t := a(k, i); a(k, i) := a(l, i); a(l, i) := t; end loop; end loop; end decomp; -- ハウスホルダー変換とQR分解で固有値を求める begin -- ハウスホルダー変換; tridiagonalize(a, d, e); -- QR分解 decomp(a, d, e); New_Line; Put_Line("eigenvalue"); disp_vector(d); New_Line; Put_Line("eigenvector"); disp_matrix(a); end Ada1106;