Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import math
- import matplotlib.pyplot as plt
- import numpy as np
- import numpy
- from mpl_toolkits.mplot3d import Axes3D
- import pylab
- def func1(j,h):
- x = j*h
- return x**2
- def func2(j,h,t):
- x = j*h
- return x*t+x**2
- def gamma1(k,t):
- t = t*k
- return -t
- def gamma2(k,t):
- t = k*t
- return 2+t/((math.cosh(t)))**2
- def realfunc(j,k,h,t):
- x = j*h
- t = k*t
- return x**2+(math.tanh(t*x))
- def square(u,r):
- x = 0.0
- for k in range (len(u)):
- for j in range (len(u[k])):
- x+=(u[k][j]-r[k][j])**2
- return x
- def f(j,k,h,t):
- x = h*j
- t1 = t*j
- f = -1 +((t1)**2-2*(x**2))*(math.tanh(x*t1)/(math.cosh(x*t1))**2)
- return f
- l = 1
- T = 1
- h = 0.05
- t = 0.05
- n = int(l/h)
- m = int(T/t)
- g = (t**2)/(h**2)/2
- u =[
- [],[]
- ]
- for j in range(0,n):
- u[0].append(func1(j,h))
- u[1].append(func2(j,h,t))
- for k in range(2,m):
- u.append([])
- for j in range(1,n-1):
- u[k].append((g**2)*u[k-1][j-1]+(2 - 2*(g**2))*u[k-1][j]+(g**2)*u[k-1][j+1]-u[k-2][j])
- u[k].insert(0,(1/(1+h)*u[k][0]+h*gamma1(k,t)/(h+1)))
- u[k].append(u[k][n-2]+h*gamma2(k,t))
- u2 = [
- [],[]
- ]
- for j in range(0,m):
- u2[0].append(func1(j,h))
- for j in range(1,m-1):
- u2[1].append(u2[0][j]+t*func2(j,h,t)/2+f(j,0,h,t)*(t**2)/2+(u2[0][j+1]-2*u2[0][j]+u2[0][j-1])*(t**2)/(4*(h**2)))
- u2[1].insert(0,((4*u2[1][0]-u2[1][1]-gamma1(0,t)*2*h)/3))
- u2[1].append(2*h*gamma2(m,t)/3+4*u[1][m-2]/3-(u[1][m-3])/3)
- for k in range(2,m):
- u2.append([])
- for j in range(1,n-1):
- u2[k].append((g**2)*u2[k-1][j-1]+(2 - 2*(g**2))*u2[k-1][j]+(g**2)*u2[k-1][j+1]-u2[k-2][j])
- u2[k].insert(0,(4*u2[k][0]-u2[k][1]+2*gamma1(k,t)*h)/(3+2*h))
- u2[k].append((2*h*gamma2(k,t)-u2[k][n-3]+4*u2[k][n-2])/3)
- u0 = []
- for k in range(0,n):
- u0.append([])
- for j in range(0,m):
- u0[k].append(realfunc(j,k,h,t))
- print(np.matrix(u0))
- print('jopa')
- print(np.matrix(u))
- print('pisya')
- print(np.matrix(u2))
- print('asd',square(u,u0))
- print('asd',square(u2,u0))
- def makeData ():
- x = numpy.arange(0,1,0.05)
- y = numpy.arange(0,1,0.05)
- xgrid, ygrid = numpy.meshgrid(x, y)
- zgrid = numpy.array(u)
- return xgrid, ygrid, zgrid
- def makeDataReal():
- x = numpy.arange(0,1,0.05)
- y = numpy.arange(0,1,0.05)
- xgrid, ygrid = numpy.meshgrid(x, y)
- zgrid = numpy.array(u0)
- return xgrid, ygrid, zgrid
- def makeDataRealSection(t):
- x = numpy.arange(0,1,0.05)
- zgrid = numpy.array(u0[t])
- return x, zgrid
- def makeData2Section(t):
- x = numpy.arange(0,1,0.05)
- zgrid = numpy.array(u2[t])
- return x, zgrid
- def makeDataSection(t):
- x = numpy.arange(0,1,0.05)
- zgrid = numpy.array(u[t])
- return x, zgrid
- def makeData2():
- x = numpy.arange(0,1,0.05)
- y = numpy.arange(0,1,0.05)
- xgrid, ygrid = numpy.meshgrid(x, y)
- zgrid = numpy.array(u2)
- return xgrid, ygrid, zgrid
- #1
- fig = pylab.figure()
- axes = Axes3D(fig)
- x, y, z = makeData()
- axes.plot_surface(x,y,z)
- axes.set_title("№1 orange - real , blue - numerical")
- x, y, z = makeDataReal()
- axes.set_xlabel('X')
- axes.set_ylabel('T')
- axes.set_zlabel('Z')
- axes.plot_surface(x,y,z)
- #2
- fig = pylab.figure()
- axes = Axes3D(fig)
- x, y, z = makeData2()
- axes.plot_surface(x,y,z)
- axes.set_title("№2 orange - real , blue - numerical")
- x, y, z = makeDataReal()
- axes.set_xlabel('X')
- axes.set_ylabel('T')
- axes.set_zlabel('Z')
- axes.plot_surface(x,y,z)
- pylab.show()
- #3 Section t = 0.5
- def generate2dData(t):
- x1, y1 = makeDataSection(t)
- x2, y2 = makeData2Section(t)
- x3, y3 = makeDataRealSection(t)
- return x1,y1,x2,y2,x3,y3
- def twoDemensionalGraphic(x1,y1,x2,y2,x3,y3,t):
- plt.xlabel(r'$x$')
- plt.ylabel(r'$x$')
- plt.title("t =" + str(t))
- plt.plot(x1, y1,"-k", label = "Первая точность")
- plt.plot(x2, y2, "-r", label ="Вторая точность")
- plt.plot(x3,y3, "-b", label="Реальное значение")
- plt.legend()
- plt.grid(True)
- plt.show()
- # меняй s(типо t), при котором будет строиться 2d график
- s = 0.1
- k = 0
- h = 0.0
- while True:
- h += t
- if (h >= s):
- break
- k += 1
- print(k,m)
- if (k>= m):
- print("too big t")
- else:
- x1,y1,x2,y2,x3,y3 = generate2dData(k)
- twoDemensionalGraphic(x1,y1,x2,y2,x3,y3,s)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement