i try to solver $u_{t}=u_{xx}+u_{yy}+u-u^2 in(-a,a)times(-b,b)times(0,T)with initial condition u_0=1 if(x,y)=0and 0 other wise and u=0 in$partialomega$ where $omega=[-a;a]times[-b;b], i do that by finite difference method
this is the code i use :
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps):
u = np.zeros((Nx, Ny))
u[Nx // 2, Ny // 2] = 1 # Set initial condition
Dxx, Dyy = 1/dx**2, 1/dy**2
i = slice(1, Nx-1)
iw = slice(0, Nx-2)
ie = slice(2, Nx)
j = slice(1, Ny-1)
js = slice(0, Ny-2)
jn = slice(2, Ny)
for _ in range(steps):
dux = Dxx * (u[ie, j] - 2*u[i, j] + u[iw, j])
duy = Dyy * (u[i, jn] - 2*u[i, j] + u[i, js])
reaction = r * u[i, j] * (1 - u[i, j])
u[i, j] += dt * (D * (dux + duy) + reaction)
return u
# Parameters
Nx = 100 # Number of mesh points in the x-direction
Ny = 100 # Number of mesh points in the y-direction
D = 0.1 # Diffusion coefficient
r = 1.0 # Reaction rate
dx = 1 / Nx
dy = 1 / Ny
dt = (dx**2 * dy**2) / (2 * D * (dx**2 + dy**2))
# Time steps for t = 5, t = 10, t = 20
steps_t5 = int(5 / dt)
steps_t10 = int(10 / dt)
steps_t20 = int(20 / dt)
# Solve the Fisher-KPP equation for t = 5, t = 10, and t = 20
solution_t5 = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_t5)
solution_t10 = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_t10)
solution_t20 = solve_fisher_kpp(Nx, Ny, D, r, dx, dy, dt, steps_t20)
# Create meshgrid for plotting
x = np.linspace(0, 1, Nx)
y = np.linspace(0, 1, Ny)
X, Y = np.meshgrid(x, y)
# Plotting the solution in 3D
fig = plt.figure(figsize=(18, 6))
# Plot at t = 5
ax1 = fig.add_subplot(131, projection='3d')
Z_t5 = solution_t5.T # Transpose for correct orientation in plot
ax1.plot_surface(X, Y, Z_t5, cmap='viridis')
ax1.set_title('Fisher-KPP Solution at t = 5')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('U')
# Plot at t = 10
ax2 = fig.add_subplot(132, projection='3d')
Z_t10 = solution_t10.T # Transpose for correct orientation in plot
ax2.plot_surface(X, Y, Z_t10, cmap='viridis')
ax2.set_title('Fisher-KPP Solution at t = 10')
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('U')
# Plot at t = 20
ax3 = fig.add_subplot(133, projection='3d')
Z_t20 = solution_t20.T # Transpose for correct orientation in plot
ax3.plot_surface(X, Y, Z_t20, cmap='viridis')
ax3.set_title('Fisher-KPP Solution at t = 20')
ax3.set_xlabel('X')
ax3.set_ylabel('Y')
ax3.set_zlabel('U')
plt.show()
my problem is that it doesnt run ,where is the issue and how ca,n i fix it? Thanks