我正试图解决偏微分方程enter image description here
$\frac{\partial \textbf{u}}{\partial t}=b\sin{\alpha}+ \nu \frac{\partial^2 \textbf{u}}{\partial z^2}$
$\frac{\partial \textbf{b}}{\partial t}=-\textbf{u} N^2\sin{\alpha}+ \frac{\nu}{Pr} \frac{\partial^2 \textbf{b}}{\partial z^2}$
我将PED转换为一个ODE系统,但它无论如何都不会收敛。我使用了每一种解决ivp的方法,但仍然没有收敛,我的代码有什么问题吗?因为它应该有解决方案
我尝试了许多方法,但大多数方法都不收敛,并且会显示溢出
`
def func(t,y):
n = len(y)
# set u,b,k,w as variables
u=y
b=y
# 2 equations initialization
f1=np.zeros(nz)
f2=np.zeros(nz)
#1st equation with BCs
f1[1]= b[1]*m.sin(alpha)+(Nu*(u[2]-2*u[1]+0)/dz**2)
f2[1]= -u[1]*N**2*m.sin(alpha)+(Nu*(b[2]-2*b[1]+bs)/dz**2)
for i in range(2,nz-2):
f1[i]=b[i]*m.sin(alpha)+((Nu)*(u[i+1]-2*u[i]+u[i-1])/dz**2)
f2[i]=-u[i]*N**2*m.sin(alpha)+((Nu)*(b[i+1]-2*b[i]+b[i-1])/dz**2)
f1[nz-2]= b[nz-2]*m.sin(alpha)+(Nu*(0-2*u[nz-2]+u[nz-3])/dz**2)
f2[nz-2]= -u[nz-2]*N**2*m.sin(alpha)+(Nu*(0-2*b[nz-2]+b[nz-3])/dz**2)
f1=f1[1:]
f1=f1[:-1]
f2=f2[1:]
f2=f2[:-1]
return np.append(f1,f2)
global Beta, alpha, Nu, dz, N, BetaS, Beta, AlphaS, SigmaS, bs #needed in equations
H = 800
T = 1200
alpha = m.radians(90) # sloping angle [rad]
bs = 8.26e-04 # surface buoyancy [m/s^2]
N = 0.006 # buoyancy frequency [Hz]
Nu = 1e-05 # kinematic viscosity [m^2/s]
Pr = 1 # Prandtl number
nz = 1200 # number of elements
nt = 6000 # number of time steps
dz = H/nz # grid stencil [m]
dt = T/nt # time step [s]
z=np.linspace(0,H,nz)
print(dt)
# I.Cs
S0 = m.sqrt(N*N*m.sin(alpha)**2/Pr)
SC = m.sqrt(S0/2/Kb)#*100.25))
u_0=np.zeros(nz)
b_0=np.zeros(nz)
"""Initial condition for U,b Prantdl model analytical solution with constant viscosity"""
for n in range(nz):
u_0[n] = (bs/N/m.sqrt(Pr))*m.exp(-SC*z[n])*m.sin(SC*z[n]);
b_0[n] = bs*m.exp(-SC*z[n])*m.cos(SC*z[n]);
u_0=u_0[1:]
u_0=u_0[:-1]
b_0=b_0[1:]
b_0=b_0[:-1]
s=z[1:]
s=s[:-1]
plt.plot(u_0, s, label='u')
plt.legend()
plt.show()
plt.plot(b_0, s, label='b')
plt.legend()
plt.show()
IC=np.append(u_0,b_0)
#Record start time
start_time=time.time()
YS=solve_ivp(func, [0,T], IC,method='LSODA', t_eval=np.arange(0,T,dt), rtol=1e-3)
#Record end time
end_time=time.time()
默认方法将立即显示溢出,我唯一可以使用的方法是'LSODA',但它将运行很长时间而不会结束
目前没有回答
相关问题 更多 >
编程相关推荐