system of non-linear algebraic equations not converging with SymPy 'nsolve'

Question

I have written the following code for finding the solution of system of three non-linear algebraic equations (Nonlinear_eq_1, Nonlinear_eq_2 and Nonlinear_eq_3).

import sympy as smp
import scipy as sp
import numpy as np

pi = np.pi

# Problem-specific given data
Number_of_shape_functions_used = 3

g_0 = 5.0e-06                           # Initial gap m
h_beam = 1.80e-06                       # Beam height m
b_beam_prismatic = 15.0e-06             # Beam width m
L_beam = 100.0e-06                      # Beam length m

epsilon_p = 8.8541878128e-12             # Permittivity of air in Farads/meter

nu = 0.23                                # Poisson's ratio 
E_Young = 166.0e9                        # Young's modulus Pa
E_Young_effective = (E_Young) / (1 - ((nu) * (nu)))              # Define iff((b / h) >= 5) i.e., ((E_Young) / (1 - ((nu) * (nu))))

# Problem-specific derived quantities
MI_prismatic = ((b_beam_prismatic * (h_beam * h_beam * h_beam)) / 12)
b_bar = (b_beam_prismatic / (0.65 * g_0))

# Defining all symbols used
x, a1, a2, a3, V = smp.symbols('x a1 a2 a3 V')

K1 = (np.sinh(1.875104069) + np.sin(1.875104069))/(np.cosh(1.875104069) + np.cos(1.875104069))
phi1 = (smp.sin(1.875104069*x) - smp.sinh(1.875104069*x)) - (K1*(smp.cos(1.875104069*x) - smp.cosh(1.875104069*x)))

K2 = (np.sinh(4.694091133) + np.sin(4.694091133))/(np.cosh(4.694091133) + np.cos(4.694091133))
phi2 = (smp.sin(4.694091133*x) - smp.sinh(4.694091133*x)) - (K2*(smp.cos(4.694091133*x) - smp.cosh(4.694091133*x)))

K3 = (np.sinh(7.854757438) + np.sin(7.854757438))/(np.cosh(7.854757438) + np.cos(7.854757438))
phi3 = (smp.sin(7.854757438*x) - smp.sinh(7.854757438*x)) - (K3*(smp.cos(7.854757438*x) - smp.cosh(7.854757438*x)))
phi = [0, 2 * phi1 / phi1.subs(x, 1), 2 * phi2 / phi2.subs(x, 1), 2 * phi3 / phi3.subs(x, 1)]

# Derivatives involved in the residue
D1 = smp.diff(phi[1], x, 4)
D2 = smp.diff(phi[2], x, 4)
D3 = smp.diff(phi[3], x, 4)

# Arrays for storing MROM coefficients
g1 = np.zeros([4, 4], dtype=float)
g2 = np.zeros([4, 10], dtype=float)
g3 = np.zeros([4, 19], dtype=float)
g4 = np.zeros([4, 31], dtype=float)
g5 = np.zeros([4, 2], dtype=float)
g6 = np.zeros([4, 2], dtype=float)
g7 = np.zeros([4, 4], dtype=float)
g8 = np.zeros([4, 7], dtype=float)

# Coefficients g1_m
for p in range(1, Number_of_shape_functions_used + 1):
    g1[p][1] = (sp.integrate.quad(smp.lambdify(x, (phi[p] * D1)), 0, 1))[0]
    g1[p][2] = (sp.integrate.quad(smp.lambdify(x, (phi[p] * D2)), 0, 1))[0]
    g1[p][3] = (sp.integrate.quad(smp.lambdify(x, (phi[p] * D3)), 0, 1))[0]

# Coefficients g2_im
for p in range(1, Number_of_shape_functions_used + 1):
    g2[p][1] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[1] * D1)), 0, 1))[0]
    g2[p][2] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[1] * D2)), 0, 1))[0]
    g2[p][3] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[1] * D3)), 0, 1))[0]

    g2[p][4] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[2] * D1)), 0, 1))[0]
    g2[p][5] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[2] * D2)), 0, 1))[0]
    g2[p][6] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[2] * D3)), 0, 1))[0]

    g2[p][7] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[3] * D1)), 0, 1))[0]
    g2[p][8] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[3] * D2)), 0, 1))[0]
    g2[p][9] = (sp.integrate.quad(smp.lambdify(x, (-3 * phi[p] * phi[3] * D3)), 0, 1))[0]

