Не запускается код который был написан на python 2.7

Всем добрый день помогите запустить этот код(2d моделирование), препод дал задание сделать gui, но после запуска программы выходит такая формавведите сюда описание изображения

но чтоб началась само моделирование надо закрыть форму и вот тут начинается вся беда, начинается открываться несколько форм и выходит такая ошибка:

RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (matplotlib.pyplot.figure) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam figure.max_open_warning). fig = plt.figure()

line 309, in plt.show(1)

line 353, in show return _backend_mod.show(*args, **kwargs) TypeError: show() takes 1 positional argument but 2 were given

Пожалуйста помогите решить эту проблему, так как в питоне я вообще новичок и не понимаю, как решить эту проблему. У самого питон 3.9 и пишу код в pycharm'е. У препода все работает, но у него питон 2.7 стоит. Возможно ли вообще запустить такой код на 3.9?

#import idlelib.PyShell
#import sys
import time
import pylab
import numpy as np
#import sys, idlelib.PyShell; idlelib.PyShell.warning_stream = sys.stderr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#idlelib.PyShell.warning_stream = sys.stderr

N1 = 40
tt =3.0
l = 2.00
h = l/N1
N2 = 50
t=tt/N2
t2=t/2.0
hh=h*h
#rr=hh/t2
rr=hh/t

a11=0.03
a12=0.01
p1=-0.05
a21=0.03
a22=0.01
p2=-0.03
q1=0.0
q2=0.0

m1=0.95
m2=1.2
d1=0.7
d2=0.5
c1=0.25
c2=0.65
h1=0.55
b1=0.1

def e1(u,v):
    return a11+a12*u*v
def e2(u,v):
    return a21+a22*u*v

def f1(u,v):
    return u*(m1-d1)+v*d1-c1*u*u-h1*u-b1
def f2(u,v):
    return v*(m2-d2)+u*d2-c2*v*v

def Rs(x,y):
    return 4
    #return 2.*np.sin(np.pi*y/3.)+0.3*np.sin(3*np.pi*y/3.)*np.sin(2.*np.pi*y/3.)+np.sin(np.pi*x/2.)
    #return 2.*np.sin((np.pi*y/3.)+0.2)+0.8*np.sin((3.*np.pi*y/3.)+1.0)*np.sin((1.*np.pi*y/3.)+2.0)+np.sin(np.pi*x/2.)
    #return 2.*np.sin((np.pi*y/3.)+0.5)+0.8*np.sin((3.*np.pi*y/3.))*np.sin((1.*np.pi*y/3.)+2.)+np.sin(np.pi*x/2.)
    #return np.exp((y/2)*np.sin(3.*np.pi*y/2.)) * np.sin((1.*np.pi*x/3.)+0.1)

def y_0(x,y):
    return 1.5 * np.exp(-30.0*(((x/1.8)-0.55)**2 + (y-0.35)**2))
def v_0(x,y):
    return 1.5 * np.exp(-15.0*(((x/2.0)-0.5)**2 + (y-1.35)**2))

x1 = np.zeros((N1+1), 'float')
x2 = np.zeros((N1+1), 'float')
y = np.zeros((N1+1, N1+1), 'float')
y0 = np.zeros((N1+1, N1+1), 'float')
v = np.zeros((N1+1, N1+1), 'float')
v0 = np.zeros((N1+1, N1+1), 'float')
a = np.zeros((N1+1), 'float')
b = np.zeros((N1+1), 'float')
z = np.zeros((N1+1), 'float')
R = np.zeros((N1+1, N1+1), 'float')

#Computation of y,v
tst = time.clock()

for i1 in range(0,N1+1):
    x1[i1]=i1*h
for i2 in range(0,N1+1):
    x2[i2]=i2*h

#for i2 in range(0,N1+1):
    #for i1 in range(0,N1+1):
        #if i2*h<0.62: y[i1,i2] = 0.0
        #elif i2*h<1.27: y[i1,i2]=-21.30*(i2*h)*(i2*h)+40.26*(i2*h)-16.77
        #else: y[i1,i2]=0.0
        #if i2*h<1.32: v[i1,i2] = 0.0
        #elif i2*h<1.96: v[i1,i2]=-21.30*(i2*h)*(i2*h)+69.65*(i2*h)-54.69
        #else: v[i1,i2]=0.0

for i1 in range(0, N1+1):
    for i2 in range(0, N1+1):
        y0[i1, i2] = y_0(x1[i1], x2[i2])
        v0[i1, i2] = v_0(x1[i1], x2[i2])
        R[i1, i2] = Rs(x1[i1], x2[i2])
        y[i1, i2] = y0[i1, i2]
        v[i1, i2] = v0[i1, i2]

