matbiz01

Untitled

Apr 29th, 2021
691
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import matplotlib.pyplot as plt
  2.  
  3. from itertools import zip_longest
  4. from math import sin, cos
  5. from math import pi
  6.  
  7. #kod dla metody Hermite'a, parametry działają analogicznie.
  8.  
  9. x0 = 1/3
  10. xn = 3.0
  11. n = 30
  12. czebyszew = 0
  13.  
  14. def getOrigValues(x):
  15.     return x * sin(3 * pi / x)
  16.  
  17. def getDerivValues(x):
  18.     return sin( 3 * pi / x) - 3 * pi * cos(3 * pi / x) / x
  19.  
  20. def calculateBasicPoints(n, czebyszew):
  21.     points = []
  22.     interval = (xn - x0) / (n - 1)
  23.     for x in range(n):
  24.         if czebyszew:
  25.             currX = 1/2 * (xn + x0) + 1/2 * (xn - x0) * cos((2 * (x + 1) - 1) / 2 / n * pi)
  26.         else:
  27.             currX = x0 + x * interval
  28.         points.append([currX, getOrigValues(currX)])
  29.     return points
  30.  
  31. def maxDev(l1, l2):
  32.     temp = 0
  33.     for x,y in zip(l1, l2):
  34.         if abs(x - y) > temp: temp = abs(x - y)
  35.     return temp
  36.  
  37. def makePyramid(xVals, yVals, derVals):
  38.     n = len(xVals)
  39.  
  40.     pyramid = []
  41.     for x,y in zip(xVals, yVals):
  42.         pyramid.append([x, y])
  43.     for x in range(int(n/2)):
  44.         pyramid[x * 2 + 1].append(derVals[x])
  45.     for x in range(int(n/2) - 1):
  46.         pyramid[x * 2 + 2].append( (pyramid[x * 2 + 2][1] - pyramid[x * 2 + 1][1]) / (xVals[x * 2 + 2] - xVals[x * 2 + 1]) )
  47.  
  48.     for x in range(2, n):
  49.         for y in range(x - 1):
  50.             pyramid[x].append( (pyramid[x][2 + y] - pyramid[x - 1][2 + y]) / ( xVals[x] - xVals[x - (2 + y)] ))
  51.  
  52.     coeffs = []
  53.     for x in range(n):
  54.         coeffs.append(pyramid[x][x + 1])
  55.     return coeffs
  56.  
  57. def getHermiteVal(coefficients, x, xVals):
  58.     points = list(reversed(xVals))
  59.     coeffs = list(reversed(coefficients))
  60.     total = coeffs[0]
  61.     for i in range(1, len(coeffs)):
  62.         total *= (x - points[i])
  63.         total += coeffs[i]
  64.     return total
  65.  
  66.  
  67. def do(n):
  68.     points = calculateBasicPoints(n, czebyszew)
  69.     xVal = [x for x, y in points]
  70.     yVal = [y for x, y in points]
  71.     xVals = []
  72.     yVals = []
  73.     derVals = [getDerivValues(x) for x in xVal]
  74.     for x in xVal:
  75.         xVals.append(x)
  76.         xVals.append(x)
  77.     for y in yVal:
  78.         yVals.append(y)
  79.         yVals.append(y)
  80.  
  81.     coefficients = makePyramid(xVals, yVals, derVals)
  82.  
  83.     point_count = 500
  84.     pointX = []
  85.     pointY = []
  86.     pointYorig = []
  87.     interval = (xn - x0) / point_count
  88.     for x in range(point_count):
  89.         currX = x0 + x * interval
  90.         pointX.append(currX)
  91.         pointY.append(getHermiteVal(coefficients, currX, xVals))
  92.         pointYorig.append(getOrigValues(currX))
  93.  
  94.     print(n, maxDev(pointY, pointYorig))
  95.  
  96.  
  97.     plt.plot(pointX, pointY, label="wielomian interpolujący")
  98.  
  99.  
  100.     plt.plot(pointX, pointYorig, label="oryginalna funkcja")
  101.     for point in calculateBasicPoints(n, czebyszew):
  102.         plt.scatter(point[0], point[1], s=50)
  103.     plt.legend(loc=1)
  104.     plt.show()
  105.  
  106. do(n)
  107.  
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×