如何在 python 中模拟耦合 PDE

Question

我正在尝试及时模拟 space 以下偏微分方程组。我为此使用 python 3。

Here is a link to the set of equations with their boundary conditions

我的想法是将所有方程转化为离散形式（正向欧拉作为最简单的起点），然后运行代码。正向欧拉意味着： Here lin to image i = 0,...,Nx - mesh for n = 0,1,...,Nt 这是我所拥有的（通过 numpy）

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
#Define exponents for PDE
m = 0
n = 2
#Define constants for PDE
a = 0.2
b= -0.4
av = 5.0
c = 0.6
d = -0.8
Du = 1
Dv = 20
Dz = 1000
u0 = 0.5
v0 = 0.5
kz = 0.001

L = 10
Nx = 100
T = 5
Nt = 100
x = np.linspace(0, L, Nx+1)
dx = x[1] - x[0]
#print(dx)
#print(dt)
t = np.linspace(0, T, Nt+1)
dt = t[1] - t[0]
if dt<=0.5*dx**2:
    print("Ok!")
else:
    print("Alert! dt is not smaller than dx^2/2")
u = np.zeros(Nx+1)
v = np.zeros(Nx+1)
z = np.zeros(Nx+1)
u_1 = np.zeros(Nx+1)
v_1 = np.zeros(Nx+1)
z_1 = np.zeros(Nx+1)
# mesh points in space
# mesh points in time
# Set initial condition u(x,0) = I(x)
for i in range(0, Nx+1):
    u_1[i] = np.random.random_sample()
    v_1[i] = np.random.random_sample()
    z_1[i] = np.random.random_sample()
for n in range(0, Nt):
    # Compute u at inner mesh points
    for i in range(1, Nx):
        u[i] = u_1[i] + dt*(a*(u_1[i]-u0) + 
b*(v_1[i]-v0)+av*(u_1[i]-u0)**3+(Du/dx**2)*(u_1[i-1] - 
2*u_1[i] + u_1[i+1]))*z_1[i]**n
    v[i] = v_1[i] + dt*(c*(u_1[i]-u0)+d*(v_1[i]-v0)+(Dv/dx**2)*(v_1[i-1] - 2*v_1[i] + v_1[i+1]))*z_1[i]**n
    z[i] = (Dz/dx**2)*((z_1[i-1] - 2*z_1[i] + z_1[i+1]) - kz * z[i])
    # Insert boundary conditions u[0]=0; u[Nx]=0
    u[0]=0; u[Nx]=1/Dz
    v[0]=0; v[Nx]=1
    z[0]=0; z[Nx]=1
    # Update u_1 before next step
    u_1[:]= u
    v_1[:]= v
    z_1[:]= z

我遇到的第一个问题是不同的警告：

/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:31: RuntimeWarning: overflow encountered in double_scalars
   /miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:31: RuntimeWarning: invalid value encountered in double_scalars
    /miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: invalid value encountered in double_scalars
    /miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:32: RuntimeWarning: overflow encountered in double_scalars
    /miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:33: RuntimeWarning: overflow encountered in double_scalars
    /miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:33: RuntimeWarning: invalid value encountered in double_scalars

我的主要问题是：我现在正在尝试的正向欧拉法可以解决这个问题吗？先谢谢大家了！

Answer 1

答案是“是”，但您的代码还需要改进。例如，您需要研究算法的稳定性（以避免它崩溃）。此外，BC 不反映您的系统；我认为您正在寻找零通量条件；如果是这样，则说明您的编码不正确。最后，您还可以考虑使用 Fipy，它可以让您的生活更轻松。请看这里https://www.ctcms.nist.gov/fipy/ I also wrote a basic example here http://biological-complexity.blogspot.pe/

如何在 python 中模拟耦合 PDE

How to simulate coupled PDE in python

python

numpy

pde