fig = plt.figure()
#ax = Axes3D(fig)
x1 = np.linspace(0., l, N1+1)
x2 = np.linspace(0., l, N1+1)
X1,X2 = np.meshgrid(x1, x2)
plt.xlabel('X')
plt.ylabel('Y')
#ax.plot_surface(X1, X2, y0, rstride=1, cstride=1, cmap='Blues')
#ax.plot_surface(X1, X2, v0, rstride=1, cstride=1, cmap='BuGn')
#ax.plot_wireframe(X1, X2, R, rstride=1, cstride=1)
#ax.plot_surface(X1, X2, R, rstride=1, cstride=1, cmap='summer')

          #Liniya urovnya ('ax' turn off before activation)

g = np.linspace(0., 2., 21)
#plt.clabel(plt.contour(x1,x2,y0,g,colors='gray'), inline=True, fontsize=8)
#plt.contourf(x1,x2,y0,g,cmap=plt.cm.Blues)
plt.clabel(plt.contour(x1,x2,v0,g,colors='gray'), inline=True, fontsize=8)
plt.contourf(x1,x2,v0,g,cmap=plt.cm.BuGn)
plt.colorbar()

plt.show()
pylab.ion()

for j in range(1,N2+1):

# y[k+1/2] to x1
    for i2 in range(0,N1+1):
        #a[0]=0.
        #b[0]=0.
        e210 = e1((y[1,i2]+y[0,i2])/2.00,(v[1,i2]+v[0,i2])/2.00)
        a[0] = e210/(e210+rr/2.0)
        FD0 = (-p1/rr)*(y[1,i2]+y[0,i2])*(v[1,i2]-v[0,i2])- \
              (q1/rr)*(y[1,i2]+y[0,i2])*(R[1,i2]-R[0,i2])+f1(y[0,i2],v[0,i2])*t2
        b[0] = (y[0,i2]+FD0)/(1+2.0*e210/rr)

        for i1 in range(1,N1):
            e11 = e1((y[i1,i2]+y[i1-1,i2])/2.00,(v[i1,i2]+v[i1-1,i2])/2.00)
            e12 = e1((y[i1,i2]+y[i1+1,i2])/2.00,(v[i1,i2]+v[i1+1,i2])/2.00)
            aa = e11/rr
            bb = e12/rr
            cc = aa+bb+1
            a[i1] = bb/(cc-aa*a[i1-1])

            r11 = (1.0/rr)*(y[i1,i2]+y[i1-1,i2])/2.00
            r12 = (1.0/rr)*(y[i1,i2]+y[i1+1,i2])/2.00

            L1v = p1*(r12*(v[i1+1,i2]-v[i1,i2]) - r11*(v[i1,i2]-v[i1-1,i2]))
            L1R = q1*(r12*(R[i1+1,i2]-R[i1,i2]) - r11*(R[i1,i2]-R[i1-1,i2]))

            d = y[i1,i2] + f1(y[i1,i2],v[i1,i2])*t2 - L1v - L1R
            b[i1]=(aa*b[i1-1]+d)/(cc-aa*a[i1-1])
        #y[N1,i2]=b[N1-1]/(1-a[N1-1])
        #y[N1,i2]=0.
        e22n = e1((y[N1,i2]+y[N1-1,i2])/2.00,(v[N1,i2]+v[N1-1,i2])/2.00)
        FDn = p1*(y[N1,i2]+y[N1-1,i2])*(v[N1,i2]-v[N1-1,i2])/2.0 + \
              q1*(y[N1,i2]+y[N1-1,i2])*(R[N1,i2]-R[N1-1,i2])/2.0 + \
              f1(y[N1,i2],v[N1,i2])*hh/4.0
        y[N1,i2]=(e22n*b[N1-1] + rr*y[N1,i2]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)
        for ii in range(1,N1+1):
            i1=N1-ii
            y[i1,i2] = a[i1]*y[i1+1,i2]+b[i1]