# Coefficients g3_ijm
for p in range(1, Number_of_shape_functions_used + 1):
    g3[p][1] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[1] * D1)), 0, 1))[0]
    g3[p][2] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[1] * D2)), 0, 1))[0]
    g3[p][3] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[1] * D3)), 0, 1))[0]

    g3[p][4] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[2] * D1)), 0, 1))[0]
    g3[p][5] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[2] * D2)), 0, 1))[0]
    g3[p][6] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[2] * D3)), 0, 1))[0]

    g3[p][7] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[3] * D1)), 0, 1))[0]
    g3[p][8] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[3] * D2)), 0, 1))[0]
    g3[p][9] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[1] * phi[3] * D3)), 0, 1))[0]

    g3[p][10] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[2] * phi[2] * D1)), 0, 1))[0]
    g3[p][11] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[2] * phi[2] * D2)), 0, 1))[0]
    g3[p][12] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[2] * phi[2] * D3)), 0, 1))[0]

    g3[p][13] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[2] * phi[3] * D1)), 0, 1))[0]
    g3[p][14] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[2] * phi[3] * D2)), 0, 1))[0]
    g3[p][15] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[2] * phi[3] * D3)), 0, 1))[0]

    g3[p][16] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[3] * phi[3] * D1)), 0, 1))[0]
    g3[p][17] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[3] * phi[3] * D2)), 0, 1))[0]
    g3[p][18] = (sp.integrate.quad(smp.lambdify(x, (3 * phi[p] * phi[3] * phi[3] * D3)), 0, 1))[0]

# Coefficients g4_ijkm
for p in range(1, Number_of_shape_functions_used + 1):
    g4[p][1] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[1] * D1)), 0, 1))[0]
    g4[p][2] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[1] * D2)), 0, 1))[0]
    g4[p][3] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[1] * D3)), 0, 1))[0]

    g4[p][4] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[2] * D1)), 0, 1))[0]
    g4[p][5] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[2] * D2)), 0, 1))[0]
    g4[p][6] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[2] * D3)), 0, 1))[0]

    g4[p][7] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[3] * D1)), 0, 1))[0]
    g4[p][8] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[3] * D2)), 0, 1))[0]
    g4[p][9] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] * phi[3] * D3)), 0, 1))[0]

    g4[p][10] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[1] * D1)), 0, 1))[0]
    g4[p][11] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[1] * D2)), 0, 1))[0]
    g4[p][12] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[1] * D3)), 0, 1))[0]

    g4[p][13] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[2] * D1)), 0, 1))[0]
    g4[p][14] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[2] * D2)), 0, 1))[0]
    g4[p][15] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[2] * D3)), 0, 1))[0]

    g4[p][16] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[3] * D1)), 0, 1))[0]
    g4[p][17] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[3] * D2)), 0, 1))[0]
    g4[p][18] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] * phi[3] * D3)), 0, 1))[0]

    g4[p][19] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[1] * D1)), 0, 1))[0]
    g4[p][20] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[1] * D2)), 0, 1))[0]
    g4[p][21] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[1] * D3)), 0, 1))[0]

    g4[p][22] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[2] * D1)), 0, 1))[0]
    g4[p][23] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[2] * D2)), 0, 1))[0]
    g4[p][24] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[2] * D3)), 0, 1))[0]

    g4[p][25] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[3] * D1)), 0, 1))[0]
    g4[p][26] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[3] * D2)), 0, 1))[0]
    g4[p][27] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] * phi[3] * D3)), 0, 1))[0]

    g4[p][28] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[2] * phi[3] * D1)), 0, 1))[0]
    g4[p][29] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[2] * phi[3] * D2)), 0, 1))[0]
    g4[p][30] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[2] * phi[3] * D3)), 0, 1))[0]

# Coefficients g5
for p in range(1, Number_of_shape_functions_used + 1):
    g5[p][1] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * (1 + 1 / b_bar))), 0, 1))[0]

# Coefficients g6
for p in range(1, Number_of_shape_functions_used + 1):
    g6[p][1] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * 1)), 0, 1))[0]

# Coefficients g7_i
for p in range(1, Number_of_shape_functions_used + 1):
    g7[p][1] = (sp.integrate.quad(smp.lambdify(x, (phi[p] * phi[1])), 0, 1))[0]
    g7[p][2] = (sp.integrate.quad(smp.lambdify(x, (phi[p] * phi[2])), 0, 1))[0]
    g7[p][3] = (sp.integrate.quad(smp.lambdify(x, (phi[p] * phi[3])), 0, 1))[0]

# Coefficients g8_ij
for p in range(1, Number_of_shape_functions_used + 1):
    g8[p][1] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[1] / b_bar)), 0, 1))[0]
    g8[p][2] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[2] / b_bar)), 0, 1))[0]
    g8[p][3] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[1] * phi[3] / b_bar)), 0, 1))[0]

    g8[p][4] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[2] / b_bar)), 0, 1))[0]
    g8[p][5] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[2] * phi[3] / b_bar)), 0, 1))[0]
    g8[p][6] = (sp.integrate.quad(smp.lambdify(x, (-1 * phi[p] * phi[3] * phi[3] / b_bar)), 0, 1))[0]

# Linear algebraic equations
Linear_eq_1 = smp.simplify(a1 * g1[1][1] + a2 * g1[1][2] + a3 * g1[1][3] + V * V * g5[1][1] + \
    a1 * V * V * g7[1][1] + a2 * V * V * g7[1][2] + a3 * V * V * g7[1][3] + 2 * V * V * a1/b_bar)

Linear_eq_2 = smp.simplify(a1 * g1[2][1] + a2 * g1[2][2] + a3 * g1[2][3] + V * V * g5[2][1] + \
    a1 * V * V * g7[2][1] + a2 * V * V * g7[2][2] + a3 * V * V * g7[2][3] + 2 * V * V * a2/b_bar)

Linear_eq_3 = smp.simplify(a1 * g1[3][1] + a2 * g1[3][2] + a3 * g1[3][3] + V * V * g5[3][1] + \
    a1 * V * V * g7[3][1] + a2 * V * V * g7[3][2] + a3 * V * V * g7[3][3] + 2 * V * V * a3/b_bar)

# Non-linear algebraic equations
Nonlinear_eq_1 = smp.simplify(a1 * g1[1][1] + a2 * g1[1][2] + a3 * g1[1][3] + a1 * a1 * g2[1][1] + a1 * a2 * g2[1][2] + \
    a1 * a3 * g2[1][3] + a2 * a1 * g2[1][4] + a2 * a2 * g2[1][5] + a2 * a3 * g2[1][6] + a3 * a1 * g2[1][7] + a3 * a2 * g2[1][8] + \
    a3 * a3 * g2[1][9] + a1 * a1 * a1 * g3[1][1] + a1 * a1 * a2 * g3[1][2] + a1 * a1 * a3 * g3[1][3] + 2 * a1 * a2 * a1 * g3[1][4] + \
    2 * a1 * a2 * a2 * g3[1][5] + 2 * a1 * a2 * a3 * g3[1][6] + 2 * a1 * a3 * a1 * g3[1][7] + 2 * a1 * a3 * a2 * g3[1][8] + \
    2 * a1 * a3 * a3 * g3[1][9] + a2 * a2 * a1 * g3[1][10] + a2 * a2 * a2 * g3[1][11] + a2 * a2 * a3 * g3[1][12] + \
    2 * a2 * a3 * a1 * g3[1][13] + 2 * a2 * a3 * a2 * g3[1][14] + 2 * a2 * a3 * a3 * g3[1][15] + a3 * a3 * a1 * g3[1][16] + \
    a3 * a3 * a2 * g3[1][17] + a3 * a3 * a3 * g3[1][18] + a1 * a1 * a1 * a1 * g4[1][1] + a1 * a1 * a1 * a2 * g4[1][2] + \
    a1 * a1 * a1 * a3 * g4[1][3] + 3 * a1 * a1 * a2 * a1 * g4[1][4] + 3 * a1 * a1 * a2 * a2 * g4[1][5] + \
    3 * a1 * a1 * a2 * a3 * g4[1][6] + 3 * a1 * a1 * a3 * a1 * g4[1][7] + 3 * a1 * a1 * a3 * a2 * g4[1][8] + \
    3 * a1 * a1 * a3 * a3 * g4[1][9] + 3 * a2 * a2 * a1 * a1 * g4[1][10] + 3 * a2 * a2 * a1 * a2 * g4[1][11] + \
    3 * a2 * a2 * a1 * a3 * g4[1][12] + a2 * a2 * a2 * a1 * g4[1][13] + a2 * a2 * a2 * a2 * g4[1][14] + \
    a2 * a2 * a2 * a3 * g4[1][15] + 3 * a2 * a2 * a3 * a1 * g4[1][16] + 3 * a2 * a2 * a3 * a2 * g4[1][17] + \
    3 * a2 * a2 * a3 * a3 * g4[1][18] + 3 * a3 * a3 * a1 * a1 * g4[1][19] + 3 * a3 * a3 * a1 * a2 * g4[1][20] + \
    3 * a3 * a3 * a1 * a3 * g4[1][21] + 3 * a3 * a3 * a2 * a1 * g4[1][22] + 3 * a3 * a3 * a2 * a2 * g4[1][23] + \
    3 * a3 * a3 * a2 * a3 * g4[1][24] + a3 * a3 * a3 * a1 * g4[1][25] + a3 * a3 * a3 * a2 * g4[1][26] + \
    a3 * a3 * a3 * a3 * g4[1][27] + 6 * a1 * a2 * a3 * a1 * g4[1][28] + 6 * a1 * a2 * a3 * a2 * g4[1][29] + \
    6 * a1 * a2 * a3 * a3 * g4[1][30] + V * V * g5[1][1] + a1 * V * V * g7[1][1] + a2 * V * V * g7[1][2] + \
    a3 * V * V * g7[1][3] + 2 * V * V * a1/b_bar + a1 * a1 * V * V * g8[1][1] + 2 * a1 * a2 * V * V * g8[1][2] + \
    2 * a1 * a3 * V * V * g8[1][3] + a2 * a2 * V * V * g8[1][4] + 2 * a2 * a3 * V * V * g8[1][5] + a3 * a3 * V * V * g8[1][6])
                 
