2017-05-02 7 views
0

私は以下の偏微分方程式のシステムを時空間でシミュレートしようとしています。私はそのためにPython 3を使用しています。Pythonで連結されたPDEをシミュレートする方法

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

私のアイデアは、離散形式(最も簡単な出発点として前進オイラー)にすべての方程式を変換してからコードを実行することでした。 前進オイラーは意味:

/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 

私の主な質問は次のとおりです:ここで私が持っているもの 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 

私が遭遇してる私の最初の問題は、異なる警告であることが前方オイラー法Iamでこのセットを解くことは現時点では可能ですか? 皆さんありがとうございます!

答えて

0

答えは「はい」ですが、コードはより多くの作業が必要です。たとえば、アルゴリズムの安定性に取り組む必要があります(爆発を避けるため)。 BCもあなたのシステムを反映していないようですが、もしそうなら、あなたはゼロフラックス条件を探していると思います。最後に、Fipyを使用することも考えられます。これはあなたの人生をより簡単にすることができます。ここをクリックしてくださいhttps://www.ctcms.nist.gov/fipy/と私はここに基本的な例を書いていますhttp://biological-complexity.blogspot.pe/

関連する問題