# v[k+1/2] to x1
        e210 = e1((y[1,i2]+y[0,i2])/2.00,(v[1,i2]+v[0,i2])/2.00)
        a[0] = e210/(e210+rr/2.0)
        FD0 = (-p2/rr)*(v[1,i2]+v[0,i2])*(y[1,i2]-y[0,i2])- \
              (q2/rr)*(v[1,i2]+v[0,i2])*(R[1,i2]-R[0,i2])+f2(y[0,i2],v[0,i2])*t2
        b[0] = (v[0,i2]+FD0)/(1+2.0*e210/rr)
        for i1 in range(1,N1):
            e21=e2((y[i1,i2]+y[i1-1,i2])/2.00,(v[i1,i2]+v[i1-1,i2])/2.00)
            e22=e2((y[i1,i2]+y[i1+1,i2])/2.00,(v[i1,i2]+v[i1+1,i2])/2.00)
            aa = e21/rr
            bb = e22/rr
            cc = aa+bb+1
            a[i1] = bb/(cc-aa*a[i1-1])

            r11 = (1.0/rr)*(v[i1,i2]+v[i1-1,i2])/2.
            r12 = (1.0/rr)*(v[i1,i2]+v[i1+1,i2])/2.

            L1y = p2*(r12*(y[i1+1,i2]-y[i1,i2]) - r11*(y[i1,i2]-y[i1-1,i2]))
            L1R = q2*(r12*(R[i1+1,i2]-R[i1,i2]) - r11*(R[i1,i2]-R[i1-1,i2]))

            d = v[i1,i2] + f2(y[i1,i2],v[i1,i2])*t2 - L1y - L1R
            b[i1]=(aa*b[i1-1]+d)/(cc-aa*a[i1-1])
        #v[N1,i2]=b[N1-1]/(1-a[N1-1])
        #v[N1,i2]=0.
        e22n = e2((y[N1,i2]+y[N1-1,i2])/2.00,(v[N1,i2]+v[N1-1,i2])/2.00)
        FDn = p2*(v[N1,i2]+v[N1-1,i2])*(y[N1,i2]-y[N1-1,i2])/2.0 + \
              q2*(v[N1,i2]+v[N1-1,i2])*(R[N1,i2]-R[N1-1,i2])/2.0 + \
              f2(y[N1,i2],v[N1,i2])*hh/4.0
        v[N1,i2]=(e22n*b[N1-1] + rr*v[N1,i2]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)

        for ii in range(1,N1+1):
            i1=N1-ii
            v[i1,i2] = a[i1]*v[i1+1,i2]+b[i1]

# y[k+1] to x2
    for i1 in range(0,N1+1):
        a[0]=0.
        b[0]=0.
        #e210 = e1((y[i1,1]+y[i1,0])/2.00,(v[i1,1]+v[i1,0])/2.00)
        #FD0 = (-p1/rr)*(y[i1,1]+y[i1,0])*(v[i1,1]-v[i1,0])- \
        #      (q1/rr)*(y[i1,1]+y[i1,0])*(R[i1,1]-R[i1,0])+f1(y[i1,0],v[i1,0])*t2
        #b[0] = (y[i1,0]+FD0)/(1+2.0*e210/rr)

        for i2 in range(1,N1):
            e21 = e1((y[i1,i2]+y[i1,i2-1])/2.00,(v[i1,i2]+v[i1,i2-1])/2.00)
            e22 = e1((y[i1,i2]+y[i1,i2+1])/2.00,(v[i1,i2]+v[i1,i2+1])/2.00)
            aa = e21/rr
            bb = e22/rr
            cc = aa+bb+1
            a[i2] = bb/(cc-aa*a[i2-1])

            r21 = (1.0/rr)*(y[i1,i2]+y[i1,i2-1])/2.
            r22 = (1.0/rr)*(y[i1,i2]+y[i1,i2+1])/2.

            L2v = p1*(r22*(v[i1,i2+1]-v[i1,i2]) - r21*(v[i1,i2]-v[i1,i2-1]))
            L2R = q1*(r22*(R[i1,i2+1]-R[i1,i2]) - r21*(R[i1,i2]-R[i1,i2-1]))

            d = y[i1,i2] + f1(y[i1,i2],v[i1,i2])*t2 - L2v - L2R
            b[i2]=(aa*b[i2-1]+d)/(cc-aa*a[i2-1])
        y[i1,N1]=0.
        #e22n = e1((y[i1,N1]+y[i1,N1-1])/2.00,(v[i1,N1]+v[i1,N1-1])/2.00)
        #FDn = p1*(y[i1,N1]+y[i1,N1-1])*(v[i1,N1]-v[i1,N1-1])/2.0 + \
        #      q1*(y[i1,N1]+y[i1,N1-1])*(R[i1,N1]-R[i1,N1-1])/2.0 + \
        #      f1(y[i1,N1],v[i1,N1])*hh/4.0
        #y[i1,N1]=(e22n*b[N1-1] + rr*y[i1,N1]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)

        for ii in range(1,N1+1):
            i2=N1-ii
            y[i1,i2] = a[i2]*y[i1,i2+1]+b[i2]

        for i2 in range(0,N1+1):
            if y[i1,i2]<0: y[i1,i2]=0.0