Nonlinear_eq_2 = smp.simplify(a1 * g1[2][1] + a2 * g1[2][2] + a3 * g1[2][3] + a1 * a1 * g2[2][1] + a1 * a2 * g2[2][2] + \
    a1 * a3 * g2[2][3] + a2 * a1 * g2[2][4] + a2 * a2 * g2[2][5] + a2 * a3 * g2[2][6] + a3 * a1 * g2[2][7] + a3 * a2 * g2[2][8] + \
    a3 * a3 * g2[2][9] + a1 * a1 * a1 * g3[2][1] + a1 * a1 * a2 * g3[2][2] + a1 * a1 * a3 * g3[2][3] + 2 * a1 * a2 * a1 * g3[2][4] + \
    2 * a1 * a2 * a2 * g3[2][5] + 2 * a1 * a2 * a3 * g3[2][6] + 2 * a1 * a3 * a1 * g3[2][7] + 2 * a1 * a3 * a2 * g3[2][8] + \
    2 * a1 * a3 * a3 * g3[2][9] + a2 * a2 * a1 * g3[2][10] + a2 * a2 * a2 * g3[2][11] + a2 * a2 * a3 * g3[2][12] + \
    2 * a2 * a3 * a1 * g3[2][13] + 2 * a2 * a3 * a2 * g3[2][14] + 2 * a2 * a3 * a3 * g3[2][15] + a3 * a3 * a1 * g3[2][16] + \
    a3 * a3 * a2 * g3[2][17] + a3 * a3 * a3 * g3[2][18] + a1 * a1 * a1 * a1 * g4[2][1] + a1 * a1 * a1 * a2 * g4[2][2] + \
    a1 * a1 * a1 * a3 * g4[2][3] + 3 * a1 * a1 * a2 * a1 * g4[2][4] + 3 * a1 * a1 * a2 * a2 * g4[2][5] + \
    3 * a1 * a1 * a2 * a3 * g4[2][6] + 3 * a1 * a1 * a3 * a1 * g4[2][7] + 3 * a1 * a1 * a3 * a2 * g4[2][8] + \
    3 * a1 * a1 * a3 * a3 * g4[2][9] + 3 * a2 * a2 * a1 * a1 * g4[2][10] + 3 * a2 * a2 * a1 * a2 * g4[2][11] + \
    3 * a2 * a2 * a1 * a3 * g4[2][12] + a2 * a2 * a2 * a1 * g4[2][13] + a2 * a2 * a2 * a2 * g4[2][14] + \
    a2 * a2 * a2 * a3 * g4[2][15] + 3 * a2 * a2 * a3 * a1 * g4[2][16] + 3 * a2 * a2 * a3 * a2 * g4[2][17] + \
    3 * a2 * a2 * a3 * a3 * g4[2][18] + 3 * a3 * a3 * a1 * a1 * g4[2][19] + 3 * a3 * a3 * a1 * a2 * g4[2][20] + \
    3 * a3 * a3 * a1 * a3 * g4[2][21] + 3 * a3 * a3 * a2 * a1 * g4[2][22] + 3 * a3 * a3 * a2 * a2 * g4[2][23] + \
    3 * a3 * a3 * a2 * a3 * g4[2][24] + a3 * a3 * a3 * a1 * g4[2][25] + a3 * a3 * a3 * a2 * g4[2][26] + \
    a3 * a3 * a3 * a3 * g4[2][27] + 6 * a1 * a2 * a3 * a1 * g4[2][28] + 6 * a1 * a2 * a3 * a2 * g4[2][29] + \
    6 * a1 * a2 * a3 * a3 * g4[2][30] + V * V * g5[2][1] + a1 * V * V * g7[2][1] + a2 * V * V * g7[2][2] + \
    a3 * V * V * g7[2][3] + 2 * V * V * a2/b_bar + a1 * a1 * V * V * g8[2][1] + 2 * a1 * a2 * V * V * g8[2][2] + \
    2 * a1 * a3 * V * V * g8[2][3] + a2 * a2 * V * V * g8[2][4] + 2 * a2 * a3 * V * V * g8[2][5] + a3 * a3 * V * V * g8[2][6])

