我不理解如何使用python或python ode求解器求解耦合的PDE方程中的eta和V。 (或者是否可以在没有求解器的情况下对这些耦合方程进行数值求解?)我已经花了几天时间,但是我仍然不知道如何开始!任何提示都会有所帮助。我了解[
中给出的示例Solve 4 coupled differential equations in MATLAB
但是我仍然需要更多提示来说明如何将这些原理应用于下面的耦合PDE。
我想绘制一个eta和V的时间序列,给定强迫输入tau的时间序列。
x是空间点,t是时间点。 h和f是根据x的值分配的。V = V(x,t)eta = eta(x,t)tau = tau(x,t)h = h(x),f = f(x)而g和rho是常数。边界值为V(0,0)= 0,eta(0,0)= 0和tau(0,0)=0。假设处于稳态条件,则通过将tau_sy和tau_by相等来找到解(V)。
好,所以这是一个稍微简单的数字方案,它显示了系统的概念属性。它类似于(显式)欧拉方法。可以很容易地将其推广为类似的隐式欧拉式方法。
您被赋予:
功能h(x)
,f(x)
,tau_sx(x, t)
,tau_sy(x, t)
和tau_by(x, t)
常数g
和rho
您正在寻找:
满足以上一对微分方程的函数V(x, t)
和eta(x, t)
。
为了能够找到此问题的解决方案,您需要提供:
[V(x, 0) = V0(x)
和eta(0, t) = eta0(t)
假设您的域是[0, L] X [0, T]
,其中x
中的[0, L]
和t
中的[0, T]
。按如下所示离散域:选择M
和N
正整数,并让dx = L / M
和dt = T / N
。然后,对于任何整数x = m dx
和t = n dt
,仅考虑点m = 0, 1, ..., M
和n = 0, 1, ..., N
的有限组。
我将把所有函数限制在上面定义的有限点上,并对任意函数funct
使用以下符号:
对于funct(x, t) = funct[m, n]
和funct(x) = funct[m]
中的任何一个,[x = m dx
和t = n dt
。
然后,微分方程组可以离散为
g*(h[m] + eta[m,n])*(eta[m+1, n] - eta[m,n])/dx = f[m]*(h[m] + eta[m,n])*V[m,n] + tau_sx[m,n]/rho
(V[m, n+1] - V[m,n])/dt = (tau_sy[m,n] - tau_by[m,n])/(rho*(h[m] + eta[m,n]))
解决eta[m+1,n]
和V[m,n+1]
eta[m+1,n] = eta[m,n] + f[m]*V[m,n]*dx/g + tau_sx[m,n]*dx/(g*rho*(h[m] + eta[m,n]))
V[m,n+1] = V[m,n] + (tau_sy[m,n] - tau_by[m,n])*dt/(rho*(h[m] + eta[m,n]))
为简单起见,我将上述等式的右边缩写为
eta[m+1,n] = F_eta(m, n, eta[m,n], V[m,n])
V[m,n+1] = F_V(m, n, eta[m,n], V[m,n])
即类似
def F_eta(m, n, eta[m,n], V[m,n]):
return eta[m,n] + f[m]*V[m,n]*dx/g + tau_sx[m,n]*dx/(g*rho*(h[m] + eta[m,n]))
def F_V(m, n, eta[m,n], V[m,n]):
return V[m,n] + (tau_sy[m,n] - tau_by[m,n])*dt/(rho*(h[m] + eta[m,n]))
根据边界条件,我们知道
eta[0,n] = eta0[n] = eta0(n*dt)
和
V[m,0] = V0[m] = V0(m*dx)
作为输入,对于m = 0,..., M
和n = 0,..., N
。
for n in range(N):
for m in range(M):
eta[m+1,n] = F_eta(m, n, eta[m,n], V[m,n])
V[m,n+1] = F_V(m, n, eta[m,n], V[m,n])
((您必须调整这些循环才能到达最右边和最上面的边界点,但是原理保持不变)
[基本上,您遵循此模式:沿水平x轴生成eta,同时向上生成一层V。然后您向上移动到下一个水平位置。
o --eta--> o --eta--> o --eta--> o --eta--> o
| | | | |
V V V V V
| | | | |
o --eta--> o --eta--> o --eta--> o --eta--> o