# v[k+1] to x2
        #e210 = e2((y[i1,1]+y[i1,0])/2.00,(v[i1,1]+v[i1,0])/2.00)
        #a[0] = e210/(e210+rr/2.0)
        #FD0 = (-p2/rr)*(v[i1,1]+v[i1,0])*(y[i1,1]-y[i1,0])- \
        #      (q2/rr)*(v[i1,1]+v[i1,0])*(R[i1,1]-R[i1,0])+f2(y[i1,0],v[i1,0])*t2
        #b[0] = (v[i1,0]+FD0)/(1+2.0*e210/rr)

        for i2 in range(1,N1):
            e21=e2((y[i1,i2]+y[i1,i2-1])/2.00,(v[i1,i2]+v[i1,i2-1])/2.00)
            e22=e2((y[i1,i2]+y[i1,i2+1])/2.00,(v[i1,i2]+v[i1,i2+1])/2.00)
            aa = e21/rr
            bb = e22/rr
            cc = aa+bb+1
            a[i2] = bb/(cc-aa*a[i2-1])

            r21 = (1.0/rr)*(v[i1,i2]+v[i1,i2-1])/2.
            r22 = (1.0/rr)*(v[i1,i2]+v[i1,i2+1])/2.

            L2y = p2*(r22*(y[i1,i2+1]-y[i1,i2])-r21*(y[i1,i2]-y[i1,i2-1]))
            L2R = q2*(r22*(R[i1,i2+1]-R[i1,i2])-r21*(R[i1,i2]-R[i1,i2-1]))

            d = v[i1,i2] + f2(y[i1,i2],v[i1,i2])*t2 - L2y - L2R
            b[i2]=(aa*b[i2-1]+d)/(cc-aa*a[i2-1])
        v[i1,N1]=0.
        #e22n = e2((y[i1,N1]+y[i1,N1-1])/2.00,(v[i1,N1]+v[i1,N1-1])/2.00)
        #FDn = p2*(v[i1,N1]+v[i1,N1-1])*(y[i1,N1]-y[i1,N1-1])/2.0 + \
        #      q2*(v[i1,N1]+v[i1,N1-1])*(R[i1,N1]-R[i1,N1-1])/2.0 + \
        #      f2(y[i1,N1],v[i1,N1])*hh/4.0
        #v[i1,N1]=(e22n*b[N1-1] + rr*v[i1,N1]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)

        for ii in range(1,N1+1):
            i2=N1-ii
            v[i1,i2] = a[i2]*v[i1,i2+1]+b[i2]

        for i2 in range(0,N1+1):
            if v[i1,i2]<0: v[i1,i2]=0.0
#graphics
    pylab.clf()
    x1 = np.linspace(0., l, N1+1)
    x2 = np.linspace(0., l, N1+1)
    fig = plt.figure()
    ax = Axes3D(fig)
    X1, X2 = np.meshgrid(x1, x2)
    plt.xlabel('X')
    plt.ylabel('Y')
    ax.plot_surface(X1, X2, y, rstride=1, cstride=1, cmap='Blues')
    ax.plot_surface(X1, X2, v, rstride=1, cstride=1, cmap='BuGn')
    #ax.plot_wireframe(X1, X2, y, rstride=1, cstride=1)
    #plt.contourf(XX1, XX2, y, cmap=plt.cm.gray)
    #plt.colorbar()
    #plt.show()
    pylab.draw()

fig = plt.figure()
#ax = Axes3D(fig)
x1 = np.linspace(0., l, N1+1)
x2 = np.linspace(0., l, N1+1)
X1,X2 = np.meshgrid(x1, x2)
plt.xlabel('X')
plt.ylabel('Y')
#ax.plot_surface(X1, X2, y, rstride=1, cstride=1, cmap='Blues')
#ax.plot_surface(X1, X2, v, rstride=1, cstride=1, cmap='BuGn')

          #Liniya urovnya ('ax' turn off before activation)

g = np.linspace(0., 2.0, 21)
plt.clabel(plt.contour(x1,x2,y,g,colors='gray'), inline=True, fontsize=8)
plt.contourf(x1,x2,y,g,cmap=plt.cm.Blues)
#plt.clabel(plt.contour(x1,x2,v,g,colors='gray'), inline=True, fontsize=8)
#plt.contourf(x1,x2,v,g,cmap=plt.cm.BuGn)
plt.colorbar()

plt.show()
plt.show(1)
f1 = open('rezy.txt','w')
for i in range(0,N1+1):
    for j in range(0,N1+1):
        f1.write("%f\n" % y[i,j])
f1.close()
f2 = open('rezv.txt','w')
for i in range(0,N1+1):
    for j in range(0,N1+1):
        f2.write("%f\n" % v[i,j])
    f2.write("\n") 
