我编写了以下代码，用于求三个非线性代数方程组**（Nonlinear_eq_1、Nonlinear_eq_2和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代码给出了以下公式，为方便起见，直接在下面生成）

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])

当电压增加时，尝试遵循以下解决方案：

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

目的是使用最终非线性答案计算捕捉值：

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)

系统在临界电压值以下是稳定的，而在临界电压值以上变得不稳定。

按照代码进行翻译如下：
*当V〈临界电压时，非线性代数方程组使用nsolve求解。
*当V超过临界电压时，非线性代数方程组使用nsolve给出'ValueError'。
我将此逻辑与 try-except 沿着使用，以迭代方式找到临界电压值。
现在我观察到的关于此代码的奇怪行为如下：
*对于电压初始值= 0.1或0.001或0.000001，代码未给予g_0 = 5.0e-06、h梁= 1.80e-06、b梁棱柱= 15.0e-06、L梁= 100.0e-06的收敛答案
*然而，对于电压初始值= 0.01或0.0001或0.00001，相同的代码给出了g_0 = 5.0e-06、h梁= 1.80e-06、b梁棱柱= 15.0e-06、L梁= 100.0e-06的收敛答案
现在，代码必须运行，并对上述所有Voltage_initial值给予收敛答案（只要Voltage_initial初始化为0.1、0.01等小值）。必须了解为什么此代码敏感，并仅对Voltage_initial的选择性值给出响应，因为肯定发生了我无法理解的事情。

该问题的实质是利用三模态降阶模型（ROM）求出微悬臂梁的静态pull-in失稳参数。