vm.fail
symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
1open Diff
2open Eval
3open Multivariate
4
5let newton_raphson expr var initial tolerance max_iter =
6 let f_prime = diff var expr in
7 let rec iterate x n =
8 if n >= max_iter then None
9 else
10 let fx = eval [(var, x)] expr in
11 if abs_float fx < tolerance then Some x
12 else
13 let fpx = eval [(var, x)] f_prime in
14 if abs_float fpx < 1e-10 then None
15 else
16 let x_next = x -. fx /. fpx in
17 if abs_float (x_next -. x) < tolerance then Some x_next
18 else iterate x_next (n + 1)
19 in
20 iterate initial 0
21
22let bisection expr var left right tolerance max_iter =
23 let rec iterate a b n =
24 if n >= max_iter then None
25 else
26 let fa = eval [(var, a)] expr in
27 let fb = eval [(var, b)] expr in
28 if fa *. fb > 0.0 then None
29 else
30 let mid = (a +. b) /. 2.0 in
31 let fmid = eval [(var, mid)] expr in
32 if abs_float fmid < tolerance || abs_float (b -. a) < tolerance then
33 Some mid
34 else if fa *. fmid < 0.0 then
35 iterate a mid (n + 1)
36 else
37 iterate mid b (n + 1)
38 in
39 iterate left right 0
40
41let trapezoidal expr var lower upper n =
42 let h = (upper -. lower) /. float_of_int n in
43 let rec sum_interior i acc =
44 if i >= n then acc
45 else
46 let x = lower +. float_of_int i *. h in
47 let fx = eval [(var, x)] expr in
48 sum_interior (i + 1) (acc +. fx)
49 in
50 let f_lower = eval [(var, lower)] expr in
51 let f_upper = eval [(var, upper)] expr in
52 let interior = sum_interior 1 0.0 in
53 h *. (f_lower /. 2.0 +. interior +. f_upper /. 2.0)
54
55let simpsons expr var lower upper n =
56 let n = if n mod 2 = 1 then n + 1 else n in
57 let h = (upper -. lower) /. float_of_int n in
58 let rec sum_terms i acc_odd acc_even =
59 if i >= n then (acc_odd, acc_even)
60 else
61 let x = lower +. float_of_int i *. h in
62 let fx = eval [(var, x)] expr in
63 if i mod 2 = 1 then
64 sum_terms (i + 1) (acc_odd +. fx) acc_even
65 else if i > 0 then
66 sum_terms (i + 1) acc_odd (acc_even +. fx)
67 else
68 sum_terms (i + 1) acc_odd acc_even
69 in
70 let f_lower = eval [(var, lower)] expr in
71 let f_upper = eval [(var, upper)] expr in
72 let (odd, even) = sum_terms 1 0.0 0.0 in
73 h /. 3.0 *. (f_lower +. 4.0 *. odd +. 2.0 *. even +. f_upper)
74
75let rec adaptive_quadrature expr var lower upper tolerance =
76 let mid = (lower +. upper) /. 2.0 in
77 let whole = simpsons expr var lower upper 10 in
78 let left_half = simpsons expr var lower mid 10 in
79 let right_half = simpsons expr var mid upper 10 in
80 let error = abs_float (whole -. (left_half +. right_half)) in
81 if error < tolerance then
82 left_half +. right_half
83 else
84 let left = adaptive_quadrature expr var lower mid (tolerance /. 2.0) in
85 let right = adaptive_quadrature expr var mid upper (tolerance /. 2.0) in
86 left +. right
87
88let gradient_descent expr vars initial learning_rate max_iter =
89 let grad_exprs = gradient vars expr in
90 let rec iterate point n =
91 if n >= max_iter then Some point
92 else
93 let env = List.combine vars point in
94 let grad_vals = List.map (eval env) grad_exprs in
95 let new_point = List.map2 (fun p g -> p -. learning_rate *. g) point grad_vals in
96 let diff = List.map2 (fun a b -> abs_float (a -. b)) new_point point in
97 let max_diff = List.fold_left max 0.0 diff in
98 if max_diff < 1e-6 then Some new_point
99 else iterate new_point (n + 1)
100 in
101 iterate initial 0
102
103let invert_matrix matrix =
104 let n = List.length matrix in
105 let augmented = List.mapi (fun i row ->
106 row @ List.init n (fun j -> if i = j then 1.0 else 0.0)
107 ) matrix in
108
109 let rec gaussian_elimination mat row =
110 if row >= n then mat
111 else
112 let pivot_row = List.nth mat row in
113 let pivot = List.nth pivot_row row in
114 if abs_float pivot < 1e-10 then mat
115 else
116 let normalized = List.map (fun x -> x /. pivot) pivot_row in
117 let updated = List.mapi (fun i r ->
118 if i = row then normalized
119 else
120 let factor = List.nth r row in
121 List.map2 (fun a b -> a -. factor *. b) r normalized
122 ) mat in
123 gaussian_elimination updated (row + 1)
124 in
125
126 let reduced = gaussian_elimination augmented 0 in
127 List.map (fun row -> List.filteri (fun i _ -> i >= n) row) reduced
128
129let newtons_method_opt expr vars initial tolerance max_iter =
130 let grad_exprs = gradient vars expr in
131 let hess_matrix = hessian vars expr in
132
133 let rec iterate point n =
134 if n >= max_iter then None
135 else
136 let env = List.combine vars point in
137 let grad_vals = List.map (eval env) grad_exprs in
138 let hess_vals = List.map (fun row ->
139 List.map (eval env) row
140 ) hess_matrix in
141
142 let hess_inv = invert_matrix hess_vals in
143 let delta = List.map (fun row ->
144 List.fold_left2 (fun acc h g -> acc +. h *. g) 0.0 row grad_vals
145 ) hess_inv in
146
147 let new_point = List.map2 (fun p d -> p -. d) point delta in
148 let diff = List.map2 (fun a b -> abs_float (a -. b)) new_point point in
149 let max_diff = List.fold_left max 0.0 diff in
150
151 if max_diff < tolerance then Some new_point
152 else iterate new_point (n + 1)
153 in
154 iterate initial 0