f2.close()

dt = time.clock() - tst
print "N2 = %i, N1 = %i, time solution = %1.3e, h = %1.3e" % (N2, N1, dt, h)






Ответы (1 шт):

Автор решения: GrAnd

Вот, поковырял чутка. Какую-то анимацию роста горбиков показывает. Python 3.7.4 под Windows 7. И Python 3.9.0 под Windows 10.

#import idlelib.PyShell
#import sys
import time
#import pylab
import numpy as np
#import sys, idlelib.PyShell; idlelib.PyShell.warning_stream = sys.stderr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#idlelib.PyShell.warning_stream = sys.stderr

N1 = 40
tt =3.0
l = 2.00
h = l/N1
N2 = 50
t=tt/N2
t2=t/2.0
hh=h*h
#rr=hh/t2
rr=hh/t

a11=0.03
a12=0.01
p1=-0.05
a21=0.03
a22=0.01
p2=-0.03
q1=0.0
q2=0.0

m1=0.95
m2=1.2
d1=0.7
d2=0.5
c1=0.25
c2=0.65
h1=0.55
b1=0.1

def e1(u,v):
    return a11+a12*u*v
def e2(u,v):
    return a21+a22*u*v

def f1(u,v):
    return u*(m1-d1)+v*d1-c1*u*u-h1*u-b1
def f2(u,v):
    return v*(m2-d2)+u*d2-c2*v*v

def Rs(x,y):
    return 4
    #return 2.*np.sin(np.pi*y/3.)+0.3*np.sin(3*np.pi*y/3.)*np.sin(2.*np.pi*y/3.)+np.sin(np.pi*x/2.)
    #return 2.*np.sin((np.pi*y/3.)+0.2)+0.8*np.sin((3.*np.pi*y/3.)+1.0)*np.sin((1.*np.pi*y/3.)+2.0)+np.sin(np.pi*x/2.)
    #return 2.*np.sin((np.pi*y/3.)+0.5)+0.8*np.sin((3.*np.pi*y/3.))*np.sin((1.*np.pi*y/3.)+2.)+np.sin(np.pi*x/2.)
    #return np.exp((y/2)*np.sin(3.*np.pi*y/2.)) * np.sin((1.*np.pi*x/3.)+0.1)

def y_0(x,y):
    return 1.5 * np.exp(-30.0*(((x/1.8)-0.55)**2 + (y-0.35)**2))
def v_0(x,y):
    return 1.5 * np.exp(-15.0*(((x/2.0)-0.5)**2 + (y-1.35)**2))

x1 = np.zeros((N1+1), 'float')
x2 = np.zeros((N1+1), 'float')
y = np.zeros((N1+1, N1+1), 'float')
y0 = np.zeros((N1+1, N1+1), 'float')
v = np.zeros((N1+1, N1+1), 'float')
v0 = np.zeros((N1+1, N1+1), 'float')
a = np.zeros((N1+1), 'float')
b = np.zeros((N1+1), 'float')
z = np.zeros((N1+1), 'float')
R = np.zeros((N1+1, N1+1), 'float')

#Computation of y,v
tst = time.perf_counter()

for i1 in range(0,N1+1):
    x1[i1]=i1*h
for i2 in range(0,N1+1):
    x2[i2]=i2*h

#for i2 in range(0,N1+1):
    #for i1 in range(0,N1+1):
        #if i2*h<0.62: y[i1,i2] = 0.0
        #elif i2*h<1.27: y[i1,i2]=-21.30*(i2*h)*(i2*h)+40.26*(i2*h)-16.77
        #else: y[i1,i2]=0.0
        #if i2*h<1.32: v[i1,i2] = 0.0
        #elif i2*h<1.96: v[i1,i2]=-21.30*(i2*h)*(i2*h)+69.65*(i2*h)-54.69
        #else: v[i1,i2]=0.0

for i1 in range(0, N1+1):
    for i2 in range(0, N1+1):
        y0[i1, i2] = y_0(x1[i1], x2[i2])
        v0[i1, i2] = v_0(x1[i1], x2[i2])
        R[i1, i2] = Rs(x1[i1], x2[i2])
        y[i1, i2] = y0[i1, i2]
        v[i1, i2] = v0[i1, i2]

fig = plt.figure()
x1 = np.linspace(0., l, N1+1)
x2 = np.linspace(0., l, N1+1)
X1,X2 = np.meshgrid(x1, x2)
plt.xlabel('X')
plt.ylabel('Y')
#ax.plot_surface(X1, X2, y0, rstride=1, cstride=1, cmap='Blues')
#ax.plot_surface(X1, X2, v0, rstride=1, cstride=1, cmap='BuGn')
#ax.plot_wireframe(X1, X2, R, rstride=1, cstride=1)
#ax.plot_surface(X1, X2, R, rstride=1, cstride=1, cmap='summer')

          #Liniya urovnya ('ax' turn off before activation)

