4Физика3 задание

import pylab as plb
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

a = np.pi
b = 1
T = 3

Nt = 50
Nx = 100
Ny = 100

hx = a / Nx
ht = T / Nt
hy = b / Ny


u = np.empty((Nt+1, Nx+1, Ny+1))

def U0(x, y):
    return 0

def mu1(y, t):
    return 1

def mu2(y, t):
    return 1

def v1(x, t):
    return np.cos(2*x)

def v2(x, t):
    return np.cos(2*x)

def f(x, y, t):
    return 0.5* np.power(t, 2/5)*np.exp(-t)



for i in range(Nx+1): # H. y.
    for k in range(Ny+1):
        u[0][i][k] = U0(i * hx, k * hy)


for j in range(Nt+1):# Г. у.
    for k in range(Ny+1):
        u[j][0][k] = mu1(k * hy, j * ht)
        u[j][Nx][k] = mu2(k * hy, j * ht)

for j in range(Nt+1):# Г. у.
    for i in range(Nx+1):
        u[j][i][0] = v1(i * hx, j * ht)
        u[j][i][Ny] = v2(i * hx, j * ht)

for j in range(1, Nt+1):
    for i in range(1, Nx):
        for k in range(1 ,Ny):
            u[j][i][k] = ht *( (u[j-1][i-1][k] - 2*u[j-1][i][k] + u[j-1][i+1][k])/(hx *hx) + (u[j-1][i][k-1] - 2*u[j-1][i][k] + u[j-1][i][k+1])/(hy *hy) + f(i * hx, k * hy, (j-1)*ht)) + u[j-1][i][k]

x = np.linspace(0, a, Nx+1)
y = np.linspace(0, b, Ny+1)

xgrid, ygrid = np.meshgrid(x, y)# для изолинии оси

zgrid=u[14]# допустим, возьмем 15 по счету время
#изолинии
cs = plb.contour(xgrid, ygrid, zgrid)
plb.contour(cs)
plb.show()