I'm trying to simulate in time and space the following system of partial differential equations. I'm using python 3 for that.
Here is a link to the set of equations with their boundary conditions
My ideas was to transform all the equations to the discrete form (forward Euler as the simplest starting point) and then run the code. Forward Euler implies: Here lin to image i = 0,...,Nx - mesh for n = 0,1,...,Nt Here what I have (by the means of 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
My first problem I'm encountering is different warnings:
/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
My main question is: it possible to solve this set with the forward Euler Method Iam trying at the moment ? Thank you everyone in advance!
The answer is "yes", but your code needs more work. For example, you will need to work on the algorithm's stability (to avoid it blowing up). Also, the BC does not reflect your system; I think you are looking for zero flux conditions; if so, you are not coding it right. Finally, you can also consider using Fipy, which could make your life easier. Please take a look here https://www.ctcms.nist.gov/fipy/ I also wrote a basic example here http://biological-complexity.blogspot.pe/