3VP_2_Laksa

import numpy
import pylab
from matplotlib import mlab
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D



# точное решение
def ft(t, x):
    return numpy.sinh(2*(x-a*t)) + (x - a*t)**3 + 0.5

# начальное условие (x,0)
def fn(x):
    return  numpy.sinh(2*x) + x**3 + 0.5


# левое граничное условие (0,t)
def fl(t):
    return numpy.sinh(-(2*a*t)) + (-a*t)**3 +0.5

a = 0.1

t_knots = 11
x_knots = 11

t = 1/(t_knots-1) #шаг по t
h = 1/(x_knots-1) #шаг по х

#начальные условия
U = [[0]*x_knots for i in range(t_knots)]
for k in range (x_knots):
    U[0][k] = fn(k*h)

#метод

for j in range (1, t_knots):
    U[j][0] = fl(j*t)
    for k in range (1, x_knots-1):
        U[j][k] = (U[j-1][k+1] + U[j-1][k-1])/2 - a*t*(U[j-1][k+1] + U[j-1][k-1])/(2*h)
    U[j][x_knots-1] = U[j-1][x_knots-1] - a*t * (U[j-1][x_knots-1] - U[j-1][x_knots-2])/h


print ('матрица приближенного решения :')
for j in range(t_knots):
    for k in range (x_knots):
        print("%.4f" %U[j][k] ,end = ' ')
    print()

#вычисление погрешности
max = 0
for j in range(t_knots):
    for k in range(x_knots):
        if max < numpy.fabs(U[j][k] - ft(j*t, k*h)):
            max = numpy.fabs(U[j][k] - ft(j*t, k*h))

print()
print("h^2/t + h^2 + t:", h**2/t + h**2 +t, " максимальная погрешность: ", max)


x = 0.5
a = 0
b = 1

#рисование приближенного решения
tlist = mlab.frange(a, b, t)
zlist = [U[t][round(1/h*x)] for t in range (t_knots)]
pylab.plot(tlist, zlist, color = "red",label = 'Приближенное решение')

#точное решение
tlist = mlab.frange(a, b, t) #массив значений t
flist = [ft(t,x) for t in tlist]  #х=0.5, t идет от 0 до 1 по 0.1
for i in range(0,t_knots):
    print(i/t_knots-1 , " : ", "%.4f" % flist[i])
print('Точное решение для х =', x, 'и для значений t от',a ,'до',b,'с шагом', t , ':')
pylab.plot(tlist,flist, label = 'Точное решение')
pylab.legend()
pylab.xlabel('t')
pylab.show()




x = numpy.arange(a, b + h, h)
y = numpy.arange(a, b + t, t)
xgrid, ygrid = numpy.meshgrid(x, y) #двумерная матрица-сетка
zgrid = ft(ygrid, xgrid) #значения в узлах сетки
fig = pylab.figure()
axes = Axes3D(fig)
axes.plot_surface(xgrid, ygrid, zgrid,  cmap=cm.jet)
pylab.xlabel('x')
pylab.ylabel('t')
pylab.show()