g = np.linspace(0., 2., 21)
plt.clabel(plt.contour(x1,x2,y0,g,colors='gray'), inline=True, fontsize=8)
plt.contourf(x1,x2,y0,g,cmap=plt.cm.Blues)
plt.clabel(plt.contour(x1,x2,v0,g,colors='gray'), inline=True, fontsize=8)
plt.contourf(x1,x2,v0,g,cmap=plt.cm.BuGn)
plt.colorbar()
plt.show()
plt.ion()
fig = plt.figure()


for j in range(1,N2+1):

# y[k+1/2] to x1
    for i2 in range(0,N1+1):
        #a[0]=0.
        #b[0]=0.
        e210 = e1((y[1,i2]+y[0,i2])/2.00,(v[1,i2]+v[0,i2])/2.00)
        a[0] = e210/(e210+rr/2.0)
        FD0 = (-p1/rr)*(y[1,i2]+y[0,i2])*(v[1,i2]-v[0,i2])- \
              (q1/rr)*(y[1,i2]+y[0,i2])*(R[1,i2]-R[0,i2])+f1(y[0,i2],v[0,i2])*t2
        b[0] = (y[0,i2]+FD0)/(1+2.0*e210/rr)

        for i1 in range(1,N1):
            e11 = e1((y[i1,i2]+y[i1-1,i2])/2.00,(v[i1,i2]+v[i1-1,i2])/2.00)
            e12 = e1((y[i1,i2]+y[i1+1,i2])/2.00,(v[i1,i2]+v[i1+1,i2])/2.00)
            aa = e11/rr
            bb = e12/rr
            cc = aa+bb+1
            a[i1] = bb/(cc-aa*a[i1-1])

            r11 = (1.0/rr)*(y[i1,i2]+y[i1-1,i2])/2.00
            r12 = (1.0/rr)*(y[i1,i2]+y[i1+1,i2])/2.00

            L1v = p1*(r12*(v[i1+1,i2]-v[i1,i2]) - r11*(v[i1,i2]-v[i1-1,i2]))
            L1R = q1*(r12*(R[i1+1,i2]-R[i1,i2]) - r11*(R[i1,i2]-R[i1-1,i2]))

            d = y[i1,i2] + f1(y[i1,i2],v[i1,i2])*t2 - L1v - L1R
            b[i1]=(aa*b[i1-1]+d)/(cc-aa*a[i1-1])
        #y[N1,i2]=b[N1-1]/(1-a[N1-1])
        #y[N1,i2]=0.
        e22n = e1((y[N1,i2]+y[N1-1,i2])/2.00,(v[N1,i2]+v[N1-1,i2])/2.00)
        FDn = p1*(y[N1,i2]+y[N1-1,i2])*(v[N1,i2]-v[N1-1,i2])/2.0 + \
              q1*(y[N1,i2]+y[N1-1,i2])*(R[N1,i2]-R[N1-1,i2])/2.0 + \
              f1(y[N1,i2],v[N1,i2])*hh/4.0
        y[N1,i2]=(e22n*b[N1-1] + rr*y[N1,i2]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)
        for ii in range(1,N1+1):
            i1=N1-ii
            y[i1,i2] = a[i1]*y[i1+1,i2]+b[i1]

# v[k+1/2] to x1
        e210 = e1((y[1,i2]+y[0,i2])/2.00,(v[1,i2]+v[0,i2])/2.00)
        a[0] = e210/(e210+rr/2.0)
        FD0 = (-p2/rr)*(v[1,i2]+v[0,i2])*(y[1,i2]-y[0,i2])- \
              (q2/rr)*(v[1,i2]+v[0,i2])*(R[1,i2]-R[0,i2])+f2(y[0,i2],v[0,i2])*t2
        b[0] = (v[0,i2]+FD0)/(1+2.0*e210/rr)
        for i1 in range(1,N1):
            e21=e2((y[i1,i2]+y[i1-1,i2])/2.00,(v[i1,i2]+v[i1-1,i2])/2.00)
            e22=e2((y[i1,i2]+y[i1+1,i2])/2.00,(v[i1,i2]+v[i1+1,i2])/2.00)
            aa = e21/rr
            bb = e22/rr
            cc = aa+bb+1
            a[i1] = bb/(cc-aa*a[i1-1])

            r11 = (1.0/rr)*(v[i1,i2]+v[i1-1,i2])/2.
            r12 = (1.0/rr)*(v[i1,i2]+v[i1+1,i2])/2.

            L1y = p2*(r12*(y[i1+1,i2]-y[i1,i2]) - r11*(y[i1,i2]-y[i1-1,i2]))
            L1R = q2*(r12*(R[i1+1,i2]-R[i1,i2]) - r11*(R[i1,i2]-R[i1-1,i2]))

            d = v[i1,i2] + f2(y[i1,i2],v[i1,i2])*t2 - L1y - L1R
            b[i1]=(aa*b[i1-1]+d)/(cc-aa*a[i1-1])
        #v[N1,i2]=b[N1-1]/(1-a[N1-1])
        #v[N1,i2]=0.
        e22n = e2((y[N1,i2]+y[N1-1,i2])/2.00,(v[N1,i2]+v[N1-1,i2])/2.00)
        FDn = p2*(v[N1,i2]+v[N1-1,i2])*(y[N1,i2]-y[N1-1,i2])/2.0 + \
              q2*(v[N1,i2]+v[N1-1,i2])*(R[N1,i2]-R[N1-1,i2])/2.0 + \
              f2(y[N1,i2],v[N1,i2])*hh/4.0
        v[N1,i2]=(e22n*b[N1-1] + rr*v[N1,i2]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)

        for ii in range(1,N1+1):
            i1=N1-ii
            v[i1,i2] = a[i1]*v[i1+1,i2]+b[i1]

# y[k+1] to x2
    for i1 in range(0,N1+1):
        a[0]=0.
        b[0]=0.
        #e210 = e1((y[i1,1]+y[i1,0])/2.00,(v[i1,1]+v[i1,0])/2.00)
        #FD0 = (-p1/rr)*(y[i1,1]+y[i1,0])*(v[i1,1]-v[i1,0])- \
        #      (q1/rr)*(y[i1,1]+y[i1,0])*(R[i1,1]-R[i1,0])+f1(y[i1,0],v[i1,0])*t2
        #b[0] = (y[i1,0]+FD0)/(1+2.0*e210/rr)

        for i2 in range(1,N1):
            e21 = e1((y[i1,i2]+y[i1,i2-1])/2.00,(v[i1,i2]+v[i1,i2-1])/2.00)
            e22 = e1((y[i1,i2]+y[i1,i2+1])/2.00,(v[i1,i2]+v[i1,i2+1])/2.00)
            aa = e21/rr
            bb = e22/rr
            cc = aa+bb+1
            a[i2] = bb/(cc-aa*a[i2-1])

            r21 = (1.0/rr)*(y[i1,i2]+y[i1,i2-1])/2.
            r22 = (1.0/rr)*(y[i1,i2]+y[i1,i2+1])/2.

            L2v = p1*(r22*(v[i1,i2+1]-v[i1,i2]) - r21*(v[i1,i2]-v[i1,i2-1]))
            L2R = q1*(r22*(R[i1,i2+1]-R[i1,i2]) - r21*(R[i1,i2]-R[i1,i2-1]))

            d = y[i1,i2] + f1(y[i1,i2],v[i1,i2])*t2 - L2v - L2R
            b[i2]=(aa*b[i2-1]+d)/(cc-aa*a[i2-1])
        y[i1,N1]=0.
        #e22n = e1((y[i1,N1]+y[i1,N1-1])/2.00,(v[i1,N1]+v[i1,N1-1])/2.00)
        #FDn = p1*(y[i1,N1]+y[i1,N1-1])*(v[i1,N1]-v[i1,N1-1])/2.0 + \
        #      q1*(y[i1,N1]+y[i1,N1-1])*(R[i1,N1]-R[i1,N1-1])/2.0 + \
        #      f1(y[i1,N1],v[i1,N1])*hh/4.0
        #y[i1,N1]=(e22n*b[N1-1] + rr*y[i1,N1]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)

        for ii in range(1,N1+1):
            i2=N1-ii
            y[i1,i2] = a[i2]*y[i1,i2+1]+b[i2]

        for i2 in range(0,N1+1):
            if y[i1,i2]<0: y[i1,i2]=0.0