Nonlinear_eq_3 = smp.simplify(a1 * g1[3][1] + a2 * g1[3][2] + a3 * g1[3][3] + a1 * a1 * g2[3][1] + a1 * a2 * g2[3][2] + \
    a1 * a3 * g2[3][3] + a2 * a1 * g2[3][4] + a2 * a2 * g2[3][5] + a2 * a3 * g2[3][6] + a3 * a1 * g2[3][7] + a3 * a2 * g2[3][8] + \
    a3 * a3 * g2[3][9] + a1 * a1 * a1 * g3[3][1] + a1 * a1 * a2 * g3[3][2] + a1 * a1 * a3 * g3[3][3] + 2 * a1 * a2 * a1 * g3[3][4] + \
    2 * a1 * a2 * a2 * g3[3][5] + 2 * a1 * a2 * a3 * g3[3][6] + 2 * a1 * a3 * a1 * g3[3][7] + 2 * a1 * a3 * a2 * g3[3][8] + \
    2 * a1 * a3 * a3 * g3[3][9] + a2 * a2 * a1 * g3[3][10] + a2 * a2 * a2 * g3[3][11] + a2 * a2 * a3 * g3[3][12] + \
    2 * a2 * a3 * a1 * g3[3][13] + 2 * a2 * a3 * a2 * g3[3][14] + 2 * a2 * a3 * a3 * g3[3][15] + a3 * a3 * a1 * g3[3][16] + \
    a3 * a3 * a2 * g3[3][17] + a3 * a3 * a3 * g3[3][18] + a1 * a1 * a1 * a1 * g4[3][1] + a1 * a1 * a1 * a2 * g4[3][2] + \
    a1 * a1 * a1 * a3 * g4[3][3] + 3 * a1 * a1 * a2 * a1 * g4[3][4] + 3 * a1 * a1 * a2 * a2 * g4[3][5] + \
    3 * a1 * a1 * a2 * a3 * g4[3][6] + 3 * a1 * a1 * a3 * a1 * g4[3][7] + 3 * a1 * a1 * a3 * a2 * g4[3][8] + \
    3 * a1 * a1 * a3 * a3 * g4[3][9] + 3 * a2 * a2 * a1 * a1 * g4[3][10] + 3 * a2 * a2 * a1 * a2 * g4[3][11] + \
    3 * a2 * a2 * a1 * a3 * g4[3][12] + a2 * a2 * a2 * a1 * g4[3][13] + a2 * a2 * a2 * a2 * g4[3][14] + \
    a2 * a2 * a2 * a3 * g4[3][15] + 3 * a2 * a2 * a3 * a1 * g4[3][16] + 3 * a2 * a2 * a3 * a2 * g4[3][17] + \
    3 * a2 * a2 * a3 * a3 * g4[3][18] + 3 * a3 * a3 * a1 * a1 * g4[3][19] + 3 * a3 * a3 * a1 * a2 * g4[3][20] + \
    3 * a3 * a3 * a1 * a3 * g4[3][21] + 3 * a3 * a3 * a2 * a1 * g4[3][22] + 3 * a3 * a3 * a2 * a2 * g4[3][23] + \
    3 * a3 * a3 * a2 * a3 * g4[3][24] + a3 * a3 * a3 * a1 * g4[3][25] + a3 * a3 * a3 * a2 * g4[3][26] + \
    a3 * a3 * a3 * a3 * g4[3][27] + 6 * a1 * a2 * a3 * a1 * g4[3][28] + 6 * a1 * a2 * a3 * a2 * g4[3][29] + \
    6 * a1 * a2 * a3 * a3 * g4[3][30] + V * V * g5[3][1] + a1 * V * V * g7[3][1] + a2 * V * V * g7[3][2] + \
    a3 * V * V * g7[3][3] + 2 * V * V * a3/b_bar + a1 * a1 * V * V * g8[3][1] + 2 * a1 * a2 * V * V * g8[3][2] + \
    2 * a1 * a3 * V * V * g8[3][3] + a2 * a2 * V * V * g8[3][4] + 2 * a2 * a3 * V * V * g8[3][5] + a3 * a3 * V * V * g8[3][6])

