Python 中的 PolyeigOrr-Sommerfeld 方程

如何解决Python 中的 PolyeigOrr-Sommerfeld 方程

我正在做一个小项目，试图使用 Python 确定 Orr-Sommerfeld 方程的特征值，这是特定的方程：

在特定的边界条件下：

注意到的数量 U 是通过解决 Blasius 方程（我求解）推导出的速度。我们还考虑以下数量：

? = 0、? = 0.125、?? = 520 和 ? = ∂/∂?。

问题的目的是找到?，为此建议将问题写成如下：

? = ?1 + ?2 ? + ?3 ?^2 + ?4 ?^3 + ?5 ?^4

知道这是一个特征值问题，我试着写了这段代码来解决它：

ny=99
long=70
y=np.linspace(-long,long,ny+2)
beta=0
I=np.identity(ny+1)
w=0.125
Re=520
dy=long/ny

# We discretize our operators

D2 = sp.diags([1,-2,1],[-1,(ny+1,ny+1) )/dy**2 
D4=sp.diags([1,-4,6,[-2,-1,1,2],ny+1) )/dy**4

# Here we use the Blasius resolution with the boundary value problem solver of Python

U=res.y[0] 
U2=res.y[2]

然后我们像这样构造矩阵：

A=U2 #alpha
B=-U #alpha**3
C=-w*D2 #alpha**0
D=D2*U #alpha
E=w*I #alpha**2
F=(1/(Re*1j))*D4 #alpha**0
G=(1/(Re*1j))*I #alpha**4
H=-(1/(Re*1j))*2*D2 #alpha**2
M1=C+F
M2=A+D
M3=H+E
M4=B
M5=G
V=A+B+C+D+E+F+G+H

通常情况下，使用 Matlab 我知道我们可以应用 Matlab 的“polyeig”求解器来实现这个问题，但我认为 Python 中没有直接的等价物，所以我试着写这样的问题：通过尝试带回问题的一维表示（两个矩阵 L 和 K 的目标）：

L1=np.array([[-M1,0],[0,0]])
Z=np.zeros((2,2))
L=np.block([[L1,Z],[Z,Z]])
B1=np.array([[M2,M3]])
B2=np.array([[M4,M5]])
K=np.block([[B1,B2]])

#V=(C+F)+(alphar +1j*alphai)*(A+D)+((alphar +1j*alphai)**2)*(H+E) +((alphar +1j*alphai)**3)*B +((alphar +1j*alphai)**4)*G

# We implement the values of the boundary conditions

V[0]=0
V[ny]=0
V[1]=V[0]
V[ny]=V[ny-1]

alpha=linalg.eigvals(V)
plt.plot(alpha.real,alpha.imag,'.')
plt.xlim(-0.1,0.4)
plt.ylim(-0.2,0.5)

不幸的是它不起作用，因为当我将“linalg.eigvals”应用于 L 和 K 时，我遇到了错误...

ValueError: object arrays are not supported 

有人可以建议我如何在我的情况下使用 Matlab 的“polyeig”替代品吗？

非常感谢您的帮助。

PS：这是用于 Blasius 方程的求解器，以便完全运行代码：

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

    
def blasius(t,u):
    F,dF,ddF = u
    return [dF,ddF,-0.5*(F*ddF)]

def cond(u0,u1): return [u0[0],u0[1],1.0-u1[1]]

x = np.linspace(0,10,100)
u = [x,np.exp(-x),np.exp(-x)]

res = solve_bvp(blasius,cond,x,u)
    
plt.plot(res.x,res.y[0],label="$f(\eta)$")
plt.xlabel("$\eta$")
plt.ylabel("f")
plt.title("Blasius solution")
plt.grid()