二次元拡散モデルのシミュレーション実行時にRuntime　warningが出る

前提・実現したいこと

Google Colabで二次元拡散方程式を可視化しようとしています。
シミュレーション実行中に以下のエラーメッセージが発生しました。

発生している問題・エラーメッセージ

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in double_scalars
  import sys

該当のソースコード

Python
1#01-02-01 二次元の拡散方程式の関数定義
2
3import numpy as np
4import matplotlib.pyplot as plt
5
6def diff_eq(u_arr,v_arr,Du,Dv,du,dv,dh,dt,k1,k2,k3):
7  new_u=u_arr[1,1]+(Du*(u_arr[0,1]+u_arr[2,1]+u_arr[1,0]+u_arr[1,2]-4*u_arr[1,1])/(dh**2)-du*u_arr[1,1]+k1*(u_arr[1,1]*u_arr[1,1]/v_arr[1,1])+k2)*dt
8  new_v=v_arr[1,1]+(Dv*(v_arr[0,1]+v_arr[2,1]+v_arr[1,0]+v_arr[1,2]-4*v_arr[1,1])/(dh**2)-dv*v_arr[1,1]+k3*u_arr[1,1]*u_arr[1,1])*dt
9  return new_u,new_v
10
11
12#01-02-02 二次元の拡散方程式に基づく更新
13def update(fieldu,fieldv,Du,Dv,du,dv,dh,dt,k1,k2,k3):
14  nu_i,nu_j=fieldu.shape
15  new_fieldu= np.zeros((nu_i,nu_j),dtype=fieldu.dtype)
16  nv_i,nv_j=fieldv.shape
17  new_fieldv= np.zeros((nv_i,nv_j),dtype=fieldv.dtype)
18  for i in range(nu_i):
19    for j in range(nu_j):
20      if i==0:
21        if j==0:
22          #境界条件処理 i=0 , j=0
23          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[-1,0,1],:][:,[-1,0,1]],fieldv[[-1,0,1],:][:,[-1,0,1]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
24        elif 0<j<nu_j-1:
25          #境界条件処理i=0,0<j<n_j-1
26          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[-1,0,1],:][:,[j-1,j,j+1]],fieldv[[-1,0,1],:][:,[j-1,j,j+1]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
27        elif j==nu_j-1:
28          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[-1,0,1],:][:,[nu_j-2,nu_j-1,0]],fieldv[[-1,0,1],:][:,[nu_j-2,nu_j-1,0]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
29      elif 0< i <nu_j-1:
30        if j==0:
31          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[i-1,i,i+1],:][:,[-1,0,1]],fieldv[[i-1,i,i+1],:][:,[-1,0,1]],Du,Dv,du,dv,dh,dt,k1,k2,k3)  
32        elif 0<j<nu_j-1:
33          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[i-1,i,i+1],:][:,[j-1,j,j+1]],fieldv[[i-1,i,i+1],:][:,[j-1,j,j+1]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
34        elif j==nu_j-1:
35          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[i-1,i,i+1],:][:,[nu_j-2,nu_j-1,0]],fieldv[[i-1,i,i+1],:][:,[nu_j-2,nu_j-1,0]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
36      elif i == nu_i-1:
37        if j ==0:
38          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[nu_i-2,nu_i-1,0],:][:,[-1,0,1]],fieldv[[nu_i-2,nu_i-1,0],:][:,[-1,0,1]],Du,Dv,du,dv,dh,dt,k1,k2,k3)          
39        elif 0<j<nu_j-1:
40          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[nu_i-2,nu_i-1,0],:][:,[j-1,j,j+1]],fieldv[[nu_i-2,nu_i-1,0],:][:,[j-1,j,j+1]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
41        elif j==nu_j-1:
42          new_fieldu[i,j],new_fieldv[i,j]=diff_eq(fieldu[[nu_i-2,nu_i-1,0],:][:,[nu_j-2,nu_j-1,0]],fieldv[[nu_i-2,nu_i-1,0],:][:,[nu_j-2,nu_j-1,0]],Du,Dv,du,dv,dh,dt,k1,k2,k3)
43
44  return new_fieldu, new_fieldv
45#01-02-03 二次元拡散モデルのシミュレーション実行
46
47Du=14.0
48Dv=3.0
49du=0.6
50dv=0.2
51dt=0.01
52dh=1
53k1=0.3
54k2=0.2
55k3=0.9
56
57u_0=25
58v_0=1800
59
60x=np.arange(-25,25,dh)
61y=np.arange(-25,25,dh)
62xmesh,ymesh=np.meshgrid(x,y)
63
64t_end=50
65
66t_list=[]
67u_list=[]
68v_list=[]
69
70
71u=np.zeros([len(x),len(y)],dtype=float)
72v=np.zeros([len(x),len(y)],dtype=float)
73u[25,25]=u_0
74v[25,25]=v_0
75
76
77t_list.append(0)
78u_list.append(np.copy(u))
79v_list.append(np.copy(v))
80
81for n in range(int(t_end/dt)):
82  t= (n+1)*dt #時刻ｔの計算
83
84  u_tmp,v_tmp=update(u,v,Du,Dv,dh,dt,du,dv,k1,k2,k3)
85
86  t_list.append(t)
87  
88  u=np.copy(u_tmp)
89  v=np.copy(v_tmp)
90  u_list.append(np.copy(u))
91  v_list.append(np.copy(v))
92