トミー PDFからテキスト出力してcsv化
Python
# ファイル名 print_pdf_textboxes_csv.py
import sys
import csv
import os

from pdfminer.converter import PDFPageAggregator
from pdfminer.layout import LAParams, LTContainer, LTTextBox
from pdfminer.pdfinterp import PDFPageInterpreter, PDFResourceManager
from pdfminer.pdfpage import PDFPage

’’’
python print_pdf_textboxes_csv.py [pdfのパス]のようにコマンドプロンプトから実行する
’’’

# 再帰的にテキストボックス(LTTextBox)を探して、テキストボックスのリストを取得する
def find_textboxes_recursively(layout_obj):

    # LTTextBoxを継承するオブジェクトの場合は1要素のリストを返す
    if isinstance(layout_obj, LTTextBox):
        return [layout_obj]

    # LTContainerを継承するオブジェクトは子要素を含むので、再帰的に探す
    if isinstance(layout_obj, LTContainer):
        boxes = []
        for child in layout_obj:
            boxes.extend(find_textboxes_recursively(child))

        return boxes

    # その他の場合は空リストを返す
    return []


# Layout Analysisのパラメーターを設定縦書きの検出を有効にする
laparams = LAParams(detect_vertical=True)

# 共有のリソースを管理するリソースマネージャーを作成
resource_manager = PDFResourceManager()

# ページを集めるPageAggregatorオブジェクトを作成
device = PDFPageAggregator(resource_manager, laparams=laparams)

# Interpreterオブジェクトを作成
interpreter = PDFPageInterpreter(resource_manager, device)

csv_data = []

if len(sys.argv) == 1:
  print('パスを指定してください')
else:
  # 本来、パラメータがファイルパスであるかとか存在確認とかしないといけないが割愛
  pdfpath = sys.argv[1]
  base = os.path.splitext(pdfpath)[0]
  newpath = base + ".csv"

  # ファイルをバイナリ形式で開く
  with open(pdfpath, "rb") as f:

      # PDFPage.get_pages()にファイルオブジェクトを指定して、PDFPageオブジェクトを順に取得する
      # 時間がかかるファイルは、キーワード引数pagenosで処理するページ番号(0始まり)のリストを指定するとよい
      for i, page in enumerate(PDFPage.get_pages(f), start=1):

          # ページを処理する
          interpreter.process_page(page)

          # LTPageオブジェクトを取得
          layout = device.get_result()

          # ページ内のテキストボックスのリストを取得する
          boxes = find_textboxes_recursively(layout)

          # テキストボックスの左上の座標の順でテキストボックスをソートする
          # y1(Y座標の値)は上に行くほど大きくなるので、正負を反転させている
          boxes.sort(key=lambda b: (-b.y1, b.x0))

          for box in boxes:
              # ページ数、X座標、Y座標、テキスト
              csv_data.append([i, box.x0, box.y1, box.get_text().strip()])

  # 文字コードがUTF-8 With BOMだと、エクセルで開いたときに文字化けしません
  # テキストで開くと文字化けする・・・
  # BOM(byte order mark) Unicodeの符号化形式で符号化したテキストの先頭につける数バイトのデータのこと
  # encoding="utf_8_sig"とすることでBOM付きのUTF-8でCSVファイルを書き出すことが出来る
  with open(newpath, "wt", encoding="utf_8_sig") as fw:
      writer = csv.writer(fw, lineterminator="\n")
      writer.writerows(csv_data)
anonymous No title
Python
def to_bin(x,ans):
    #print(x, ans)
    if x == 0:
        return ''.join(reversed(list(map(str,ans))))
    else:
        ans.append(x%2)
        return to_bin(x//2, ans)

while True:
    v = int(input())
    print(to_bin(v, []))
anonymous No title
Python
import json
import re
from pprint import pprint

# MEMO
# CapsLock → F21
# F21 → 132


def configure(keymap):
    # *---jsonファイルの読み込み
    with open('keymaps.jsonc', encoding="utf-8") as f:
        raw_data = f.read()
    raw_data_comments_ignored = re.sub(r'/\*[\s\S]*?\*/|//.*', '', raw_data)

    df :dict = json.loads(raw_data_comments_ignored)
    modifiers :dict = df["modifier"]
    hotkeys :dict = df["hotkeys"]


    # *---Modifierの設定
    for key, value in modifiers.items():
        keymap.defineModifier(key, value)


    # *---hotkeyの設定
    keymap_global = keymap.defineWindowKeymap()

    def blind(key :str, value :str) -> None:
        default_modifiers_list = ("", "S-", "C-", "C-S-", "A-", "A-S-", "A-C-", "A-C-S-", "W-", "W-S-", "W-C-", "W-C-S-", "W-A-", "W-A-S-", "W-A-C-", "W-A-C-S-")
        for modifier in default_modifiers_list:
            keymap_global[modifier + key] = modifier + value

    def send(key: str, value: str) -> None:
        keymap_global[key] = value

    for genre in hotkeys.keys():
        for key, value in hotkeys[genre].items():
            if genre == "IME":
                send(key, value)
            else:
                blind(key, value)

anonymous set_logger and kill_logger
Python
#------------------------------------------------------------------------------------------
def set_logger():
    #======================================================================================
    # set logger
    #--------------------------------------------------------------------------------------
    # ロガーにファイルハンドラ・コンソールハンドラをセットします
    # プログラムを終了する場合は、ロガーの終了処理(関数:kill_logger)を挟んでください
    # INPUT:なし
    # OUTPUT:なし
    #======================================================================================

    #ロガーを定義
    logger = logging.getLogger(__name__)
    logger.setLevel(DEBUG)

    #出力フォーマット定義
    #20220212 12:00:00 CRITICAL >> Excution has failed.
    formatter = logging.Formatter('[%(asctime)s] %(levelname)s >> %(message)s')

    #コンソール出力用ハンドラ定義
    s_handler = logging.StreamHandler()
    s_handler.setLevel(DEBUG)
    s_handler.setFormatter(formatter)

    #ログファイル出力用ハンドラ定義
    f_handler = logging.FileHandler('.\\ExcelKeyReadLog.log','a+')
    f_handler.setLevel(INFO)
    f_handler.setFormatter(formatter)

    #ハンドラをロガーに登録
    logger.addHandler(s_handler)
    logger.addHandler(f_handler)

    logger.debug('set the logger has completed.')

    #デバッグ:ハンドラ一覧を表示
    print([logging.getLogger(name) for name in logging.root.manager.loggerDict])

    return
#------------------------------------------------------------------------------------------

def kill_logger():

    #======================================================================================
    # kill logger
    #--------------------------------------------------------------------------------------
    # set loggerで起動したロガーを終了します
    # INPUT:なし
    # OUTPUT:なし
    #======================================================================================

    #ロガーを削除
    del logging.Logger.manager.loggerDict[__name__]

    #デバッグ:ハンドラ一覧を表示
    print([logging.getLogger(name) for name in logging.root.manager.loggerDict])
   
    return
#------------------------------------------------------------------------------------------
anonymous No title
Python
import pygame.mixer
import schedule
import time
import tkinter as tk  # Tkinterをインポート

fonts = ("", 20)  # フォントの書式とサイズを指定


# 音再生処理
def Sound():
    pygame.mixer.init()  # 初期化
    pygame.mixer.music.load('alerm1.mp3')  # 読み込み
    pygame.mixer.music.play(-1)  # ループ再生(引数を1にすると1回のみ再生)
    input()
    pygame.mixer.music.stop()  # 終了


def Alarm():
    print("時間です")
    print("\007")  # ビープ音
    Sound()
    exit()  # これがないと無限ループになるので注意


class Calculator(tk.Frame):
    # ウィンドウの作成
    def __init__(self, master=None):
        tk.Frame.__init__(self, master)
        self.master.geometry('1024x768')
        self.master.title('アラーム時刻設定')  # ウィンドウタイトルを指定
        self.entry = tk.Entry(self.master, justify='right', font=fonts)  # テキストボックスを作成
        self.creat_widgets()

    def input(self, num):
        def n():
            self.entry.insert(tk.END, num)
            return n

    def schedule_check(self):  # scheduleモジュールにて定期的にタスク状態を確認
        schedule.run_pending()
        self.after(1, self.schedule_check)

    def btn_click(self):
        self.ala = self.entry.get()
        target = ':'
        idx = self.ala.find(target)
        r = self.ala[idx + 1:]
        s = self.ala[:idx]
        setal = f"{s.zfill(2)}:{r.zfill(2)}"
        schedule.every().day.at(setal).do(Alarm)
        print(f"アラーム時刻・・・{setal}")
        self.schedule_check()

    # 1文字削除
    def clear_one(self):
        text = self.entry.get()
        self.entry.delete(0, tk.END)
        self.entry.insert(0, text[:-1])

    def creat_widgets(self):
        Buttons = [
            ['1', '2', '3', '4', '5'],
            ['6', '7', '8', '9', '0']]

        # 各ボタンを配置
        for i, ro in enumerate(Buttons):
            for j, co in enumerate(ro):
                tk.Button(self.master, font=fonts, text=co, width=7, height=2, command=self.input(co)).grid(row=i + 2,
                                                                                                            column=j)

        # テキストボックスを配置
        self.entry.grid(row=0, column=1, ipady=10, columnspan=3, sticky='nsew')

        # 各ボタンの処理
        tk.Button(self.master, text=':', font=fonts, width=1, height=2, command=self.input(':')).grid(row=9, column=0)
        tk.Button(self.master, text='←', font=fonts, width=1, height=2, command=lambda: self.clear_one()).grid(row=9,
                                                                                                           column=1)
        tk.Button(self.master, text='O K', font=fonts, width=20, height=2, command=lambda: self.btn_click()).grid(row=10,
                                                                                                              column=0,
                                                                                                              columnspan=10)


# 実行
calc = Calculator(tk.Tk())
calc.mainloop()
anonymous No title
Python
"><xmp>
ごるだっく 3node_origin
Python
# (1)拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

# (2)時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T(日)(250,150,150,700と値を変えてシミュレーションを行う)
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T (T=250のときはn=2500)
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個

# (3)SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口(人)(東京都1400万人に匹敵するエリアを想定) N=S+I+R=一定
N2 = 100000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数(人)(10,50,100,5と値を変えてシミュレーションを行う)
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数(日)
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/100
beta23=(beta2+beta3)/100
beta31=(beta3+beta1)/100

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率(隔離率)
# 3-2初期値
Im1_0 = 100                # 初期感染者数(人)100人
Is1_0=0
R1_0 = 0                  # 初期回復者数(人)0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数(人)100人
Is2_0=0
R2_0 = 0                  # 初期回復者数(人)0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数(人)100人
Is3_0=0
R3_0 = 0                  # 初期回復者数(人)0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]=Is2_0
R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0

# (4)数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1](I[j],R[j]も同じ) 要素数n個

  kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

  kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs22 =h * dIs2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR22 = h * dR2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS32 = h * dS3dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm32 =h * dIm3dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm1
ごるだっく 3node_case2検証
Python
# (1)拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化

# (2)時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T(日)(250,150,150,700と値を変えてシミュレーションを行う)
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T (T=250のときはn=2500)
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個

# (3)SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口(人)(東京都1400万人に匹敵するエリアを想定) N=S+I+R=一定
N2 = 5000000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数(人)(10,50,100,5と値を変えてシミュレーションを行う)
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数(日)
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/10
beta23=(beta2+beta3)/10
beta31=(beta3+beta1)/10

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率(隔離率)
# 3-2初期値
Im1_0 = 100                # 初期感染者数(人)100人
Is1_0=0
R1_0 = 0                  # 初期回復者数(人)0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数(人)100人
Is2_0=0
R2_0 = 0                  # 初期回復者数(人)0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数(人)100人
Is3_0=0
R3_0 = 0                  # 初期回復者数(人)0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)



# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]=Is2_0
R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0


#c12をきる場合
c12=0
c23=1
c31=1
for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1](I[j],R[j]も同じ) 要素数n個

  kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

  kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIs22 =h * dIs2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kR22 = h * dR2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
  kS32 = h * dS3dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
  kIm32 =h * dIm3dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm
ごるだっく 3node_rockdown_封鎖期間変化
Python
# (1)拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化
tha=np.arange
# (2)時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T(日)(250,150,150,700と値を変えてシミュレーションを行う)
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T (T=250のときはn=2500)
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個
tha=np.arange(0,30000,100)
lima=np.arange(0,300,1)
f=0 #flag
cnt=0
th=1000 #1日の感染者がこの値を超えたらlim日間ロックダウンする
lim=30 #ロックダウンの日数
# (3)SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口(人)(東京都1400万人に匹敵するエリアを想定) N=S+I+R=一定
N2 = 5000000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数(人)(10,50,100,5と値を変えてシミュレーションを行う)
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数(日)
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/100
beta23=(beta2+beta3)/100
beta31=(beta3+beta1)/100

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率(隔離率)
# 3-2初期値
Im1_0 = 100                # 初期感染者数(人)100人
Is1_0=0
R1_0 = 0                  # 初期回復者数(人)0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数(人)100人
Is2_0=0
R2_0 = 0                  # 初期回復者数(人)0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数(人)100人
Is3_0=0
R3_0 = 0                  # 初期回復者数(人)0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)

peak=np.empty(300)
sumkekka=np.empty(300)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]= Is2_0

R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0

# (4)数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for i in range(300):
  #lim=lima[i]
  lim=lima[i]
  Im1_0 = 100                # 初期感染者数(人)100人
  Is1_0=0
  R1_0 = 0                  # 初期回復者数(人)0人
  S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

  Im2_0 = 0                # 初期感染者数(人)100人
  Is2_0=0
  R2_0 = 0                  # 初期回復者数(人)0人
  S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

  Im3_0 = 0                # 初期感染者数(人)100人
  Is3_0=0
  R3_0 = 0                  # 初期回復者数(人)0人
  S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

  for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1](I[j],R[j]も同じ) 要素数n個

    if f==0 and j>10 and Is1[j]-Is1[j-10]>th:
      f=1
    
    if f==1:
      cnt=cnt+1
    
    if cnt>0 and cnt<lim:
      c23=1
      c12=0
      c31=0
    else:
      c23=1
      c12=1
      c31=1

    kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

    kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/
ごるだっく 3node_rockdown_閾値変化
Python
# (1)拡張モジュールのインポート
import numpy as np                  # 配列を扱う数値計算ライブラリNumPy
import matplotlib.pyplot as plt     # グラフ描画ライブラリmatplotlib
import japanize_matplotlib          # matplotlibの日本語化
tha=np.arange
# (2)時間変数tの導入
T = 700                  # 変数tの範囲 0≦t<T(日)(250,150,150,700と値を変えてシミュレーションを行う)
n = 10*T                 # 変数tの範囲をn等分   n=T/h=T/0.1=10*T (T=250のときはn=2500)
h = 0.1                  # 等差数列の公差:0.1 固定
t = np.arange(0,T,h)     # 0から公差dtでTを超えない範囲で等差数列を生成 t[0],...,t[n-1] 要素数n個
tha=np.arange(0,30000,100)
lima=np.arange(0,300,1)
f=0 #flag
cnt=0
th=1000 #1日の感染者がこの値を超えたらlim日間ロックダウンする
lim=30 #ロックダウンの日数
# (3)SIRモデル
# 3-1パラメータ
N1 = 10000000             # モデルエリアの人口(人)(東京都1400万人に匹敵するエリアを想定) N=S+I+R=一定
N2 = 5000000
N3 = 1000000
m1 = 10
m2 = 10                  # 1日1人あたり接触する人数(人)(10,50,100,5と値を変えてシミュレーションを行う)
m3 = 10
p = 0.02               #5接触ごとに感染が生じる1日あたりの確率
d = 14                   # 感染者の回復平均日数(日)
nu1 = 0                 #ワクチン接種率
nu2 = 0
nu3 = 0
alpha1 = 0.01               #重症化率
alpha2 = 0.01
alpha3 = 0.01
beta1 = m1*p / N1           # 接触あたりの感染率
beta2 = m2*p / N2           # 接触あたりの感染率
beta3 = m3*p / N3           # 接触あたりの感染率
beta12=(beta1+beta2)/100
beta23=(beta2+beta3)/100
beta31=(beta3+beta1)/100

c12=1
c23=1
c31=1

gamma = 1/d              # 回復率(隔離率)
# 3-2初期値
Im1_0 = 100                # 初期感染者数(人)100人
Is1_0=0
R1_0 = 0                  # 初期回復者数(人)0人
S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im2_0 = 0                # 初期感染者数(人)100人
Is2_0=0
R2_0 = 0                  # 初期回復者数(人)0人
S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

Im3_0 = 0                # 初期感染者数(人)100人
Is3_0=0
R3_0 = 0                  # 初期回復者数(人)0人
S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

# 3-3微分方程式
dS1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta1*S1*Im1 - nu1*S1-c12*beta12* S1*Im2-c31*beta31* S1*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta1*S1*Im1 - alpha*Im1 - gamma*Im1+c12*beta12* S1*Im2+c31*beta31* S1*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im1-gamma*Is1
dR1dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im1+Is1)+nu1*S1                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta2*S2*Im2 - nu2*S2-c12*beta12* S2*Im1-c23*beta23* S2*Im3             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta2*S2*Im2 - alpha*Im2 - gamma*Im2+c12*beta12* S2*Im1+c23*beta23* S2*Im3       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im2-gamma*Is2
dR2dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im2+Is2)+nu2*S2                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

dS3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : - beta3*S3*Im3 - nu3*S3 - c31*beta31* S3*Im1-c23*beta23* S3*Im2             # dSdt ≡ dS/dt  dSdt(S, I, R, t)
dIm3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: beta3*S3*Im3 - alpha*Im3 - gamma*Im3+c31*beta31* S3*Im1+c23*beta23* S3*Im2       # dIdt ≡ dI/dt  dIdt(S, I, R, t)
dIs3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t : alpha*Im3-gamma*Is3
dR3dt = lambda S1,S2,S3, Im1,Im2,Im3,Is1,Is2,Is3, R1,R2,R3,alpha,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3,c12,c23,c31, t: gamma*(Im3+Is3)+nu3*S3                  # dRdt ≡ dR/dt  dRdt(S ,I ,R, t)

# 3-4数値積分変数S,I,Rをリストとして生成
S1 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im1 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is1=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R1 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S2 = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im2 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is2=np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R2 = np.empty(n)          # R[0],...,R[n-1] 要素数n個

S3  = np.empty(n)          # S[0],...,S[n-1] 要素数n個
Im3 = np.empty(n)          # Im[0],...,Im[n-1] 要素数n個
Is3 = np.empty(n)            # Is[0],...,Is[n-1] 要素数n個
R3  = np.empty(n)          # R[0],...,R[n-1] 要素数n個

sum1=np.empty(n)
sum2=np.empty(n)
sum3=np.empty(n)

peak=np.empty(300)
sumkekka=np.empty(300)
# 3-5初期値代入
S1[0] = S1_0
Im1[0] = Im1_0
Is1[0]=Is1_0
R1[0] = R1_0

S2[0] = S2_0
Im2[0] = Im2_0
Is2[0]= Is2_0

R2[0] = R2_0

S3[0] = S3_0
Im3[0] = Im3_0
Is3[0]=Is3_0
R3[0] = R3_0
sum1[0]=0
sum2[0]=0
sum3[0]=0

# (4)数値積分 4次ルンゲ-クッタ法 4th-Order Runge–Kutta Methods
for i in range(300):
  #lim=lima[i]
  th=tha[i]
  Im1_0 = 100                # 初期感染者数(人)100人
  Is1_0=0
  R1_0 = 0                  # 初期回復者数(人)0人
  S1_0 = N1 - Im1_0 - Is1_0 - R1_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

  Im2_0 = 0                # 初期感染者数(人)100人
  Is2_0=0
  R2_0 = 0                  # 初期回復者数(人)0人
  S2_0 = N2 - Im2_0 - Is2_0 - R2_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

  Im3_0 = 0                # 初期感染者数(人)100人
  Is3_0=0
  R3_0 = 0                  # 初期回復者数(人)0人
  S3_0 = N3 - Im3_0 - Is3_0 - R3_0      # 初期未感染者数(人)S_0 = N - Im_0- Is_0 - R_0

  for j in range(n-1):     # j=0,...,n-2 -> S[j]=S[0],...,S[n-1](I[j],R[j]も同じ) 要素数n個

    if f==0 and j>10 and Is1[j]-Is1[j-10]>th:
      f=1
    
    if f==1:
      cnt=cnt+1
    
    if cnt>0 and cnt<lim:
      c23=1
      c12=0
      c31=0
    else:
      c23=1
      c12=1
      c31=1

    kS11 = h * dS1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm11 = h * dIm1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs11 = h * dIs1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR11 = h * dR1dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS21 = h * dS2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm21 = h * dIm2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs21 = h * dIs2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR21 = h * dR2dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS31 = h * dS3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm31 = h * dIm3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs31 = h * dIs3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR31 = h * dR3dt( S1[j],S2[j],S3[j] ,Im1[j],Im2[j],Im3[j],Is1[j],Is2[j],Is3[j],R1[j] ,R2[j],R3[j],alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])

    kS12 = h * dS1dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm12 =h * dIm1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIs12 =h * dIs1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kR12 = h * dR1dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j])
    kS22 = h * dS2dt(  S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,Is1[j]+kIs11/2,Is2[j]+kIs21/2,Is3[j]+kIs31/2,R1[j]+kR11/2 ,R2[j]+kR21/2,R3[j]+kR31/2,alpha1,beta1,beta2,beta3,beta12,beta23,beta31,nu1,nu2,nu3 ,c12,c23,c31,t[j] )
    kIm22 =h * dIm2dt(S1[j]+kS11/2,S2[j]+kS21/2,S3[j]+kS31/2 ,Im1[j]+kIm11/2,Im2[j]+kIm21/2,Im3[j]+kIm31/2,
Don't you submit code?
Submit