# v[k+1] to x2
        #e210 = e2((y[i1,1]+y[i1,0])/2.00,(v[i1,1]+v[i1,0])/2.00)
        #a[0] = e210/(e210+rr/2.0)
        #FD0 = (-p2/rr)*(v[i1,1]+v[i1,0])*(y[i1,1]-y[i1,0])- \
        #      (q2/rr)*(v[i1,1]+v[i1,0])*(R[i1,1]-R[i1,0])+f2(y[i1,0],v[i1,0])*t2
        #b[0] = (v[i1,0]+FD0)/(1+2.0*e210/rr)

        for i2 in range(1,N1):
            e21=e2((y[i1,i2]+y[i1,i2-1])/2.00,(v[i1,i2]+v[i1,i2-1])/2.00)
            e22=e2((y[i1,i2]+y[i1,i2+1])/2.00,(v[i1,i2]+v[i1,i2+1])/2.00)
            aa = e21/rr
            bb = e22/rr
            cc = aa+bb+1
            a[i2] = bb/(cc-aa*a[i2-1])

            r21 = (1.0/rr)*(v[i1,i2]+v[i1,i2-1])/2.
            r22 = (1.0/rr)*(v[i1,i2]+v[i1,i2+1])/2.

            L2y = p2*(r22*(y[i1,i2+1]-y[i1,i2])-r21*(y[i1,i2]-y[i1,i2-1]))
            L2R = q2*(r22*(R[i1,i2+1]-R[i1,i2])-r21*(R[i1,i2]-R[i1,i2-1]))

            d = v[i1,i2] + f2(y[i1,i2],v[i1,i2])*t2 - L2y - L2R
            b[i2]=(aa*b[i2-1]+d)/(cc-aa*a[i2-1])
        v[i1,N1]=0.
        #e22n = e2((y[i1,N1]+y[i1,N1-1])/2.00,(v[i1,N1]+v[i1,N1-1])/2.00)
        #FDn = p2*(v[i1,N1]+v[i1,N1-1])*(y[i1,N1]-y[i1,N1-1])/2.0 + \
        #      q2*(v[i1,N1]+v[i1,N1-1])*(R[i1,N1]-R[i1,N1-1])/2.0 + \
        #      f2(y[i1,N1],v[i1,N1])*hh/4.0
        #v[i1,N1]=(e22n*b[N1-1] + rr*v[i1,N1]/2.0 + FDn)/(e22n*(1-a[N1-1])+rr/2.0)

        for ii in range(1,N1+1):
            i2=N1-ii
            v[i1,i2] = a[i2]*v[i1,i2+1]+b[i2]

        for i2 in range(0,N1+1):
            if v[i1,i2]<0: v[i1,i2]=0.0
#graphics
    plt.clf()
    x1 = np.linspace(0., l, N1+1)
    x2 = np.linspace(0., l, N1+1)
    #fig = plt.figure()
    ax = Axes3D(fig, auto_add_to_figure=False)
    fig.add_axes(ax)
    X1, X2 = np.meshgrid(x1, x2)
    plt.xlabel('X')
    plt.ylabel('Y')
    ax.plot_surface(X1, X2, y, rstride=1, cstride=1, cmap='Blues')
    ax.plot_surface(X1, X2, v, rstride=1, cstride=1, cmap='BuGn')
    ax.plot_wireframe(X1, X2, y, rstride=1, cstride=1)
    #plt.contourf(XX1, XX2, y, cmap=plt.cm.gray)
    #plt.colorbar()
    #plt.show()
    #plt.draw()
    plt.pause(0.01)

plt.ioff()
fig = plt.figure()
#ax = Axes3D(fig)
x1 = np.linspace(0., l, N1+1)
x2 = np.linspace(0., l, N1+1)
X1,X2 = np.meshgrid(x1, x2)
plt.xlabel('X')
plt.ylabel('Y')
#ax.plot_surface(X1, X2, y, rstride=1, cstride=1, cmap='Blues')
#ax.plot_surface(X1, X2, v, rstride=1, cstride=1, cmap='BuGn')

          #Liniya urovnya ('ax' turn off before activation)

g = np.linspace(0., 2.0, 21)
plt.clabel(plt.contour(x1,x2,y,g,colors='gray'), inline=True, fontsize=8)
plt.contourf(x1,x2,y,g,cmap=plt.cm.Blues)
#plt.clabel(plt.contour(x1,x2,v,g,colors='gray'), inline=True, fontsize=8)
#plt.contourf(x1,x2,v,g,cmap=plt.cm.BuGn)
plt.colorbar()

plt.show()
#plt.show()
f1 = open('rezy.txt','w')
for i in range(0,N1+1):
    for j in range(0,N1+1):
        f1.write("%f\n" % y[i,j])
f1.close()
f2 = open('rezv.txt','w')
for i in range(0,N1+1):
    for j in range(0,N1+1):
        f2.write("%f\n" % v[i,j])
    f2.write("\n") 
f2.close()

dt = time.perf_counter() - tst
print(f"N2 = {N2}, N1 = {N1}, time solution = {dt:1.3e}, h = {h:1.3e}")

-

→ Ссылка