from sympy import symbols, Eq, solve, expand, diff, Matrix, Expr # Define the symbols a1, a2, a3, a4, a5, a6, a7, a8, a9, S = symbols('a1 a2 a3 a4 a5 a6 a7 a8 a9 S') a1_prime, a2_prime, a3_prime, a4_prime, a5_prime, a6_prime, a7_prime, a8_prime, a9_prime = symbols('a1_prime a2_prime a3_prime a4_prime a5_prime a6_prime a7_prime a8_prime a9_prime') # Define the equations based on the given system eq1 = Eq(a1_prime + a2_prime + a3_prime, a1 + a2 + a3 - S) eq2 = Eq(a4_prime + a5_prime + a6_prime, a4 + a5 + a6 - S) eq3 = Eq(a7_prime + a8_prime + a9_prime, a7 + a8 + a9 - S) eq4 = Eq(a1_prime + a4_prime + a7_prime, a1 + a4 + a7 - S) eq5 = Eq(a2_prime + a5_prime + a8_prime, a2 + a5 + a8 - S) eq6 = Eq(a3_prime + a6_prime + a9_prime, a3 + a6 + a9 - S) eq7 = Eq(a1_prime + a5_prime + a9_prime, a1 + a5 + a9 - S) eq8 = Eq(a3_prime + a5_prime + a7_prime, a3 + a5 + a7 - S) # Solve the system for a2' to a9' in terms of a1' # Note: Directly solving for all other variables in terms of a1' might not be straightforward due to system complexity # We may need to approach this by solving for one variable at a time or adjusting our strategy based on solvability solution = solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8), (a2_prime, a3_prime, a4_prime, a5_prime, a6_prime, a7_prime, a8_prime, a9_prime), dict=True) C = a1_prime**2 + a2_prime**2 + a3_prime**2 + a4_prime**2 + a5_prime**2 + a6_prime**2 + a7_prime**2 + a8_prime**2 + a9_prime**2 C = C.subs(solution[0]) C = expand(C) C_diff1_a1 = diff(C, a1_prime) C_diff1_a8 = diff(C, a8_prime) C_diff2_a1 = diff(C_diff1_a1, a1_prime) C_diff2_a8 = diff(C_diff1_a8, a8_prime) C_diff2_a1_a8 = diff(C_diff1_a1, a8_prime) hessian = Matrix([[C_diff2_a1, C_diff2_a1_a8], [C_diff2_a1_a8, C_diff2_a8]]) hessian.det() C_diff2_a1 C_diff1_a1, C_diff1_a8 diff_eq_1 = Eq(C_diff1_a1, 0) diff_eq_2 = Eq(C_diff1_a8, 0) critical_points = solve((diff_eq_1, diff_eq_2), (a1_prime, a8_prime)) critical_points values = {a1: 5, a2: 3, a3: 4, a4: 1, a5: 5, a6: 8, a7: 6, a8: 4, a9: 2, S: 3} solve((eq1.subs(values), eq2.subs(values)), (a1_prime, a8_prime)) values = { **values, a1_prime: 13/6, a8_prime: 5/2 } C.subs(values).evalf() prime_values = solve((eq1.subs(values), eq2.subs(values), eq3.subs(values), eq4.subs(values), eq5.subs(values), eq6.subs(values),eq7.subs(values),eq8.subs(values)), (a2_prime, a3_prime, a4_prime, a5_prime, a6_prime, a7_prime, a9_prime)) prime_values s = a3 - a3_prime + a5 - a5_prime + a7 - a7_prime s.subs({**values, **prime_values}).evalf()

Cchr