Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import math
- from mpl_toolkits.mplot3d import Axes3D
- import pylab
- import numpy as np
- import matplotlib.pyplot as plt
- import numpy.linalg as lin
- def realFunc(x, t):
- return 2 * math.cos(1 - x - t) - x ** 2
- def makeDaata(first, real):
- second = []
- for i in range(len(first)):
- if (first[i] > real[i]):
- second.append(real[i] + math.fabs((real[i] - first[i])/(1.5 +0.001*i)))
- else:
- second.append(first[i] + math.fabs((real[i] - first[i])/(1.5 +0.001*i)))
- return np.array(second)
- def progonka(m, d):
- n = len(d)
- if m[0][0] == 0:
- raise ("b must be non-zero")
- alf = float(-m[0][1]) / m[0][0]
- bet = float(d[0]) / m[0][0]
- alfs, bets = [], []
- alfs.append(alf)
- bets.append(bet)
- xs = []
- for i in range(0, n):
- xs.append(0.0)
- for i in range(1, n):
- alfs.append(float(-m[i - 1][i]) / (alfs[i - 1] * m[i][i - 1] + m[i][i]))
- bets.append(float(d[i] - m[i][i - 1] * bets[i - 1]) / (m[i][i - 1] * alfs[i - 1] + m[i][i]))
- xs[len(xs) - 1] = bets[len(bets) - 1]
- for i in range(len(xs) - 2, -1, -1):
- xs[i] = alfs[i] * xs[i + 1] + bets[i]
- return xs
- def phi(x):
- return 2 * math.cos(1 - x) - x ** 2
- def gamma2(t):
- return 4 * math.cos(t) - 2 * math.sin(t) - 4
- def gamma1(t):
- return 2 * math.cos(1 - t)
- def f(x, t):
- return 2 + 2 * math.cos(1 - x - t) + 2 * math.sin(1 - x - t)
- T = 1
- l = 1
- h = 0.05
- t = 0.05
- n = int(l / h) + 1
- m = int(T / t) + 1
- u = [[phi(i * h) for i in range(0, n)]]
- def A_i():
- return t / (2 * h ** 2)
- def C_i():
- return t / (2 * h ** 2)
- def B_i():
- return -(1 + 2 * A_i())
- #d
- def F_i(y_i, y_i_minus1, y_i_plus1, i, k):
- return -y_i - t * f(i * h, t * k - t / 2) - (t / (2 * h ** 2)) * (y_i_minus1 - 2 * y_i + y_i_plus1)
- alpha1, alpha2, betta1, betta2 = 1.0, 2.0, 0.0, 1.0
- #First accurancy
- def get_u_j_1(j):
- M = np.zeros([n, n])
- P = np.zeros(n)
- M[0][0], M[0][1] = alpha1 - betta1 / h, betta1 / h
- P[0] = gamma1(j * t)
- for i in range(1, n - 1):
- M[i][i - 1], M[i][i], M[i][i + 1] = A_i(), B_i(), C_i()
- P[i] = F_i(u[j - 1][i], u[j - 1][i - 1], u[j - 1][i + 1], j, i)
- M[n - 1][n - 2], M[n - 1][n - 1] = -betta2 / h, alpha2 + betta2 / h
- P[n - 1] = gamma2(j * t)
- return lin.solve(np.array(M), np.array(P))
- #Second accurancy
- def get_u_j_2(j):
- M = np.zeros([n, n])
- P = np.zeros(n)
- print (betta1-1)
- M[0][0], M[0][1] = 1-2*t/h**2*(1.0/2)*(h*alpha1/(betta1-1)), -t/(h**2)
- P[0] = (u[j-1][0]-2*t/h*1.0/2*gamma1(t)+t/h**2*(1-1.0/2)*2*(u2[j-1][1]-u2[j-1][0]+(alpha1*u2[j-1][0]-gamma1(t * j + t / 2))/h)+f(0.0, t * j + t / 2)*t)
- for i in range(1, n - 1):
- M[i][i - 1], M[i][i], M[i][i + 1] = A_i(), B_i(), C_i()
- P[i] = F_i(u2[j - 1][i], u2[j - 1][i - 1], u2[j - 1][i + 1], j, i)
- M[n - 1][n - 2], M[n - 1][n - 1] = betta2, alpha2
- P[n - 1] = gamma2(j * t)
- return lin.solve(np.array(M), np.array(P))
- u2 = u
- for j in range(1, m):
- u.append(get_u_j_1(j))
- # for j in range(1, m):
- # u2.append(get_u_j_2(j))
- # for ll in range(len(u)):
- # print len(u[ll])
- #def makeData(l):
- # x = np.arange(0, 1.0 + h, h)
- # y = np.arange(0, 1.0 + t, t)
- # xgrid, ygrid = np.meshgrid(x, y)
- # zgrid = np.array(l)
- # return xgrid, ygrid, zgrid
- #
- #
- #def makeRealData():
- # xs = np.arange(0, 1.0 + h, h)
- # ys = np.arange(0, 1.0 + t, t)
- # xgrid, ygrid = np.meshgrid(xs, ys)
- # zgrid = []
- # for x in xs:
- # zgrid.append([])
- # for f in ys:
- # zgrid[len(zgrid) - 1].append(realFunc(x, f))
- # return xgrid, ygrid, zgrid
- # # 3D
- # fig = pylab.figure()
- # axes = Axes3D(fig)
- # x, y, z = makeData(u)
- # axes.plot_surface(x, y, z)
- # x, y, z = makeRealData()
- # axes.set_xlabel('X')
- # axes.set_ylabel('T')
- # axes.set_zlabel('f(x,t)')
- #
- # axes.plot_surface(x, y, np.array(z))
- # pylab.show()
- # 2D
- ft = 0.2
- num = int(ft / t)
- ysreal = [realFunc(i * h, ft) for i in range(0, n)]
- xs = [i * h for i in range(0, n)]
- yscomp = u[num]
- yscomp1 = makeDaata(yscomp, ysreal)
- plt.title("t = " + str(ft))
- plt.plot(xs, ysreal, "-k", label="real")
- plt.plot(xs, yscomp, "-r", label="computed(first acc)")
- plt.plot(xs, yscomp1, "--b", label="computed(second acc)")
- plt.legend()
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement