我尝试使用以下代码为我的一个项目制作几个图：

import pprint
import scipy
import scipy.linalg   # SciPy Linear Algebra Library
import numpy as np
from scipy.linalg import lu , lu_factor, lu_solve
from scipy.integrate import quad
import matplotlib.pyplot as plt

# Solving the equations for the Prandtl case

K = 100
alpha = 0.1 
visc = 5  
diff = 5   
N = 0.01     
L = 5000
height = 250 

subdivisions = 100
tick = 10
points = np.arange(0,L/2+tick,tick)

def H(y):
    return ( height * (1 + np.cos(2 * np.pi * y/L)) )

def Bsfc(y):
    return 0.1

final_system = []
b=[]
for q in range(-K,K+1):
    equation1 = []
    equation2 = []
    equation3 = []
    Aki = []
    Cki = []
    Dki = []
    for k in range(-K,K+1):

        R = 2 * N**2 * np.cos(alpha)**2 / (visc * diff) * (k * np.pi / L)**2

        Q = N**2 * np.sin(alpha)**2 / (3 * visc * diff)

        S1 = abs(R + np.sqrt(Q**3 + R**2) )**(1/3)
        S2 = - abs( np.sqrt(Q**3 + R**2) -R )**(1/3)

        phi = np.sqrt(S1**2 + S2**2 - S1*S2)
        Lk = np.arccos(- (S1 + S2)/ (2 * phi) )

        m1 = - np.sqrt(S1 + S2)
        m2 = - np.sqrt(phi) * np.exp(1j * Lk/2)
        m3 = m2.conjugate()

        def f1r(y):
            return (np.exp(m1 * H(y)) * np.cos(2 * (q - k) * np.pi * y / L) ).real
        def f1i(y):
            return (np.exp(m1 * H(y)) * np.cos(2 * (q - k) * np.pi * y / L) ).imag
        gamma1 = 2/L * (quad(f1r,0,L/2,limit=subdivisions)[0] + quad(f1i,0,L/2,limit=subdivisions)[0]*1j)

        def f2r(y):
            return (np.exp(m2 * H(y)) * np.cos(2 * (q - k) * np.pi * y / L) ).real
        def f2i(y):
            return (np.exp(m2 * H(y)) * np.cos(2 * (q - k) * np.pi * y / L) ).imag
        gamma2 = 2/L * (quad(f2r,0,L/2,limit=subdivisions)[0] + quad(f2i,0,L/2,limit=subdivisions)[0]*1j)

        if k == 0:
            equation1.append(2 * gamma2.real)
            Cki.append(k)
            equation1.append(-2 * gamma2.imag)
            Dki.append(k)
        else:
            equation1.append(gamma1)
            Aki.append(k)
            equation1.append(2 * gamma2.real)
            Cki.append(k)
            equation1.append(-2 * gamma2.imag)
            Dki.append(k)

        if q != 0:

            if k == 0:
                equation2.append(0)
                equation2.append(0)
            else:
                equation2.append(k * gamma1 / (m1**3) )
                equation2.append(2 * k * (gamma2 / (m2**3) ).real)
                equation2.append(-2 * k * (gamma2 / (m2**3) ).imag)

        if k == 0:
            equation3.append(2 * (m2**2 * gamma2).real)
            equation3.append(-2 * (m2**2 * gamma2).imag)
        else:
            equation3.append(m1**2 * gamma1)
            equation3.append(2 * (m2**2 * gamma2).real)
            equation3.append(-2 * (m2**2 * gamma2).imag)

    final_system.append(equation1)
    def f4r(y):
        return (Bsfc(y) * np.cos(2 * q * np.pi * y / L) ).real
    def f4i(y):
        return (Bsfc(y) * np.cos(2 * q * np.pi * y / L) ).imag
    b.append(2/L * (quad(f4r,0,L/2,limit=subdivisions)[0] + quad(f4i,0,L/2,limit=subdivisions)[0]*1j))

    if q != 0:
        final_system.append(equation2)
        b.append(0)

    final_system.append(equation3)
    b.append(0)

final_system = np.array(final_system)
b=np.array(b)

# LU solver

P, Ls, U = scipy.linalg.lu(final_system)

Bl = np.linalg.inv(P) @ b 

Z = np.linalg.solve(Ls,Bl)

X = np.linalg.solve(U,Z)

print (np.allclose(final_system @ X, b))

# Getting the values for Ak, Ck and Dk

strings = []

for k in range(-K,K+1):
    if k != 0:
        strings.append('A')     

    strings.append('R')           
    strings.append('I')             

Ak = []
Rk = []
Ik = []
for k in range(0,len(X)):
     if 'A' in strings[k]:
         Ak.append(X[k])

     if 'R' in strings[k]:
         Rk.append(X[k])

     if 'I' in strings[k]:
         Ik.append(X[k])

Ck=[]

for k in range(0,len(Rk)):
    Ck.append(Rk[k] + Ik[k] * 1j)

Ck = np.array(Ck)

Dk = Ck.conjugate()

Ak = np.array(Ak)

# Getting the Buoyancy value

z = np.arange(0,2010,10) 
y = np.arange(-L,L+10,10) 
Y,Z = np.meshgrid(y,z)
B = np.ones_like(Y)*[0]

for k in range(-K,K+1):

    R = 2 * N**2 * np.cos(alpha)**2 / (visc * diff) * (k * np.pi / L)**2

    Q = N**2 * np.sin(alpha)**2 / (3 * visc * diff)

    S1 = abs(R + np.sqrt(Q**3 + R**2) )**(1/3)
    S2 = - abs( np.sqrt(Q**3 + R**2) -R )**(1/3)

    phi = np.sqrt(S1**2 + S2**2 - S1*S2)
    Lk = np.arccos(- (S1 + S2)/ (2 * phi) )

    m1 = - np.sqrt(S1 + S2)
    m2 = -np.sqrt(phi) * np.exp(1j * Lk/2)
    m3 = m2.conjugate()

    if k != 0:
        B = B + ( Ak[Aki.index(k)] * np.exp(m1 * Z) * np.exp(2j * (k) * np.pi * Y / L)  )
    B = B + ( ( Ck[Cki.index(k)] * np.exp(m2 * Z) + Dk[Dki.index(k)] * np.exp(m3 * Z) )  * np.exp(2j * (k) * np.pi * Y / L) )

for k in range(0,B.shape[0]):
    for t in range(0,B.shape[1]):
        if Z[k][t] < H(Y[k][t]):
            B[k][t] = np.nan
        if Z[k][t] == H(Y[k][t]):
            print (B[k][t], "B value at the ground")
        if abs(Z[k][t] - H(Y[k][t])) < 0.1:
            if B[k][t] > 0.101:
                print (B[k][t],'error -------------------------------------------------')

# print (B[k][t], Z[k][t], H(Y[k][t]), Y[k][t], '-----------------------------------------------------------------------------' )

Bp = Bsfc(Y) * np.exp(-Z * np.sqrt(N * np.sin(alpha) ) / (4*visc*diff)**(1/4) ) * np.cos(np.sqrt(N*np.sin(alpha)) /((4*visc*diff)**(1/4))*Z )

## Plotting the buoyancy

fig = plt.figure(figsize=(10,10)) # create a figure
plt.rcParams.update({'font.size':16})
plt.title('Buoyancy')
plt.contourf(Y,Z,B,np.arange(-0.2,0.201,0.001),cmap='seismic')

# plt.contourf(Y,Z,B,cmap='seismic')

plt.colorbar(label='1/s')
plt.xlabel("Y axis")
plt.ylabel("Height")
plt.xlim([-L,L])
plt.ylim([0,1500])
plt.show()

下图显示了产生良好结果的运行：
Buoyancy
然而，当我增加“height”参数时，我开始得到不稳定的结果，我怀疑这是由于数值不稳定性造成的：
Buoyancy unstable
有没有办法提高python中的数值精度呢？我已经用numpy.double做了一些实验，但是到目前为止没有成功的结果。
谢谢