{this code gives the following equations which are generated directly below for convenience}

from sympy import *
from sympy.abc import V
A = var('a1:4')
NL=Tuple(-0.320195959022147*V**2*a1**2 - 0.172708561630121*V**2*a1*a2 - 0.0312953011133775*V**2*a1*a3 + 1.43333333337127*V**2*a1 - 0.206766860454431*V**2*a2**2 - 0.195889019764585*V**2*a2*a3 + 6.46442621654586e-12*V**2*a2 - 0.180494407384852*V**2*a3**2 - 1.21986386270034e-10*V**2*a3 - 0.952639969874236*V**2 - 48.2204330265944*a1**4 - 1082.06396855855*a1**3*a2 - 3775.39653054958*a1**3*a3 + 87.1048397690087*a1**3 - 3088.21338141745*a1**2*a2**2 - 20421.7770534064*a1**2*a2*a3 + 1537.59803291102*a1**2*a2 - 20039.9749965823*a1**2*a3**2 + 4109.9738798034*a1**2*a3 - 54.8083218036971*a1**2 - 1649.72664616273*a1*a2**3 - 15710.621320911*a1*a2**2*a3 + 3788.30585004389*a1*a2**2 - 19165.1563796871*a1*a2*a3**2 + 22754.1216551676*a1*a2*a3 - 595.303911613836*a1*a2 - 5869.02732603284*a1*a3**3 + 25405.9129186628*a1*a3**2 - 827.403891746658*a1*a3 + 12.3623633763912*a1 - 947.881908357562*a2**4 - 5711.63059308329*a2**3*a3 + 380.116996933847*a2**3 - 17156.2234266921*a2**2*a3**2 + 8614.34353685311*a2**2*a3 - 1390.00433258433*a2**2 - 11164.1097134021*a2*a3**3 + 4645.59003446727*a2*a3**2 - 5820.70446108739*a2*a3 + 3.16025960955812e-9*a2 - 5737.95738495622*a3**4 + 2806.71958369875*a3**3 - 9513.14278780819*a3**2 - 4.69857241114369e-7*a3, -0.0863542808150603*V**2*a1**2 - 0.41353372090882*V**2*a1*a2 - 0.195889019764585*V**2*a1*a3 + 6.46442621654586e-12*V**2*a1 + 0.0984182662400039*V**2*a2**2 - 0.245492930006062*V**2*a2*a3 + 1.43333333329724*V**2*a2 + 0.0874525904012477*V**2*a3**2 - 1.41143041698655e-10*V**2*a3 + 0.527955339021191*V**2 - 25.5964728945001*a1**4 - 1080.52457693316*a1**3*a2 - 6061.12619989492*a1**3*a3 + 37.2534758840929*a1**3 - 1677.49471613007*a1**2*a2**2 - 14158.7079007895*a1**2*a2*a3 + 1965.58735680544*a1**2*a2 - 9023.57517892548*a1**2*a3**2 + 10126.4578527536*a1**2*a3 - 14.7813645946492*a1**2 - 2867.78085669168*a1*a2**3 - 15594.4535467757*a1*a2**2*a3 + 769.912597995959*a1*a2**2 - 32421.1479796308*a1*a2*a3**2 + 15522.4136375446*a1*a2*a3 - 1425.39686114951*a1*a2 - 10720.4053335747*a1*a3**3 + 4374.17425494042*a1*a3**2 - 5179.03108589664*a1*a3 + 7.99098565096301e-11*a1 + 158.90983790103*a2**4 - 7135.71086984839*a2**3*a3 + 2593.58482445554*a2**3 - 1338.21512086288*a2**2*a3**2 + 3879.11895709259*a2**2*a3 + 661.623512483564*a2**2 - 14929.6351109554*a2*a3**3 + 29348.7783732883*a2*a3**2 - 7294.64976939622*a2*a3 + 485.518818506518*a2 + 1489.34391821515*a3**4 + 642.195069912758*a3**3 + 4609.27843528281*a3**2 - 5.40270065130244e-7*a3, -0.0156476505566887*V**2*a1**2 - 0.195889019764585*V**2*a1*a2 - 0.360988814770438*V**2*a1*a3 - 1.21986386270034e-10*V**2*a1 - 0.122746465003048*V**2*a2**2 + 0.174905180801629*V**2*a2*a3 - 1.41143041698655e-10*V**2*a2 - 0.117555443821002*V**2*a3**2 + 1.43333333371768*V**2*a3 - 0.309550777862571*V**2 - 12.1428917577821*a1**4 - 824.110957453294*a1**3*a2 - 6723.23987716781*a1**3*a3 + 13.2616540737603*a1**3 - 1633.00844658897*a1**2*a2**2 - 10200.0024546644*a1**2*a2*a3 + 1348.62810950915*a1**2*a2 - 5881.72065030481*a1**2*a3**2 + 12764.7384379054*a1**2*a3 - 2.67842689138796*a1**2 - 1594.10182719871*a1*a2**3 - 19146.3518680383*a1*a2**2*a3 + 1773.144252372*a1*a2**2 - 22818.2952318094*a1*a2*a3**2 + 4938.2800063943*a1*a2*a3 - 675.203930870838*a1*a2 - 17232.5070825099*a1*a3**3 + 5622.5544349222*a1*a3**2 - 9544.03822943416*a1*a3 - 1.50804213561173e-9*a1 - 658.266063611188*a2**4 - 546.991056446739*a2**3*a3 + 394.212909957935*a2**3 - 16147.3625924508*a2**2*a3**2 + 17313.6165030967*a2**2*a3 - 825.171489223677*a2**2 + 4657.99516216124*a2*a3**3 + 1366.30108217397*a2*a3**2 + 5197.1844434713*a2*a3 - 6.85251144716403e-8*a2 - 5178.70585826686*a3**4 + 19331.7005096786*a3**3 - 6195.88018684571*a3**2 + 3806.54626739517*a3)

