Игровая модель поведения на рынке двух конкурирующих фирм на Python

# -*- coding: utf8 -*-    
import numpy as np
import matplotlib.pyplot as plt
a=10
def f(q): Функция цены от объёма товара
         return a*np.e**(-0.5*q**2)
def h(q): #Функция издержек от объёма товара
         return np.sqrt(q)
plt.figure()
q= np.arange(0, 2.01, 0.1)#Массив значений аргумента
plt.title(r'$y=f(q)$') #Заголовок в формате TeX
plt.ylabel(r'$f(q)$') #Метка по оси y в формате TeX
plt.xlabel(r'$q$') #Метка по оси x в формате TeX
plt.grid(True) #Сетка
plt.plot(q,f(q)) #Построение графика
plt.figure()
plt.title(r'$y=h(q)$') #Заголовок в формате TeX
plt.ylabel(r'$h(q)$') #Метка по оси y в формате TeX
plt.xlabel(r'$q$') #Метка по оси x в формате TeX
plt.grid(True) #Сетка
plt.plot(q,h(q)) #Построение графика
plt.show() #Показать график

# -*- coding: utf8 -*-
import pylab
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
xmax=1;ymax=1;N=20; h1=xmax/N;h2=ymax/N# Параметры модели
def L1(x,y):# Функция прибыли первой фирмы
         return 10*x*np.e**(-0.5*(x+y)**2)-np.sqrt(x)
def L2(x,y):# Функция прибыли второй фирмы
        return  10*y*np.e**(-0.5*(x+y)**2)-np.sqrt(y)
def fi(a,b):# Функция определения индекса списка
         for w in range(0,len(a)+1):
                  if a[w]==b:
                           return w         
L= np.zeros([N+1,N+1])# Создания матрицы прибылей
for i in range(0,N+1):
         for j in range(0,N+1):
             L[i,j]=L1(i*h1,j*h2)  
A=[];e=[]
rows, cols = L.shape
for i in range(rows):
         e=[]
         for j in range(cols):
                  e.append(L[i,j])               
         A.append(min(e))
a=max(A);ia=fi(A,a);xa=ia*h1# Итоги работы первой фирмы
B=[];e=[]
rows, cols = L.shape
for j in range(cols):
         e=[]
         for i in range(rows):
                  e.append(L[i,j])               
         B.append(max(e))        
b=min(B);ib=fi(B,b);yb=ib*h2# Итоги работы второй фирмы
p1=round(L1(xa,yb),3)
print("Прибыль первой фирмы -"+str(p1))
p2=round(L2(xa,yb),3)
print("Прибыль второй фирмы -"+str(p2))
def makeData_L1 ():# Поверхность прибыль- z, объёмы закупок - x,y для первой фирмы
         x=[h1*i for i in np.arange(0,N+1)]
         y=[h2*i for i in np.arange(0,N+1)]
         x,y= np.meshgrid(x, y)
         z=[]
         for i in range(0,N+1):
                  z.append(L1(x[i],y[i]))
         return x, y, z
fig = pylab.figure()
axes = Axes3D(fig)
x, y, z = makeData_L1()
axes.plot_surface(x, y, z)
def makeData_L2 ():# Поверхность прибыль- z объёмы закупок- x,y для второй фирмы
         x=[h1*i for i in np.arange(0,N+1)]
         y=[h2*i for i in np.arange(0,N+1)]
         x,y= np.meshgrid(x, y)
         z=[]
         for i in range(0,N+1):
                  z.append(L2(x[i],y[i]))
         return x, y, z
fig = pylab.figure()
axes = Axes3D(fig)
x, y, z = makeData_L2()
axes.plot_surface(x, y, z)
pylab.show()

# -*- coding: utf8 -*-
import pylab
import numpy as np
xmax=1;ymax=1;N=20; h1=xmax/N;h2=ymax/N
def L1(x,y):
         return 10*x*np.e**(-0.5*(x+y)**2)-np.sqrt(x)
print("Прибыль каждого при стратегии В -"+str(round(L1(xmax,ymax),3)))

# -*- coding: utf8 -*-
import matplotlib.pyplot as plt
import numpy as np
xmax=1;ymax=1;N=20; h1=xmax/N;h2=ymax/N
def L1(x,y):
         return 10*x*np.e**(-0.5*(x+y)**2)-np.sqrt(x)
z=[L1(i*h1,i*h1)for i in np.arange(0,N+1)]
x=[h1*i for i in np.arange(0,N+1)]
xopt=x[z.index(max(z))]
print("Объём товара каждой фирмы при котором достигается -"+str(xopt))
print("Прибыль каждой фирмы -"+str(round(max(z),3)))
plt.title('Зависимость прибыли от объёма товара для стратегии С')
plt.grid(True)
plt.plot(x,z) 
plt.show()