module Fs0101 printfn "%d" (3 + 5) printfn "%d" (3 - 5) printfn "%d" (3 * 5) printfn "%f" (System.Math.Pow(float 3, float 5)) printfn "%d" (5 / 3) printfn "%f" (5.0 / 3.0) printfn "%d" (5 % 3) printf "%d\n" (3 * 5) printfn "%3d" (3 * 5) printfn "%23.20f" (5.0 / 3.0) exit 0
module Fs0102 let i = 3 * 5 printfn "3 * 5 = %d" i exit 0
module Fs0103 for i in [1..9] do printf "%d, " i printfn "" exit 0
module Fs0104 for i in [1..9] do if (i % 3 = 0) then printf "%d, " i printfn "" exit 0
module Fs0105 let mutable sm = 0 for i in [1..99] do if (i % 3 = 0) then sm <- sm + i printfn "%d" sm exit 0
[1..9]
[1..9] |> List.filter (fun x -> x % 3 = 0)
[1..99] |> List.filter (fun x -> x % 3 = 0) |> List.reduce(+) [1..99] |> List.filter (fun x -> x % 3 = 0) |> List.sum
// 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 let sn a lim = let n = lim / a in // 項数:n = 上限:lim / 公差:a let l = a * n in // 末項:l = 項数:n * 公差:a (a + l) * n / 2 // 総和:sn = (初項:a + 末項:l) * 項数:n / 2 // 3 の倍数の合計 (sn 3 999)
let n = 10000
n * (n + 1) / 2
let n = 5000
n * (n + 1)
let n = 5000 System.Math.Pow(double n, 2.0)
let n = 1000
n * (n + 1) * (2 * n + 1) / 6
let n = 100 System.Math.Pow(double n, 2.0) * System.Math.Pow((double n + 1.0), 2.0) / 4.0
let n = 10 let a = 2 let r = 3 (a * (int (System.Math.Pow(double r, double n)) - 1)) / (r - 1)
[0..9] [0..9] |> List.map (fun n -> n * 3 + 5) [0L..9L] |> List.map (fun n -> n * 3L + 5L) |> List.reduce(*)
// 等差数列の積 let rec product (m:int64) (d:int64) (n:int64) :int64 = match n with | 0L -> 1L | _ -> m * (product (m + d) d (n - 1L)) // 初項 5, 公差 3, 項数 10 の数列の積 product 5L 3L 10L
// 階乗を求める関数 let rec Fact = function | 0 -> 1 | n -> n * Fact(n - 1) // 10の階乗 Fact 10
// 下降階乗冪 let rec FallingFact (x:int) (n:int):int = match n with | 1 -> x | _ -> x * (FallingFact (x - 1) (n - 1)) // 10 から 6 までの 総乗 FallingFact 10 5
// 上昇階乗冪 let rec RisingFact (x:int) (n:int):int = match n with | 1 -> x | _ -> x * (RisingFact (x + 1) (n - 1)) // 10 から 14 までの 総乗 RisingFact 10 5
// 階乗 let rec Fact = function | 0 -> 1 | n -> n * Fact(n - 1) // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) let n = 10 let r = 5 Fact n / Fact (n - r) // 下降階乗冪 let rec FallingFact (x:int) (n:int):int = match n with | 1 -> x | _ -> x * (FallingFact (x - 1) (n - 1)) // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) let n = 10 let r = 5 FallingFact n r
// 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) let n = 10 let r = 5 System.Math.Pow(float n, float r)
// 組合せ let rec Comb(n:int) (r:int):int = match n, r with | _, 0 -> 1 | _, 1 -> n | _, _ when n = r -> 1 | _, _ -> Comb (n - 1) (r - 1) + Comb(n - 1) r // 異なる 10 個のものから 5 個取ってできる組合せの総数 let n = 10 let r = 5 Comb n r
// 組合せ let rec Comb(n:int) (r:int):int = match n, r with | _, 0 -> 1 | _, 1 -> n | _, _ when n = r -> 1 | _, _ -> Comb (n - 1) (r - 1) + Comb(n - 1) r // 異なる 10 個のものから重複を許して 5 個取る組合せの総数 let n = 10 let r = 5 Comb (n + r - 1) r
// 自作の正弦関数 let rec mySin (x:double) (n:int) (nega:bool) (numerator:double) (denominator:double) (y:double):double = let m = 2 * n let denom = denominator * (double (m + 1)) * (double m) let num = numerator * x * x let a = num / denom // 十分な精度になったら処理を抜ける if a <= 0.00000000001 then y else y + (mySin x (n + 1) (not nega) num denom (if nega then a else -a) ) [0..24] |> List.map (fun n -> n * 15) |> List.filter (fun n -> (n % 30 = 0) || (n % 45 = 0)) |> List.iter (fun degree -> let radian = (double degree) * System.Math.PI / 180.0 // 自作の正弦関数 let d1 = (mySin radian 1 false radian 1.0 radian) // 標準の正弦関数 let d2 = System.Math.Sin(radian) // 標準関数との差異 printfn "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (d1 - d2) )
// 自作の余弦関数 let rec myCos (x:double) (n:int) (nega:bool) (numerator:double) (denominator:double) (y:double):double = let m = 2 * n let denom = denominator * (double m) * (double (m - 1)) let num = numerator * x * x let a = num / denom // 十分な精度になったら処理を抜ける if a <= 0.00000000001 then y else y + (myCos x (n + 1) (not nega) num denom (if nega then a else -a) ) [0..24] |> List.map ((*) 15) |> List.filter (fun n -> (n % 30 = 0) || (n % 45 = 0)) |> List.iter (fun degree -> let radian = (double degree) * System.Math.PI / 180.0 // 自作の余弦関数 let d1 = (myCos radian 1 false 1.0 1.0 1.0) // 標準の余弦関数 let d2 = System.Math.Cos(radian) // 標準関数との差異 printfn "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (d1 - d2) )
// 自作の正接関数 let rec myTan (x:double) (x2:double) (n:int) (t:double):double = let denom = x2 / ((double n) - t) let num = n - 2 if num <= 1 then x / (1.0 - denom) else myTan x x2 num denom [0..12] |> List.map ((*) 15) |> List.filter (fun n -> n % 180 <> 0) |> List.map ((-) 90) |> List.iter (fun degree -> let radian = (double degree) * System.Math.PI / 180.0 let x2 = radian * radian // 自作の正接関数 let d1 = (myTan radian x2 15 0.0) // 15:必要な精度が得られる十分大きな奇数 // 標準の正接関数 let d2 = System.Math.Tan(radian) // 標準関数との差異 printfn "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (d1 - d2) )
// 自作の指数関数 let rec myExp (x:double) (n:int) (numerator:double) (denominator:double) (y:double):double = let denom = denominator * (double n) let num = numerator * x let a = num / denom // 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (myExp x (n + 1) num denom a) [0..20] |> List.map (fun n -> (double (n - 10)) / 4.0) |> List.iter (fun x -> // 標準の指数関数 let d1 = System.Math.Exp(x) // 自作の指数関数 let d2 = (myExp x 1 1.0 1.0 1.0) // 標準関数との差異 printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2) )
// 自作の指数関数 let rec myExp (x:double) (x2:double) (n:int) (t:double):double = let denom = x2 / ((double n) + t) let num = n - 4; if num < 6 then 1.0 + ((2.0 * x) / (2.0 - x + denom)) else myExp x x2 num denom [0..20] |> List.map (fun n -> (double (n - 10)) / 4.0) |> List.iter (fun x -> // 標準の指数関数 let d1 = System.Math.Exp(x) // 自作の指数関数 let x2 = x * x let d2 = (myExp x x2 30 0.0) // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる // 標準関数との差異 printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2) )
// 自作の対数関数 let rec myLog (x2:double) (numerator:double) (denominator:double) (y:double):double = let denom = denominator + 2.0 let num = numerator * x2 * x2 let a = num / denom // 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (myLog x2 num denom a) [1..20] |> List.map (fun n -> (double n) / 5.0) |> List.iter (fun x -> // 標準の対数関数 let d1 = System.Math.Log(x) // 自作の対数関数 let x2 = (x - 1.0) / (x + 1.0) let d2 = 2.0 * (myLog x2 x2 1.0 x2) // 標準関数との差異 printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2) )
// 自作の対数関数 let rec myLog (x:double) (n:int) (t:double):double = let mutable n2 = n let mutable x2 = x if n > 3 then if n % 2 = 0 then n2 <- 2 x2 <- x * (double (n / 2)) let t2 = x2 / ((double n2) + t) if n <= 2 then x / (1.0 + t2) else myLog x (n - 1) t2 [1..20] |> List.map (fun n -> (double n) / 5.0) |> List.iter (fun x -> // 標準の対数関数 let d1 = System.Math.Log(x) // 自作の対数関数 let d2 = (myLog (x - 1.0) 27 0.0) // 27:必要な精度が得られる十分大きな奇数 // 標準関数との差異 printfn "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (d1 - d2) )
// 自作の双曲線正弦関数 let rec mySinh (x:double) (n:int) (numerator:double) (denominator:double) (y:double):double = let m = 2 * n let denom = denominator * (double (m + 1)) * (double m) let num = numerator * x * x let a = num / denom // 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (mySinh x (n + 1) num denom a) [0..20] |> List.map((-) 10) |> List.iter (fun i -> // 自作の双曲線正弦関数 let x = (double i) let d1 = (mySinh x 1 x 1.0 x) // 標準の双曲線正弦関数 let d2 = System.Math.Sinh(x) // 標準関数との差異 printfn "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (d1 - d2) )
// 自作の双曲線余弦関数 let rec myCosh (x:double) (n:int) (numerator:double) (denominator:double) (y:double):double = let m = 2 * n let denom = denominator * (double m) * (double (m - 1)) let num = numerator * x * x let a = num / denom // 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (myCosh x (n + 1) num denom a) [0..20] |> List.map((-) 10) |> List.iter (fun i -> let x = (double i) // 自作の双曲線余弦関数 let d1 = (myCosh x 1 1.0 1.0 1.0) // 標準の双曲線余弦関数 let d2 = System.Math.Cosh(x) // 標準関数との差異 printfn "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (d1 - d2) )
module Fs0601 open System let f(x:double):double = 4.0 / (1.0 + x * x) let a:double = 0.0 let b:double = 1.0 // 台形則で積分 let mutable n = 2 for j in [1..10] do let h:double = (b - a) / (double n) let mutable s:double = 0.0 let mutable x:double = a for i in [1..(n - 1)] do x <- x + h s <- s + (f x) s <- h * (((f a) + (f b)) / 2.0 + s) n <- n * 2 // 結果を π と比較 printfn "%3d : %13.10f, %13.10f" j s (s - Math.PI) exit 0
module Fs0602 open System let f(x:double):double = 4.0 / (1.0 + x * x) let a:double = 0.0 let b:double = 1.0 // 中点則で積分 let mutable n = 2 for j in [1..10] do let h:double = (b - a) / (double n) let mutable s:double = 0.0 let mutable x:double = a + (h / 2.0) for i in [1..n] do s <- s + (f x) x <- x + h s <- s * h n <- n * 2 // 結果を π と比較 printfn "%3d : %13.10f, %13.10f" j s (s - Math.PI) exit 0
module Fs0603 open System let f(x:double):double = 4.0 / (1.0 + x * x) let a:double = 0.0 let b:double = 1.0 // Simpson則で積分 let mutable n = 2 for j in [1..5] do let h:double = (b - a) / (double n) let mutable s2:double = 0.0 let mutable s4:double = 0.0 let mutable x:double = a + h for i in [1..(n / 2)] do s4 <- s4 + (f x) x <- x + h s2 <- s2 + (f x) x <- x + h s2 <- (s2 - (f b)) * 2.0 + (f a) + (f b) s4 <- s4 * 4.0 let s = (s2 + s4) * h / 3.0 n <- n * 2 // 結果を π と比較 printfn "%3d : %13.10f, %13.10f" j s (s - System.Math.PI) exit 0
module Fs0604 open System let f(x:double):double = 4.0 / (1.0 + x * x) let a:double = 0.0 let b:double = 1.0 let t = Array2D.zeroCreate<double> 7 7 // 台形則で積分 let mutable n = 2 for i in [1..6] do let h:double = (b - a) / (double n) let mutable s:double = 0.0 let mutable x:double = a for j in [1..(n - 1)] do x <- x + h s <- s + (f x) // 結果を保存 t.[i,1] <- h * (((f a) + (f b)) / 2.0 + s) n <- n * 2 // Richardsonの補外法 n <- 4 for j in [2..6] do for i in [j..6] do t.[i,j] <- t.[i, j - 1] + (t.[i, j - 1] - t.[i - 1, j - 1]) / (double (n - 1)) if i = j then // 結果を π と比較 printfn "%3d : %13.10f, %13.10f" j t.[i,j] (t.[i,j] - Math.PI) n <- n * 4 exit 0
module Fs0701 open System // データ点の数 - 1 let N = 6 // 元の関数 let f (x:double):double = x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2)) // Lagrange (ラグランジュ) 補間 let lagrange(d:double) (x:double list) (y:double list) = let mutable sum_list = [0.0] for i in [0..N] do let mutable prod_list = [y.[i]] for j in [0..N] do if i <> j then prod_list <- ((d - x.[j]) / (x.[i] - x.[j]))::prod_list sum_list <- (prod_list |> List.reduce(*))::sum_list sum_list |> List.sum // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5) let y = x |> List.map(f) // 0.5刻みで 与えられていない値を補間 let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5) let d2 = d1 |> List.map(f) let d3 = d1 |> List.map(fun d -> (lagrange d x y)) (List.zip (List.zip d1 d2) d3) |> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3)) exit 0
module Fs0702 open System // データ点の数 - 1 let N = 6 // 元の関数 let f (x:double):double = x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2)) // Neville (ネヴィル) 補間 let neville(d:double) (x:double list) (y:double list) = let w = Array2D.zeroCreate<double> (N+1) (N+1) for i in [0..N] do w.[0,i] <- y.[i] for j in [1..N] do for i in [0..N-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]) w.[N,0] // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5) let y = x |> List.map(f) // 0.5刻みで 与えられていない値を補間 let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5) let d2 = d1 |> List.map(f) let d3 = d1 |> List.map(fun d -> (neville d x y)) (List.zip (List.zip d1 d2) d3) |> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3)) exit 0
module Fs0703 open System // データ点の数 - 1 let N = 6 // 元の関数 let f (x:double):double = x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2)) // Newton (ニュートン) 補間 let newton(d:double) (x:double list) (a:double list) = let mutable sum_list = [a.[0]] for i in [1..N] do let mutable prod_list = [a.[i]] for j in [0..i-1] do prod_list <- (d - x.[j])::prod_list sum_list <- (prod_list |> List.reduce(*))::sum_list sum_list |> List.sum // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5) let y = x |> List.map(f) // 差分商の表を作る let d = Array2D.zeroCreate<double> (N+1) (N+1) for j in [0..N] do d.[0,j] <- y.[j] for i in [1..N] do for j in [0..N-i] do d.[i,j] <- (d.[i-1,j+1] - d.[i-1,j]) / (x.[j+i] - x.[j]) // n階差分商 let a = [0..N] |> List.map(fun i -> d.[i,0]) // 0.5刻みで 与えられていない値を補間 let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5) let d2 = d1 |> List.map(f) let d3 = d1 |> List.map(fun d -> (newton d x a)) (List.zip (List.zip d1 d2) d3) |> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3)) exit 0
module Fs0704 open System // データ点の数 - 1 let N = 6 let Nx2 = 13 // 元の関数 let f (x:double):double = x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2)) // 導関数 let fd (x:double):double = 1.0 - Math.Pow(x,2.0) / (float 2) + Math.Pow(x,4.0) / (float (4 * 3 * 2)) // Hermite (エルミート) 補間 let hermite(d:double) (z:double list) (a:double list) = let mutable sum_list = [a.[0]] for i in [1..Nx2] do let mutable prod_list = [a.[i]] for j in [0..i-1] do prod_list <- (d - z.[j])::prod_list sum_list <- (prod_list |> List.reduce(*))::sum_list sum_list |> List.sum // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5) let y = x |> List.map(f) let yd = x |> List.map(fd) // 差分商の表を作る let z = [0..Nx2] |> List.map(fun i -> i / 2) |> List.map(fun j -> x.[j]) let d = Array2D.zeroCreate<double> (Nx2+1) (Nx2+1) for i in [0..Nx2] do d.[0,i] <- y.[i / 2] for i in [1..Nx2] do for j in [0..Nx2-i] do if (i = 1 && j % 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]) // n階差分商 let a = [0..Nx2] |> List.map(fun i -> d.[i,0]) // 0.5刻みで 与えられていない値を補間 let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5) let d2 = d1 |> List.map(f) let d3 = d1 |> List.map(fun d -> (hermite d z a)) (List.zip (List.zip d1 d2) d3) |> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3)) exit 0
module Fs0705 open System // データ点の数 - 1 let N = 6 // 元の関数 let f (x:double):double = x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2)) // Spline (スプライン) 補間 let spline(d:double) (x:double list) (y:double list) (z:double list) = // 補間関数値がどの区間にあるか let mutable k = -1 for i in [N..(-1)..1] do if d <= x.[i] then k <- i - 1 if k < 0 then k <- N let d1 = x.[k+1] - d let d2 = d - x.[k] let d3 = x.[k+1] - x.[k] (z.[k] * Math.Pow(d1,3.0) + z.[k+1] * Math.Pow(d2,3.0)) / (6.0 * d3) + (y.[k] / d3 - z.[k] * d3 / 6.0) * d1 + (y.[k+1] / d3 - z.[k+1] * d3 / 6.0) * d2 // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5) let y = x |> List.map(f) // 3項方程式の係数の表を作る let a = Array.zeroCreate<double> N let b = Array.zeroCreate<double> N let c = Array.zeroCreate<double> N let d = Array.zeroCreate<double> N for i in [1..N - 1] do a.[i] <- x.[i] - x.[i-1] b.[i] <- 2.0 * (x.[i+1] - x.[i-1]) c.[i] <- x.[i+1] - x.[i] d.[i] <- 6.0 * ((y.[i+1] - y.[i]) / (x.[i+1] - x.[i]) - (y.[i] - y.[i-1]) / (x.[i] - x.[i-1])) // 3項方程式を解く (ト−マス法) let g = Array.zeroCreate<double> N let s = Array.zeroCreate<double> N g.[1] <- b.[1] s.[1] <- d.[1] for i in [2..N - 1] do g.[i] <- b.[i] - a.[i] * c.[i-1] / g.[i-1] s.[i] <- d.[i] - a.[i] * s.[i-1] / g.[i-1] let z = Array.zeroCreate<double> (N+1) z.[0] <- 0.0 z.[N] <- 0.0 z.[N-1] <- s.[N-1] / g.[N-1] for i in [(N-2)..(-1)..1] do z.[i] <- (s.[i] - c.[i] * z.[i+1]) / g.[i] // 0.5刻みで 与えられていない値を補間 let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5) let d2 = d1 |> List.map(f) let d3 = d1 |> List.map(fun d -> (spline d x y (Array.toList(z)))) (List.zip (List.zip d1 d2) d3) |> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3)) exit 0
module Fs0801 open System // 重力加速度 let g = -9.8 // 空気抵抗係数 let k = -0.01 // 時間間隔(秒) let h = 0.01 // 角度 let degree = 45.0 let radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 let v = float(250 * 1000 / 3600) // 水平方向の速度 let vx = v * Math.Cos(radian) // 鉛直方向の速度 let vy = v * Math.Sin(radian) // 位置 let x = 0.0 let y = 0.0 // 空気抵抗による水平方向の減速分 let fx(vx:Double) (vy:Double) = k * Math.Sqrt(vx * vx + vy * vy) * vx // 重力と空気抵抗による鉛直方向の減速分 let fy(vx:Double) (vy:Double) = g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) // Euler法 let rec euler(i:int) (vx:double) (vy:double) (x:double) (y:double):unit = // 経過秒数 let t = float(i) * h // 位置 let wx = x + h * vx let wy = y + h * vy printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t vx vy wx wy // 速度 let wvx = vx + h * (fx vx vy) let wvy = vy + h * (fy vx vy) if wy >= 0.0 then (euler (i+1) wvx wvy wx wy) else () // Euler法 (euler 1 vx vy x y) exit 0
module Fs0802 open System // 重力加速度 let g = -9.8 // 空気抵抗係数 let k = -0.01 // 時間間隔(秒) let h = 0.01 // 角度 let degree = 45.0 let radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 let v = float(250 * 1000 / 3600) // 水平方向の速度 let vx = v * Math.Cos(radian) // 鉛直方向の速度 let vy = v * Math.Sin(radian) // 位置 let x = 0.0 let y = 0.0 // 空気抵抗による水平方向の減速分 let fx(vx:Double) (vy:Double) = k * Math.Sqrt(vx * vx + vy * vy) * vx // 重力と空気抵抗による鉛直方向の減速分 let fy(vx:Double) (vy:Double) = g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) // Heun法 let rec heun(i:int) (vx:double) (vy:double) (x:double) (y:double):unit = // 経過秒数 let t = float(i) * h // 位置・速度 let wx2 = x + h * vx let wy2 = y + h * vy let wvx2 = vx + h * (fx vx vy) let wvy2 = vy + h * (fy vx vy) let wx = x + h * ( vx + wvx2 ) / 2.0 let wy = y + h * ( vy + wvy2 ) / 2.0 let wvx = vx + h * ((fx vx vy) + (fx wvx2 wvy2)) / 2.0 let wvy = vy + h * ((fy vx vy) + (fy wvx2 wvy2)) / 2.0 printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy if wy >= 0.0 then (heun (i+1) wvx wvy wx wy) else () // Heun法 (heun 1 vx vy x y) exit 0
module Fs0803 open System // 重力加速度 let g = -9.8 // 空気抵抗係数 let k = -0.01 // 時間間隔(秒) let h = 0.01 // 角度 let degree = 45.0 let radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 let v = float(250 * 1000 / 3600) // 水平方向の速度 let vx = v * Math.Cos(radian) // 鉛直方向の速度 let vy = v * Math.Sin(radian) // 位置 let x = 0.0 let y = 0.0 // 空気抵抗による水平方向の減速分 let fx(vx:Double) (vy:Double) = k * Math.Sqrt(vx * vx + vy * vy) * vx // 重力と空気抵抗による鉛直方向の減速分 let fy(vx:Double) (vy:Double) = g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) // 中点法 let rec midpoint(i:int) (vx:double) (vy:double) (x:double) (y:double):unit = // 経過秒数 let t = float(i) * h // 位置・速度 let wvx1 = h * (fx vx vy) let wvy1 = h * (fy vx vy) let wvx2 = vx + wvx1 / 2.0 let wvy2 = vy + wvy1 / 2.0 let wvx = vx + h * (fx wvx2 wvy2) let wvy = vy + h * (fy wvx2 wvy2) let wx = x + h * wvx2 let wy = y + h * wvy2 printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy if wy >= 0.0 then (midpoint (i+1) wvx wvy wx wy) else () // 中点法 (midpoint 1 vx vy x y) exit 0
module Fs0804 open System // 重力加速度 let g = -9.8 // 空気抵抗係数 let k = -0.01 // 時間間隔(秒) let h = 0.01 // 角度 let degree = 45.0 let radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 let v = float(250 * 1000 / 3600) // 水平方向の速度 let vx = v * Math.Cos(radian) // 鉛直方向の速度 let vy = v * Math.Sin(radian) // 位置 let x = 0.0 let y = 0.0 // 空気抵抗による水平方向の減速分 let fx(vx:Double) (vy:Double) = k * Math.Sqrt(vx * vx + vy * vy) * vx // 重力と空気抵抗による鉛直方向の減速分 let fy(vx:Double) (vy:Double) = g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) // Runge-Kutta法 let rec rungekutta(i:int) (vx:double) (vy:double) (x:double) (y:double):unit = // 経過秒数 let t = float(i) * h // 位置・速度 let wx1 = h * vx let wy1 = h * vy let wvx1 = h * (fx vx vy) let wvy1 = h * (fy vx vy) let wvx5 = vx + wvx1 / 2.0 let wvy5 = vy + wvy1 / 2.0 let wx2 = h * wvx5 let wy2 = h * wvy5 let wvx2 = h * (fx wvx5 wvy5) let wvy2 = h * (fy wvx5 wvy5) let wvx6 = vx + wvx2 / 2.0 let wvy6 = vy + wvy2 / 2.0 let wx3 = h * wvx6 let wy3 = h * wvy6 let wvx3 = h * (fx wvx6 wvy6) let wvy3 = h * (fy wvx6 wvy6) let wvx7 = vx + wvx3 let wvy7 = vy + wvy3 let wx4 = h * wvx7 let wy4 = h * wvy7 let wvx4 = h * (fx wvx7 wvy7) let wvy4 = h * (fy wvx7 wvy7) let wx = x + ( wx1 + wx2 * 2.0 + wx3 * 2.0 + wx4) / 6.0 let wy = y + ( wy1 + wy2 * 2.0 + wy3 * 2.0 + wy4) / 6.0 let wvx = vx + (wvx1 + wvx2 * 2.0 + wvx3 * 2.0 + wvx4) / 6.0 let wvy = vy + (wvy1 + wvy2 * 2.0 + wvy3 * 2.0 + wvy4) / 6.0 printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy if wy >= 0.0 then (rungekutta (i+1) wvx wvy wx wy) else () // Runge-Kutta法 (rungekutta 1 vx vy x y) exit 0
module Fs0805 open System // 重力加速度 let g = -9.8 // 空気抵抗係数 let k = -0.01 // 時間間隔(秒) let h = 0.01 // 角度 let degree = 45.0 let radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 let v = float(250 * 1000 / 3600) // 水平方向の速度 let vx = v * Math.Cos(radian) // 鉛直方向の速度 let vy = v * Math.Sin(radian) // 位置 let x = 0.0 let y = 0.0 // 空気抵抗による水平方向の減速分 let fx(vx:Double) (vy:Double) = k * Math.Sqrt(vx * vx + vy * vy) * vx // 重力と空気抵抗による鉛直方向の減速分 let fy(vx:Double) (vy:Double) = g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) // Runge-Kutta-Gill法 let rec rungekuttagill(i:int) (vx:double) (vy:double) (x:double) (y:double):unit = // 経過秒数 let t = float(i) * h // 位置・速度 let wx1 = h * vx let wy1 = h * vy let wvx1 = h * (fx vx vy) let wvy1 = h * (fy vx vy) let wvx5 = vx + wvx1 / 2.0 let wvy5 = vy + wvy1 / 2.0 let wx2 = h * wvx5 let wy2 = h * wvy5 let wvx2 = h * (fx wvx5 wvy5) let wvy2 = h * (fy wvx5 wvy5) let wvx6 = vx + wvx1 * ((Math.Sqrt(2.0) - 1.0) / 2.0) + wvx2 * (1.0 - 1.0 / Math.Sqrt(2.0)) let wvy6 = vy + wvy1 * ((Math.Sqrt(2.0) - 1.0) / 2.0) + wvy2 * (1.0 - 1.0 / Math.Sqrt(2.0)) let wx3 = h * wvx6 let wy3 = h * wvy6 let wvx3 = h * (fx wvx6 wvy6) let wvy3 = h * (fy wvx6 wvy6) let wvx7 = vx - wvx2 / Math.Sqrt(2.0) + wvx3 * (1.0 + 1.0 / Math.Sqrt(2.0)) let wvy7 = vy - wvy2 / Math.Sqrt(2.0) + wvy3 * (1.0 + 1.0 / Math.Sqrt(2.0)) let wx4 = h * wvx7 let wy4 = h * wvy7 let wvx4 = h * (fx wvx7 wvy7) let wvy4 = h * (fy wvx7 wvy7) let wx = x + ( wx1 + wx2 * (2.0 - Math.Sqrt(2.0)) + wx3 * (2.0 + Math.Sqrt(2.0)) + wx4) / 6.0 let wy = y + ( wy1 + wy2 * (2.0 - Math.Sqrt(2.0)) + wy3 * (2.0 + Math.Sqrt(2.0)) + wy4) / 6.0 let wvx = vx + (wvx1 + wvx2 * (2.0 - Math.Sqrt(2.0)) + wvx3 * (2.0 + Math.Sqrt(2.0)) + wvx4) / 6.0 let wvy = vy + (wvy1 + wvy2 * (2.0 - Math.Sqrt(2.0)) + wvy3 * (2.0 + Math.Sqrt(2.0)) + wvy4) / 6.0 printfn "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t wvx wvy wx wy if wy >= 0.0 then (rungekuttagill (i+1) wvx wvy wx wy) else () // Runge-Kutta-Gill法 (rungekuttagill 1 vx vy x y) exit 0
module Fs0901 open System let f (x:double):double = x * x - 2.0 let rec bisection (a:double) (b:double):double = // 区間 (a, b) の中点 c = (a + b) / 2 let c = (a + b) / 2.0 printfn "%12.10f\t%13.10f" c (c - Math.Sqrt(2.0)) let fc = f c if abs(fc) < 0.0000000001 then c else if fc < 0.0 then // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 bisection c b else // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 bisection a c let a = 1.0 let b = 2.0 printfn "%12.10f" (bisection a b) exit 0
module Fs0902 open System let f (x:double):double = x * x - 2.0 let rec falseposition (a:double) (b:double):double = // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 let c = (a * f(b) - b * f(a)) / (f(b) - f(a)) printfn "%12.10f\t%13.10f" c (c - Math.Sqrt(2.0)) let fc = f c if abs(fc) < 0.0000000001 then c else if fc < 0.0 then // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 falseposition c b else // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 falseposition a c let a = 1.0 let b = 2.0 printfn "%12.10f" (falseposition a b) exit 0
module Fs0903 open System let g (x:double):double = (x / 2.0) + (1.0 / x) let rec iterative (x0:double):double = let x1 = g(x0) printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0)) if abs(x1 - x0) < 0.0000000001 then x1 else iterative x1 let x = 1.0 printfn "%12.10f" (iterative x) exit 0
module Fs0904 open System let f0 (x:double):double = x * x - 2.0 let f1 (x:double):double = 2.0 * x let rec newton (x0:double):double = let x1 = x0 - (f0(x0) / f1(x0)) printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0)) if abs(x1 - x0) < 0.0000000001 then x1 else newton x1 let x = 2.0 printfn "%12.10f" (newton x) exit 0
module Fs0905 open System let f0 (x:double):double = x * x - 2.0 let f1 (x:double):double = 2.0 * x let f2 (x:double):double = 2.0 let rec bailey (x0:double):double = let x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2.0 * f1(x0))))) printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0)) if abs(x1 - x0) < 0.0000000001 then x1 else bailey x1 let x = 2.0 printfn "%12.10f" (bailey x) exit 0
module Fs0906 open System let f (x:double):double = x * x - 2.0 let rec secant (x0:double)(x1:double):double = let x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) printfn "%12.10f\t%13.10f" x2 (x2 - Math.Sqrt(2.0)) if abs(x2 - x1) < 0.0000000001 then x2 else secant x1 x2 let x0 = 1.0 let x1 = 2.0 printfn "%12.10f" (secant x0 x1) exit 0
module Fs1001 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // ヤコビの反復法 let jacobi (a:float[][]) (b:float[]) (x0:float[]) = let mutable finish:bool = false while not finish do let x1:float[] = [| 0.0;0.0;0.0;0.0 |] finish <- true for i in [0..N] do x1.[i] <- 0.0 for j in [0..N] do if j <> i then x1.[i] <- x1.[i] + a.[i].[j] * x0.[j] x1.[i] <- (b.[i] - x1.[i]) / a.[i].[i] if Math.Abs(x1.[i] - x0.[i]) > 0.0000000001 then finish <- false for i in [0..N] do x0.[i] <- x1.[i] if not finish then disp_vector x0 let a:float[][] = [| [| 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 |] |] let b:float[] = [| 20.0;16.0;8.0;17.0 |] let c:float[] = [| 0.0;0.0;0.0;0.0 |] // ヤコビの反復法 jacobi a b c printfn "X" disp_vector c exit 0
module Fs1002 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // ガウス・ザイデル法 let gauss (a:float[][]) (b:float[]) (x0:float[]) = let mutable finish:bool = false while not finish do finish <- true for i in [0..N] do let mutable x1:float = 0.0 for j in [0..N] do if j <> i then x1 <- x1 + a.[i].[j] * x0.[j] x1 <- (b.[i] - x1) / a.[i].[i] if Math.Abs(x1 - x0.[i]) > 0.0000000001 then finish <- false x0.[i] <- x1 if not finish then disp_vector x0 let a:float[][] = [| [| 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 |] |] let b:float[] = [| 20.0;16.0;8.0;17.0 |] let c:float[] = [| 0.0;0.0;0.0;0.0 |] // ガウス・ザイデル法 gauss a b c printfn "X" disp_vector c exit 0
module Fs1003 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = matrix |> Array.iter (fun row -> row |> Array.iter (fun col -> printf "%14.10f" col) printfn "" ) // 前進消去 let forward_elimination (a:float[][]) (b:float[]) = for pivot in [0..N - 1] do for row in [pivot + 1 .. N] do let s:float = a.[row].[pivot] / a.[pivot].[pivot] for col in [pivot..N] do a.[row].[col] <- a.[row].[col] - a.[pivot].[col] * s b.[row] <- b.[row] - b.[pivot] * s // 後退代入 let backward_substitution (a:float[][]) (b:float[]) = for row in (List.rev [0..N]) do for col in (List.rev [row + 1 .. N]) do b.[row] <- b.[row] - a.[row].[col] * b.[col] b.[row] <- b.[row] / a.[row].[row] // ピボット選択 let pivoting (a:float[][]) (b:float[]) = for pivot in [0..N] do // 各列で 一番値が大きい行を 探す let mutable max_row:int = pivot let mutable max_val:float = 0.0 for row in [pivot..N] do if (Math.Abs(a.[row].[pivot]) > max_val) then // 一番値が大きい行 max_val <- Math.Abs(a.[row].[pivot]) max_row <- row // 一番値が大きい行と入れ替え if (max_row <> pivot) then let mutable tmp:float = 0.0 for col in [0..N] do tmp <- a.[max_row].[col] a.[max_row].[col] <- a.[pivot].[col] a.[pivot].[col] <- tmp tmp <- b.[max_row] b.[max_row] <- b.[pivot] b.[pivot] <- tmp let mutable a:float[][] = [| [|-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 |] |] let mutable b:float[] = [| 8.0;17.0;20.0;16.0 |] // ピボット選択 pivoting a b printfn "pivoting" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 前進消去 forward_elimination a b printfn "forward elimination" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 後退代入 backward_substitution a b printfn "X" disp_vector b exit 0
module Fs1004 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = matrix |> Array.iter (fun row -> row |> Array.iter (fun col -> printf "%14.10f" col) printfn "" ) // 前進消去 let forward_elimination (a:float[][]) (b:float[]) = for pivot in [0..N] do for row in [0..N] do if row <> pivot then let s:float = a.[row].[pivot] / a.[pivot].[pivot] for col in [pivot..N] do a.[row].[col] <- a.[row].[col] - a.[pivot].[col] * s b.[row] <- b.[row] - b.[pivot] * s // 後退代入 let backward_substitution (a:float[][]) (b:float[]) = for pivot in [0..N] do b.[pivot] <- b.[pivot] / a.[pivot].[pivot] // ピボット選択 let pivoting (a:float[][]) (b:float[]) = for pivot in [0..N] do // 各列で 一番値が大きい行を 探す let mutable max_row:int = pivot let mutable max_val:float = 0.0 for row in [pivot..N] do if (Math.Abs(a.[row].[pivot]) > max_val) then // 一番値が大きい行 max_val <- Math.Abs(a.[row].[pivot]) max_row <- row // 一番値が大きい行と入れ替え if (max_row <> pivot) then let mutable tmp:float = 0.0 for col in [0..N] do tmp <- a.[max_row].[col] a.[max_row].[col] <- a.[pivot].[col] a.[pivot].[col] <- tmp tmp <- b.[max_row] b.[max_row] <- b.[pivot] b.[pivot] <- tmp let mutable a:float[][] = [| [|-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 |] |] let mutable b:float[] = [| 8.0;17.0;20.0;16.0 |] // ピボット選択 pivoting a b printfn "pivoting" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 前進消去 forward_elimination a b printfn "forward elimination" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 後退代入 backward_substitution a b printfn "X" disp_vector b exit 0
module Fs1005 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = matrix |> Array.iter (fun row -> row |> Array.iter (fun col -> printf "%14.10f" col) printfn "" ) // 前進消去 let forward_elimination (a:float[][]) (b:float[]) = for pivot in [0..N - 1] do for row in [pivot + 1 .. N] do let s:float = a.[row].[pivot] / a.[pivot].[pivot] for col in [pivot..N] do a.[row].[col] <- a.[row].[col] - a.[pivot].[col] * s // これが 上三角行列 a.[row].[pivot] <- s // これが 下三角行列 // b.[row] <- b.[row] - b.[pivot] * s // この値は変更しない // 前進代入 let forward_substitution (a:float[][]) (b:float[]) (y:float[]) = for row in [0..N] do for col in 0..(row - 1) do b.[row] <- b.[row] - a.[row].[col] * y.[col] y.[row] <- b.[row] // 後退代入 let backward_substitution (a:float[][]) (y:float[]) (x:float[]) = for row in (List.rev [0..N]) do for col in (List.rev [row + 1 .. N]) do y.[row] <- y.[row] - a.[row].[col] * x.[col] x.[row] <- y.[row] / a.[row].[row] // ピボット選択 let pivoting (a:float[][]) (b:float[]) = for pivot in [0..N] do // 各列で 一番値が大きい行を 探す let mutable max_row:int = pivot let mutable max_val:float = 0.0 for row in [pivot..N] do if (Math.Abs(a.[row].[pivot]) > max_val) then // 一番値が大きい行 max_val <- Math.Abs(a.[row].[pivot]) max_row <- row // 一番値が大きい行と入れ替え if (max_row <> pivot) then let mutable tmp:float = 0.0 for col in [0..N] do tmp <- a.[max_row].[col] a.[max_row].[col] <- a.[pivot].[col] a.[pivot].[col] <- tmp tmp <- b.[max_row] b.[max_row] <- b.[pivot] b.[pivot] <- tmp let mutable a:float[][] = [| [|-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 |] |] let mutable b:float[] = [| 8.0;17.0;20.0;16.0 |] // ピボット選択 pivoting a b printfn "pivoting" printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 前進消去 forward_elimination a b printfn "LU" disp_matrix a // Ly=b から y を求める (前進代入) let mutable y:float[] = [| 0.0;0.0;0.0;0.0 |] forward_substitution a b y printfn "Y" disp_vector y // Ux=y から x を求める (後退代入) let mutable x:float[] = [| 0.0;0.0;0.0;0.0 |] backward_substitution a y x printfn "X" disp_vector x exit 0
module Fs1006 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = matrix |> Array.iter (fun row -> row |> Array.iter (fun col -> printf "%14.10f" col) printfn "" ) // 前進消去 let forward_elimination (a:float[][]) (b:float[]) = for pivot in [0..N] do let mutable s:float = 0.0 for col in [0..(pivot - 1)] do s <- s + a.[pivot].[col] * a.[pivot].[col] // ここで根号の中が負の値になると計算できない! a.[pivot].[pivot] <- Math.Sqrt(a.[pivot].[pivot] - s) for row in [pivot + 1 .. N] do s <- 0.0 for col in [0..(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] // 前進代入 let forward_substitution (a:float[][]) (b:float[]) (y:float[]) = for row in [0..N] do for col in 0..(row - 1) do b.[row] <- b.[row] - a.[row].[col] * y.[col] y.[row] <- b.[row] / a.[row].[row] // 後退代入 let backward_substitution (a:float[][]) (y:float[]) (x:float[]) = for row in (List.rev [0..N]) do for col in (List.rev [row + 1 .. N]) do y.[row] <- y.[row] - a.[row].[col] * x.[col] x.[row] <- y.[row] / a.[row].[row] let mutable a:float[][] = [| [|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 |] |] let mutable b:float[] = [| 34.0;68.0;96.0;125.0 |] printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 前進消去 forward_elimination a b printfn "LL^T" disp_matrix a // Ly=b から y を求める (前進代入) let mutable y:float[] = [| 0.0;0.0;0.0;0.0 |] forward_substitution a b y printfn "Y" disp_vector y // L^Tx=y から x を求める (後退代入) let mutable x:float[] = [| 0.0;0.0;0.0;0.0 |] backward_substitution a y x printfn "X" disp_vector x exit 0
module Fs1007 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = matrix |> Array.iter (fun row -> row |> Array.iter (fun col -> printf "%14.10f" col) printfn "" ) // 前進消去 let forward_elimination (a:float[][]) (b:float[]) = for pivot in [0..N] do let mutable s:float = 0.0 // pivot < k の場合 for col in 0..(pivot - 1) do s <- a.[pivot].[col] for k in 0..(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] // pivot == k の場合 s <- a.[pivot].[pivot] for k in 0..(pivot - 1) do s <- s - a.[pivot].[k] * a.[pivot].[k] * a.[k].[k] a.[pivot].[pivot] <- s // 前進代入 let forward_substitution (a:float[][]) (b:float[]) (y:float[]) = for row in [0..N] do for col in 0..(row - 1) do b.[row] <- b.[row] - a.[row].[col] * y.[col] y.[row] <- b.[row] // 後退代入 let backward_substitution (a:float[][]) (y:float[]) (x:float[]) = for row in (List.rev [0..N]) do for col in (List.rev [row + 1 .. N]) do y.[row] <- y.[row] - a.[row].[col] * a.[row].[row] * x.[col] x.[row] <- y.[row] / a.[row].[row] let mutable a:float[][] = [| [|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 |] |] let mutable b:float[] = [| 34.0;68.0;96.0;125.0 |] printfn "A" disp_matrix a printfn "B" disp_vector b printfn "" // 前進消去 forward_elimination a b printfn "LDL^T" disp_matrix a // Ly=b から y を求める (前進代入) let mutable y:float[] = [| 0.0;0.0;0.0;0.0 |] forward_substitution a b y printfn "Y" disp_vector y // DL^Tx=y から x を求める (後退代入) let mutable x:float[] = [| 0.0;0.0;0.0;0.0 |] backward_substitution a y x printfn "X" disp_vector x exit 0
module Fs1101 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 正規化 (ベクトルの長さを1にする) let normarize(x:float[]):float[] = let s = x |> Array.map (fun n -> n * n) |> Array.sum x |> Array.map(fun n -> n / Math.Sqrt(s)) // ベキ乗法 let power (a:float[][]) (x:float[]) : float * float[] = let mutable x0:float[] = normarize x let mutable x1:float[] = [|0.0 ;0.0 ;0.0 ;0.0|] let mutable lambda = 0.0 let mutable e = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 100 && Math.Abs(e - (x1 |> Array.sum)) >= 0.00000000001) do if sw then sw <- false e <- x0 |> Array.sum // 1次元配列を表示 printf "%3d\t" k disp_vector x0 k <- k + 1 // 行列の積 x1 = A × x0 x1 <- [|0.0 ;0.0 ;0.0 ;0.0|] for i in [0..N] do for j in [0..N] do x1.[i] <- x1.[i] + a.[i].[j] * x0.[j] // 内積 let t1, t2 = Array.zip x0 x1 |> Array.map(fun t -> match t with | (f, s) -> (s * s, s * f)) |> Array.unzip // 固有値 lambda <- (t2 |> Array.sum) / (t1 |> Array.sum) // 正規化 (ベクトル x1 の長さを1にする) x1 <- normarize x1 x0 <- Array.copy x1 (lambda, x0) // ベキ乗法で最大固有値を求める let mutable a:float[][] = [| [|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|] |] let mutable x:float[] = [|1.0 ;0.0 ;0.0 ;0.0|] // ベキ乗法 let (lambda, x1) = power a x printfn "" printfn "eigenvalue" printfn "%14.10f" lambda printfn "eigenvector" disp_vector x1 exit 0
module Fs1102 open System let N = 3 // LU分解 let forward_elimination (a:float[][]) = for pivot in [0..N - 1] do for row in [pivot + 1 .. N] do let s:float = a.[row].[pivot] / a.[pivot].[pivot] for col in [pivot..N] do a.[row].[col] <- a.[row].[col] - a.[pivot].[col] * s // これが 上三角行列 a.[row].[pivot] <- s // これが 下三角行列 // 前進代入 let forward_substitution (a:float[][]) (y:float[]) (b:float[]) = for row in [0..N] do for col in 0..(row - 1) do b.[row] <- b.[row] - a.[row].[col] * y.[col] y.[row] <- b.[row] // 後退代入 let backward_substitution (a:float[][]) (x:float[]) (y:float[]) = for row in (List.rev [0..N]) do for col in (List.rev [row + 1 .. N]) do y.[row] <- y.[row] - a.[row].[col] * x.[col] x.[row] <- y.[row] / a.[row].[row] // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 正規化 (ベクトルの長さを1にする) let normarize(x:float[]):float[] = let s = x |> Array.map (fun n -> n * n) |> Array.sum x |> Array.map(fun n -> n / Math.Sqrt(s)) // 逆ベキ乗法 let inverse (a:float[][]) (x:float[]) : float * float[] = let mutable x0:float[] = normarize x let mutable x1:float[] = [|0.0 ;0.0 ;0.0 ;0.0|] let mutable lambda = 0.0 let mutable e = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 100 && Math.Abs(e - (x1 |> Array.sum)) >= 0.00000000001) do if sw then sw <- false e <- x0 |> Array.sum // 1次元配列を表示 printf "%3d\t" k disp_vector x0 k <- k + 1 // Ly = b から y を求める (前進代入) let mutable b = Array.copy x0 let mutable y:float[] = [|0.0 ;0.0 ;0.0 ;0.0|] forward_substitution a y b // Ux = y から x を求める (後退代入) let mutable x1 = [|0.0 ;0.0 ;0.0 ;0.0|] backward_substitution a x1 y // 内積 let t1, t2 = Array.zip x0 x1 |> Array.map(fun t -> match t with | (f, s) -> (s * s, s * f)) |> Array.unzip // 固有値 lambda <- (t2 |> Array.sum) / (t1 |> Array.sum) // 正規化 (ベクトル x1 の長さを1にする) x1 <- normarize x1 x0 <- Array.copy x1 (lambda, x0) // 逆ベキ乗法で最小固有値を求める let mutable a:float[][] = [| [|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|] |] let mutable x:float[] = [|1.0 ;0.0 ;0.0 ;0.0|] // LU分解 forward_elimination a // 逆ベキ乗法 let (lambda, x1) = inverse a x printfn "" printfn "eigenvalue" printfn "%14.10f" lambda printfn "eigenvector" disp_vector x1 exit 0
module Fs1103 open System let N = 3 // LU分解 let decomp (a:float[][]) : float[][] * float[][] = let l:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] let u:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] l.[0].[0] <- 1.0 for j in [0..N] do u.[0].[j] <- a.[0].[j] for i in [1..N] do u.[i].[0] <- 0.0 l.[0].[i] <- 0.0 l.[i].[0] <- a.[i].[0] / u.[0].[0] for i in [1..N] do l.[i].[i] <- 1.0 u.[i].[i] <- a.[i].[i] - ([| for k in 0..i -> l.[i].[k] * u.[k].[i] |] |> Array.sum) for j in [(i + 1)..N] do l.[j].[i] <- (a.[j].[i] - ([| for k in 0..i -> l.[j].[k] * u.[k].[i] |] |> Array.sum)) / u.[i].[i] u.[i].[j] <- a.[i].[j] - ([| for k in 0..i -> l.[i].[k] * u.[k].[j] |] |> Array.sum) (l, u) // 行列の積 let multiply (a:float[][]) (b:float[][]) : float[][] = [| for i in 0..N -> [| for j in 0..N -> [| for k in 0..N -> a.[i].[k] * b.[k].[j] |] |> Array.sum |] |] // 対角要素を表示 let disp_eigenvalue (a:float[][]) = for i in [0..N] do printf "%14.10f\t" a.[i].[i] printfn "" // LR分解で固有値を求める let mutable a:float[][] = [| [|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|] |] let mutable e = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 200 && e >= 0.00000000001) do if sw then sw <- false // LU分解 let (l, u) = decomp a // 行列の積 a <- multiply u l // 対角要素を表示 printf "%3d\t" k disp_eigenvalue a k <- k + 1 // 収束判定 e <- [| for i in 0..N -> [| for j in 0..(i - 1) -> Math.Abs(a.[i].[j]) |] |] |> Array.collect id |> Array.sum printfn "" printfn "eigenvalue" disp_eigenvalue a exit 0
module Fs1104 open System let N = 3 // QR分解 let decomp (a:float[][]) : float[][] * float[][] = let q:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] let r:float[][] = [| for i in 0..N -> [| for j in 0..N -> 0.0 |] |] for k in [0..N] do let x:float[] = [| for i in 0..N -> a.[i].[k] |] for j in [0..(k - 1)] do r.[j].[k] <- [| for i in 0..N -> a.[i].[k] * q.[i].[j] |] |> Array.sum r.[k].[j] <- 0.0 for i in [0..N] do x.[i] <- x.[i] - r.[j].[k] * q.[i].[j] r.[k].[k] <- Math.Sqrt(x |> Array.map(fun n -> n * n) |> Array.sum) for i in [0..N] do q.[i].[k] <- x.[i] / r.[k].[k] (q, r) // 行列の積 let multiply (a:float[][]) (b:float[][]) : float[][] = [| for i in 0..N -> [| for j in 0..N -> [| for k in 0..N -> a.[i].[k] * b.[k].[j] |] |> Array.sum |] |] // 対角要素を表示 let disp_eigenvalue (a:float[][]) = for i in [0..N] do printf "%14.10f\t" a.[i].[i] printfn "" // QR分解で固有値を求める let mutable a:float[][] = [| [|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|] |] let mutable e = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 200 && e >= 0.00000000001) do if sw then sw <- false // QR分解 let (l, u) = decomp a // 行列の積 a <- multiply u l // 対角要素を表示 printf "%3d\t" k disp_eigenvalue a k <- k + 1 // 収束判定 e <- [| for i in 0..N -> [| for j in 0..(i - 1) -> Math.Abs(a.[i].[j]) |] |] |> Array.collect id |> Array.sum printfn "" printfn "eigenvalue" disp_eigenvalue a exit 0
module Fs1105 open System let N = 3 // 1次元配列を表示 let disp_vector (row:float[]) = row |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 正規化 (ベクトルの長さを1にする) let normarize (x:float[]):float[] = let s = x |> Array.map (fun n -> n * n) |> Array.sum x |> Array.map(fun n -> n / Math.Sqrt(s)) // 対角要素を表示 let disp_eigenvalue (matrix:float[][]) = [| for i in 0..N -> matrix.[i].[i] |] |> Array.iter (fun x -> printf "%14.10f" x) printfn "" // 固有ベクトルを表示 let disp_eigenvector (matrix:float[][]) = for i in [0..N] do disp_vector ( normarize [| for j in 0..N -> matrix.[i].[j] |] ) // ヤコビ法 let jacobi (a:float[][]) (v:float[][]) : float[][] * float[][] = let mutable max_val = 0.0 let mutable k = 1 let mutable sw = true while sw || (k <= 100 && max_val >= 0.00000000001) do if sw then sw <- false // 最大値を探す let mutable p = 0 let mutable q = 0 max_val <- 0.0 for i in [0..N] do for j in [(i + 1)..N] do if max_val < Math.Abs(a.[i].[j]) then max_val <- Math.Abs(a.[i].[j]) p <- i q <- j // θ を求める let mutable t = 0.0 if Math.Abs(a.[p].[p] - a.[q].[q]) < 0.00000000001 then // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t <- Math.PI / 4.0 if a.[p].[p] < 0.0 then t <- -t else // a_{pp} ≠ a_{qq} のとき t <- Math.Atan(2.0 * a.[p].[q] / (a.[p].[p] - a.[q].[q])) / 2.0 // θ を使って 行列 U を作成し、A = U^t × A × U let c = Math.Cos(t) let s = Math.Sin(t) // U^t × A let mutable t1 = 0.0 let mutable t2 = 0.0 for i in [0..N] do t1 <- a.[p].[i] * c + a.[q].[i] * s t2 <- -a.[p].[i] * s + a.[q].[i] * c a.[p].[i] <- t1 a.[q].[i] <- t2 // 固有ベクトル t1 <- v.[p].[i] * c + v.[q].[i] * s t2 <- -v.[p].[i] * s + v.[q].[i] * c v.[p].[i] <- t1 v.[q].[i] <- t2 // A × U for i in [0..N] do 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 // 対角要素を表示 printf "%3d\t" k disp_eigenvalue a k <- k + 1 (a, v) // ヤコビ法で固有値を求める let mutable a:float[][] = [| [|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|] |] let mutable v:float[][] = [| [|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|] |] // ヤコビ法 let (a1, v1) = jacobi a v printfn "" printfn "eigenvalue" disp_eigenvalue a1 printfn "" printfn "eigenvector" disp_eigenvector v1 exit 0
module Fs1106 open System let N = 3 // 1次元配列を表示 let disp_vector (vector:float[]) = vector |> Array.iter (fun col -> printf "%14.10f" col) printfn "" // 2次元配列を表示 let disp_matrix (matrix:float[][]) = for row in [0..N] do for col in [0..N] do printf "%14.10f" matrix.[row].[col] printfn "" // ハウスホルダー変換 let tridiagonalize (a:float[][]) : float[] * float[] = let v:float[] = [| for i in 0..N -> 0.0 |] let d:float[] = [| for i in 0..N -> 0.0 |] let e:float[] = [| for i in 0..N -> 0.0 |] for k in [0..(N - 2)] do d.[k] <- a.[k].[k] let w = [| for i in [0..N] -> if (i < k + 1) then 0.0 else a.[i].[k] |] let mutable s = Math.Sqrt(w |> Array.map (fun n -> n * n) |> Array.sum) if (w.[k + 1] < 0.0) then s <- -s if (Math.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 / Math.Sqrt(w.[k + 1] * s)) for i in [(k + 1)..N] do w.[i] <- w.[i] * s for i in [(k + 1)..N] do v.[i] <- [| for j in [(k + 1)..N] -> if (j <= i) then a.[i].[j] * w.[j] else a.[j].[i] * w.[j] |] |> Array.sum s <- ( [| for i in[(k + 1)..N] -> w.[i] * v.[i] |] |> Array.sum ) / 2.0 for i in [(k + 1)..N] do v.[i] <- v.[i] - s * w.[i] for i in [(k + 1)..N] do for j in [(k + 1)..i] do a.[i].[j] <- a.[i].[j] - (w.[i] * v.[j] + w.[j] * v.[i]) // 次の行は固有ベクトルを求めないなら不要 for i in [(k + 1)..N] do a.[i].[k] <- w.[i] 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 if k < N - 1 then let w = [| for i in [0..N] -> if (i < k + 1) then 0.0 else a.[i].[k] |] for i in [(k + 1)..N] do v.[i] <- [| for j in [(k + 1)..N] -> a.[i].[j] * w.[j] |] |> Array.sum for i in [(k + 1)..N] do for j in [(k + 1)..N] do a.[i].[j] <- a.[i].[j] - v.[i] * w.[j] for i in [0..N] do a.[i].[k] <- 0.0 a.[k].[k] <- 1.0 (d, e) // QR分解 let decomp (a:float[][]) (d:float[]) (e:float[]) = e.[0] <- 1.0 let mutable h = N while (Math.Abs(e.[h]) < 0.00000000001) do h <- h - 1 while (h > 0) do e.[0] <- 0.0 let mutable l = h - 1 while (Math.Abs(e.[l]) >= 0.00000000001) do l <- l - 1 let mutable j = 1 let mutable sw = true while sw || (j <= 100 && Math.Abs(e.[h]) >= 0.00000000001) do if sw then sw <- false let mutable w = (d.[h - 1] - d.[h]) / 2.0 let mutable s = Math.Sqrt(w * w + e.[h] * e.[h]) if (w < 0.0) then s <- -s let mutable x = d.[l] - d.[h] + e.[h] * e.[h] / (w + s) let mutable y = e.[l + 1] let mutable z = 0.0 let mutable t = 0.0 let mutable u = 0.0 for k in [l..(h - 1)] do if (Math.Abs(x) >= Math.Abs(y)) then t <- -y / x u <- 1.0 / Math.Sqrt(t * t + 1.0) s <- t * u else t <- -x / y s <- 1.0 / Math.Sqrt(t * t + 1.0) * if (t < 0.0) then -1.0 else 1.0 u <- t * s 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] do x <- a.[k ].[i] y <- a.[k + 1].[i] a.[k ].[i] <- u * x - s * y a.[k + 1].[i] <- s * x + u * y if k < N - 2 then x <- e.[k + 1] y <- -s * e.[k + 2] z <- y e.[k + 2] <- u * e.[k + 2] printf "%3d\t" j j <- j + 1 disp_vector d e.[0] <- 1.0 while (Math.Abs(e.[h]) < 0.00000000001) do h <- h - 1 // 次の行は固有ベクトルを求めないなら不要 for k in [0..(N - 1)] do let mutable l = k for i in [(k + 1)..N] do if d.[i] > d.[l] then l <- i let t = d.[k] d.[k] <- d.[l] d.[l] <- t for i in [0..N] do let t = a.[k].[i] a.[k].[i] <- a.[l].[i] a.[l].[i] <- t // ハウスホルダー変換とQR分解で固有値を求める let mutable a:float[][] = [| [|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|] |] // ハウスホルダー変換 let (d, e) = tridiagonalize a // QR分解 decomp a d e printfn "" printfn "eigenvalue" disp_vector d printfn "" printfn "eigenvector" disp_matrix a exit 0