def linearize(eq, v):
    l = eq.replace(
        lambda x: x.is_Pow and x.base in v and x.exp > 1,
        lambda x: 0)
    return l.replace(
        lambda x: x.is_Mul and sum(
            i.as_base_exp()[1] for i in x.args if
            i.as_base_exp()[0] in v)>1, 
        lambda x: 0)

LIN = Tuple(*[linearize(i, A) for i in NL])

An attempt is made to follow a solution as the voltage is increased:

Voltage_initial = 0.000001 
Voltage_increment = 0.1 
Voltage_accuracy = 0.0000001 

while Voltage_increment > Voltage_accuracy: 
    # Solving the linear algebraic equations
    lsol = smp.solve(LIN.subs(V, Voltage_initial), (a1, a2, a3))
    linear_answers = [lsol[i] for i in A]
    
    # Solving the non-linear algebraic equations
    try:
        ans_nonlinear = smp.nsolve(NL.subs(V, Voltage_initial), A, linear_answers)
        print(Voltage_initial)
        Voltage_initial += Voltage_increment
    except ValueError:
        print('Voltage value is more than the pull-in voltage') 
        Voltage_initial = Voltage_initial - Voltage_increment 
        Voltage_increment = Voltage_increment / 10 

The intent is to use the final nonlinear answer to compute the pull-in values:

amp_phi_1 = float(ans_nonlinear[0]) 
amp_phi_2 = float(ans_nonlinear[1]) 
amp_phi_3 = float(ans_nonlinear[2]) 

#Pull-in parameters 
pull_in_voltage = ((Voltage_initial) / ((((epsilon_p) * (b_beam_prismatic) * ((L_beam)**(4))) / 
((2) * (E_Young_effective) * (MI_prismatic) * ((g_0)**(3))))**(0.5))) 

pull_in_displacement = (-1 * g_0) * ((smp.simplify(amp_phi_1 * phi[1] + amp_phi_2 * phi[2] + amp_phi_3 * phi[3])).subs(x, 1)) 

print(pull_in_voltage)

print(pull_in_displacement)

The system is stable below the critical voltage value and it becomes unstable above the critical voltage value.

This traslates in terms of the code as follows:

• When V < critical voltage, the system of non-linear algebraic equations gives solution using nsolve.
• When V exceeds the critical voltage, the system of non-linear algebraic equations gives 'ValueError' using nsolve.

I have used this logic along with try-except to find the critical voltage value in an iterative manner.

Now the strange behaviour I am observing with regard to this code is as follows:

• For Voltage_initial = 0.1 or 0.001 or 0.000001, the code does not give converged answer for g_0 = 5.0e-06, h_beam = 1.80e-06, b_beam_prismatic = 15.0e-06, L_beam = 100.0e-06
• However, for Voltage_initial = 0.01 or 0.0001 or 0.00001, the same code gives converged answer for g_0 = 5.0e-06, h_beam = 1.80e-06, b_beam_prismatic = 15.0e-06, L_beam = 100.0e-06

Now, the code must run and give converged answer for all the Voltage_initial values mentioned above (as long as Voltage_initial is initialised with a small value like 0.1, 0.01 etc). Why this code is sensitive and gives response only for selective values of Voltage_initial must be known as there is certainly something happening that I am not able to understand.

This problem is basically on finding out static pull-in instability parameters of micro-cantilevers using three-mode reduced order model (ROM).

Answer 1

Solving nonlinear equation(s) is always prone to problems with the initial guess. While you can use bisection in 1-D, there is not an equivalent procedure in multi-dimension and you are at the mercy of the function's gradient at the initial guess. If it points you in the direction of the solution, all is well. But if it points to a place where the solution is lost (or to another root away from your desired root) you have to try something else.

Although using a linearization of the problem is a nice way to get an initial guess, after you might consider using the obtained solution of the nonlinear equations be your next guess for the changed value of voltage.

How do you know what the solution is supposed to be? Do you have a physical idea of the range of the solution? If so, you might consider probing randomly around that region. Here, for example, for V=5 I sample and collect the solutions using random() to generate initial guesses. (And I prefer to normalize the equations using factor which will make the coefficients range between -1 and 1.)

>>> nl = (Nonlinear_eq_1, Nonlinear_eq_2, Nonlinear_eq_3)
>>> NL = Tuple(*[factor(i).as_coeff_Mul()[1] for i in nl]) # make it well posed
>>> from random import random
>>> g = lambda : dict(zip(A,[random() for i in range(3)]))
>>> A = (a1,a2,a3)
>>> got = set()
>>> for i in range(100):
...    d = g()
...    try: sol = nsolve(NL.subs(V,5),A,[d[i] for i in A])
...    except:continue
...    got.add(tuple([i.n(4) for i in sol]))
>>> for i in got: print(i)
(0.6199, -0.05296, 0.003206)
(0.5670, 0.8258, 0.1807)
(1.552, -0.04830, 0.0005004)

Since all your equations are in the form ci*V**2 -bi, I wonder if the critical voltage is reached when there is no longer any values for all a_i than can keep all b_i/c_i positive. Maybe there would be a way to state the conditions that that ratio must be true that could help you identify the critical voltage.

The solution space (a1,a2,a3) contains "noodles", points on which are the solutions at different values of V. Say you find 4 points on 4 noodles at V=0.01. As you follow a solution while changing V, you may come to the end of a noodle and the solver is too far from other solutions to find one of them. With two functions scan(v) and follow(v1, v2, dv, g) you can scan with random values of a for a given V=v; with follow you could trace from V=v1 where you know the solution is g to V=v2 in steps of dv, using the new solution at each step as a new g:

def go(v, g):
    try: return tuple([i.n(4) for i in nsolve(
        NL.subs(V, v), A, g)])
    except: pass

def scan(v, n, *r):
    rv = set()
    for r in r:
        rvi = {go(v, [r*(2*random()-1) for _ in range(3)]) for i in range(n) }
        rv.update(rvi)
    return sorted(rv - {None})

def follow(v1, v2, dv, g):
    noodle = []
    for v in frange(v1, v2+dv, dv): # https://stackoverflow.com/a/67053708/1089161
        gg = go(v, g)
        if gg is None:
            break
        noodle.append((v, gg))
        g = gg
        v += dv
    return noodle

e.g.

>>> scan(10,20,.01,.1,1)
[(0.7634, -0.3434, 0.07579), (0.8307, -0.5453, 0.2540), (0.9054, -0.1020, 0.07110)]
>>> for v, a in follow(10,1,-1,_[0]):print(v,a)
...
10 (0.7634, -0.3434, 0.07579)
9 (0.7595, -0.3278, 0.05184)
8 (0.7559, -0.3139, 0.03266)
7 (0.7525, -0.3015, 0.01774)
6 (0.7490, -0.2896, 0.006241)
5 (0.7445, -0.2761, -0.002454)
4 (0.7360, -0.2544, -0.008533)
3 (0.6541, -0.09098, -0.004382)
2 (0.8146, -0.3254, -0.02020)
1 (0.8293, -0.2820, -0.02004)

Note that the solution jumped in going between V=3 and V=2.