I need to calculate the vertical natural frequency and the horizontal natural frequency of a beam. In order to do this I have been given two values of the vertical and horizontal natural frequency which I need to approximate as closey as possible. Now in order to do this I have to use the moment of inertia in the y and z direction and the surface area. So the values that I need to calculate in this instance are the width(B), the height(H) and the thickness(t_h and t_v). For this I have used two solvers. The first one computes the values for B and H and the second one computes the values for t_h and t_v. Now the first solver works just fine and gives me numerical values. However, the second one somehow does not and only gives me []. Below I have put the code in that I have used. I cannot seem to understand what I am doing wrong in this code and I'm hoping that somebody can help me with his, because I have been at it for a couple of days and cannot seem to figure it out. import numpy as np from sympy import symbols, Eq, solve w_2_vb_g=4.341 w_2_hb_g=2.601 E=210E9 rho=7700 L=100 kL_2=7.85 kL_3=11.00 B,H,t_h,t_v=symbols('B H t_h t_v', real=True, positive=True) w_2_vb_rec=Eq(((((kL_2)**2)/(L**2))*(((E*((1/12)*H*B**3))/(rho*(B*H)))**(1/2))), w_2_vb_g) w_2_hb_rec=Eq(((((kL_2)**2)/(L**2))*(((E*((1/12)*B*H**3))/(rho*(H*B)))**(1/2))), w_2_hb_g) Solution=solve((w_2_hb_rec,w_2_vb_rec), (B,H)) print(Solution) print("Solutions:") for var, val in Solution.items(): print(f"{var}: {val}") list_from_dict = list(Solution.items()) print("List from dict:", list_from_dict) B=round(float(list_from_dict[1][1]*1.03),3) H=round(float(list_from_dict[0][1]*1.03),3) B_list=np.array([t_h, (H-2*t_h), B, t_v]) H_list=np.array([B, t_v, t_h, (H-2*t_h)]) D0_list=np.array([0, (0.5*B-0.5*t_v), (0.5*H-0.5*t_h),0]) w_2_vb=Eq(((((kL_2)**2)/(L**2))*(((E*(2*(((1/12)*B_list[0]*H_list[0]**3)+B_list[0]*H_list[0]*D0_list[0]**2)+2*(((1/12)*B_list[1]*H_list[1]**3)+B_list[1]*H_list[1]*D0_list[1]**2))/(rho*(((2*(B_list[0]*H_list[0])))+(2*(B_list[1]*H_list[1]))))**(1/2))))), w_2_vb_g) w_2_hb=Eq(((((kL_2)**2)/(L**2))*(((E*(2*(((1/12)*B_list[2]*H_list[2]**3)+B_list[2]*H_list[2]*D0_list[2]**2)+2*(((1/12)*B_list[3]*H_list[3]**3)+B_list[3]*H_list[3]*D0_list[3]**2))/(rho*(((2*(H_list[2]*B_list[2])))+(2*(H_list[3]*H_list[3]))))**(1/2))))), w_2_hb_g) Solution_beam=solve((w_2_hb,w_2_vb), (t_h,t_v)) print(Solution_beam) Below You can see the output that my console gives: {H: 0.279980237079456, B: 0.467279588297547} Solutions: H: 0.279980237079456 B: 0.467279588297547 List from dict: [(H, 0.279980237079456), (B, 0.467279588297547)] [] I have already tried different solver such as linsolve() or solve_undetermined_coeffs but neither has worked. This could be attributed to my code and that it doesn't work in this case but as far as I can see it should. I have also considered by splitting the equations in multiple parts so as to make it more comprehensive for the solver. However, this is not possible because I have not been given a starting value for the moment of inertia or the surface area which are the only two options in which it could've worked. Hopefully somebody can help, many thanks in advance! Continue reading...