API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Nov 14th, 2019
138
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 9.46 KB | None | 0 0
raw download clone embed print report
  1. import math
  2. from matplotlib import pyplot as plot
  3.  
  4.  
  5. class CauchySolver:
  6.     def __init__(self, phi, derivative, x0, xn, k):
  7.         self.phi = phi
  8.         self.derivative = derivative
  9.         self.x0 = x0
  10.         self.xn = xn
  11.         self.k = k
  12.  
  13.     def get_setka(self, n):
  14.         h = (self.xn - self.x0) / n
  15.         x = [self.x0 + i * h for i in range(n + 1)]
  16.         n = n + 1
  17.         return x, n, h
  18.  
  19.     def f(self, x, y):
  20.         return self.derivative(x) + self.k * (y - self.phi(x))
  21.  
  22.     def y_euler(self, n):
  23.         x, n, h = self.get_setka(n)
  24.         y = []
  25.         for i in range(n):
  26.             if i == 0:
  27.                 y.append(self.phi(x[i]))
  28.             else:
  29.                 y.append(y[i - 1] + h * self.f(x[i - 1], y[i - 1]))
  30.         return y
  31.  
  32.  
  33.     def y_adams(self, n):  
  34.         y = self.y_runge_cutte_b(n)[:4]
  35.         x, n, h = self.get_setka(n)
  36.         for i in range(4, n):
  37.             f1 = self.f(x[i - 1], y[i - 1])
  38.             f2 = self.f(x[i - 2], y[i - 2])
  39.             f3 = self.f(x[i - 3], y[i - 3])
  40.             f4 = self.f(x[i - 4], y[i - 4])
  41.             y.append(y[i - 1] + h / 24 * (55 * f1 - 59 * f2 + 37 * f3 - 9 * f4))
  42.         return y
  43.  
  44.     def y_pred_corr(self, n, l):
  45.         x, n, h = self.get_setka(n)
  46.         y = []
  47.         for i in range(n):
  48.             if i == 0:
  49.                 y.append(self.phi(x[i]))
  50.             else:
  51.                 yp = y[i - 1] + h * self.f(x[i - 1], y[i - 1])
  52.                 y.append(yp)
  53.                 # while True:
  54.                     # yc = y[i - 1] + (self.f(x[i - 1], y[i - 1]) + self.f(x[i], yp)) * h / 2
  55.                     # if abs(yc - yp) / abs(yp) <= l:
  56.                         # y.append(yc)
  57.                         # break
  58.                     # yp = yc
  59.  
  60.         return y
  61.  
  62.     def y_abm(self, n):
  63.         y = self.y_runge_cutte_a(n)[:4]
  64.         x, n, h = self.get_setka(n)
  65.         for i in range(4, n):
  66.             f1 = self.f(x[i - 1], y[i - 1])
  67.             f2 = self.f(x[i - 2], y[i - 2])
  68.             f3 = self.f(x[i - 3], y[i - 3])
  69.             f4 = self.f(x[i - 4], y[i - 4])
  70.             yc = y[i - 1] + h / 24 * (55 * f1 - 59 * f2 + 37 * f3 - 9 * f4)
  71.             y.append(y[i - 1] + h  / 24 * (9 * self.f(x[i], yc) + 19 * f1 - 5 * f2 + f3))
  72.         return y
  73.  
  74.     def y_runge_cutte_a(self, n):
  75.         x, n, h = self.get_setka(n)
  76.         y = []
  77.         for i in range(n):
  78.             if i == 0:
  79.                 y.append(self.phi(x[i]))
  80.             else:
  81.                 k1 = self.f(x[i - 1], y[i - 1])
  82.                 k2 = self.f(x[i - 1] + h / 2, y[i - 1] + k1 * h / 2)
  83.                 k3 = self.f(x[i - 1] + h, y[i - 1] - k1 * h + 2 * k2 * h)
  84.                 y.append(y[i - 1] + h * (k1 + 4 * k2 + k3) / 6)
  85.         return y
  86.  
  87.     def y_runge_cutte_b(self, n):
  88.         x, n, h = self.get_setka(n)
  89.         y = []
  90.         for i in range(n):
  91.             if i == 0:
  92.                 y.append(self.phi(x[i]))
  93.             else:
  94.                 k1 = self.f(x[i - 1], y[i - 1])
  95.                 k2 = self.f(x[i - 1] + h / 2, y[i - 1] + k1 * h / 2)
  96.                 k3 = self.f(x[i - 1] + h / 2, y[i - 1] +  k2 * h / 2)
  97.                 k4 = self.f(x[i - 1] + h, y[i - 1] + k3 * h)
  98.                 y.append(y[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6)
  99.         return y
  100.  
  101.     def y_gira(self, n, eps=1e-3):
  102.         y = self.y_runge_cutte_b(n)[:4]
  103.         x, n, h = self.get_setka(n)
  104.         for i in range(4, n):
  105.             yy = 4 * h * self.f(x[i - 1], y[i - 1]) + (y[i - 4] - 10 * y[i - 1]) / 3 - 2 * y[i - 3] + 6 * y[i - 2]
  106.             y.append((48 * y[i - 1] - 36 * y[i - 2] + 16 * y[i - 3] - 3 * y[i - 4] + 12 * h * self.f(x[i], yy)) / 25)
  107.         return y
  108.  
  109.     def local_regarding_eps(self, y):
  110.         x, n, h = self.get_setka(len(y) - 1)
  111.         eps = []
  112.         for i in range(n):
  113.             phi_i = self.phi(x[i])
  114.             eps.append(abs(phi_i - y[i]) / max(abs(phi_i), abs(y[i])))
  115.         return eps
  116.  
  117.     def absolute_eps(self, y):
  118.         x, n, h = self.get_setka(len(y) - 1)
  119.         eps = []
  120.         for i in range(n):
  121.             eps.append(abs(self.phi(x[i]) - y[i]))
  122.         return eps
  123.  
  124.     def valuation_runge(self, yh1, yh2, p):
  125.         diff = []
  126.         for i in range(len(yh1)):
  127.             diff.append((2 ** p / (2 ** p - 1)) * abs(- yh2[2 * i] + yh1[i]))
  128.         return diff
  129.  
  130.     def plot(self, y):
  131.         x, n, h = self.get_setka(len(y) - 1)
  132.         phi_x = [self.x0 + i * (self.xn - self.x0) / 1000 for i in range(1000)]
  133.         phi_y = [self.phi(xi) for xi in phi_x]
  134.         plot.plot(phi_x, phi_y, color='red', label='phi')
  135.         plot.plot(x, y, color='blue', label='y')
  136.         plot.show()
  137.  
  138. def plot_emax(f, error, p, **kwargs):
  139.     a = 100
  140.     b = 1000
  141.     step = 1
  142.  
  143.     label = kwargs['label']
  144.     del kwargs['label']
  145.  
  146.     c = 20
  147.     x = []
  148.     e_max = []
  149.     sf = []
  150.     for n in range(a, b, step):
  151.         x.append(n)
  152.         e_max.append(max(error(f(n=n, **kwargs))))
  153.         # sf.append(c * n ** (-p))
  154.  
  155.     plot.xscale('log')
  156.     plot.yscale('log')
  157.     plot.plot(x, e_max, label = label)
  158.     # plot.plot(x, sf, color='red')
  159.     # plot.show()
  160.  
  161.  
  162. def phi(x):
  163.     return math.cos(x ** 3)
  164.  
  165. def phi_derivative(x):
  166.     return -math.sin(x ** 3) * 3 * x ** 2
  167.  
  168. def main():
  169.     n = 10
  170.     p = 2
  171.    
  172.     solver = CauchySolver(phi=phi, derivative=phi_derivative, x0=-math.pi, xn=math.pi, k=0.8)
  173.  
  174.     # y = solver.y_runge_cutte_b(n=n)
  175.     # y = solver.y_adams(n=n)
  176.     # y2= solver.y_adams(n=2*n)
  177.     # y = solver.y_pred_corr(l = 0.01, n=n)
  178.     # y2 = solver.y_pred_corr(l = 0.01, n=n)
  179.     # y = solver.y_abm(n=n)
  180.     # y2 = solver.y_abm(n=2* n)
  181.     # y = solver.y_euler(n=n)
  182.     # y2 = solver.y_euler(n = 2*n)
  183.     # y = solver.y_runge_cutte_a(n=n)
  184.     # y = solver.y_gira(n=n)
  185.     # y2 = solver.y_gira(n=2*n)
  186.     # y2 = solver.y_runge_cutte_a(n=2 * n)
  187.     # pz = math.log( abs(max(solver.local_regarding_eps(y)) / max(solver.local_regarding_eps(y2))), 2)
  188.     # pz2 = math.log( abs(max(solver.absolute_eps(y)) / max(solver.absolute_eps(y2))), 2)
  189.  
  190.  
  191.     # n = 10
  192.     # e_max = None
  193.  
  194.     # while True:
  195.         # y = solver.y_euler(n=n)
  196.         # e = max(solver.absolute_eps(y))
  197.         # if e_max is not None:
  198.             # if abs(e - e_max) < 10e-3:
  199.                 # print(n, (solver.xn - solver.x0) / n)
  200.                 # break
  201.        
  202.         # n = n * 2
  203.         # e_max = e
  204. # графики е_max
  205.     # aa = [(solver.y_euler, 1, 'эйлер'),(solver.y_adams, 4, 'адамс'),(solver.y_pred_corr, 2, 'итерац'),(solver.y_abm, 4, 'АБМ'),(solver.y_runge_cutte_a, 3, 'РК(3)'),(solver.y_runge_cutte_b, 4, 'РК(4)'), (solver.y_gira, 4, 'Гира)')]
  206.     # for f, p, l in aa:
  207.         # plot_emax(f=f, error=solver.absolute_eps, p=p, label=l)
  208.     # plot.legend()
  209.     # plot.show()
  210. # итерац.метод
  211.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=10000, label='10000')
  212.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-10, label='-10')
  213.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-8, label='-8')
  214.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-6, label='-6')
  215.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-4, label='-4')
  216.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-2, label='-2')
  217.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-1, label='-1')
  218.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-0, label='-0')
  219.     plot.legend()
  220.     plot.show()
  221.  
  222.     # for i in range(20):
  223.         # y = solver.y_abm(n=n)
  224.         # y2 = solver.y_abm(n=2*n)
  225.         # y = solver.y_gira(n=n)
  226.         # y2 = solver.y_gira(n=2*n)
  227.         # y = solver.y_runge_cutte_b(n=n)
  228.         # y2 = solver.y_runge_cutte_b(n=2*n)
  229.         # y = solver.y_runge_cutte_a(n=n)
  230.         # y2 = solver.y_runge_cutte_a(n=2*n)
  231.         # y = solver.y_euler(n=n)
  232.         # y2 = solver.y_euler(n= 2 * n)
  233.         # y = solver.y_adams(n=n)
  234.         # y2 = solver.y_adams(n=2*n)
  235.         # y = solver.y_pred_corr(l = 10 ** (-1), n=n)
  236.         # y2 = solver.y_pred_corr(l = 10 ** (-1), n=2 * n)
  237.         # print(solver.local_regarding_eps(y))
  238.         # pz = math.log( abs(max(solver.local_regarding_eps(y)) / max(solver.local_regarding_eps(y2))), 2)
  239.         # pz2 = math.log( abs(max(solver.absolute_eps(y)) / max(solver.absolute_eps(y2))), 2)
  240.         # print('p* (loc) =  %.2E, p* (abs) = %.2E' % (pz, pz2))
  241.         # print('max локальная относительная %.2E, max абсолютная %.6E' % (max(solver.local_regarding_eps(y)), max(solver.absolute_eps(y))))
  242.         # print('%d оценка по правилу Рунге %.2E, %.2E' % (n, max(solver.valuation_runge(y, y2, p)), max(solver.absolute_eps(y))))
  243.         # plot_emax(f=solver.y_gira, error=solver.absolute_eps, p=p)
  244.         # n = n * 2
  245.     # print('оценка по правилу Рунге %.6E' % (max(solver.valuation_runge(y, y2, p))))
  246.     # print('p* (loc) =  %.6f, p* (abs) = %.6f' % (pz, pz2))
  247.     # print('max локальная относительная %.2E, max абсолютная %.2E' % (max(solver.local_regarding_eps(y)), max(solver.absolute_eps(y))))
  248.  
  249.     # solver.plot(y)
  250.     # solver.plot(solver.absolute_eps)
  251.     # plot_emax(f=solver.y_gira, error=solver.absolute_eps, p=p)
  252.  
  253.     # l10 = [3.14856056E+02, 9.61512887E+02, 3.60834374E+02, 6.84264503E+01, 1.60062067E+01, 3.93924985E+00, 9.81012935E-01, 2.45017209E-01, 6.12395731E-02, 1.53089731E-02, 3.82718587E-03, 9.56793005E-04, 2.39198143E-04, 5.97994648E-05, 1.49499250E-05, 3.73755566E-06, 9.34497514E-07, 2.33762336E-07, 5.83910977E-08, 1.47224043E-08]
  254.     # l6 = [3.14855847E+02, 9.61511494E+02, 3.60834328E+02, 6.84264467E+01, 1.60062119E+01, 3.93924970E+00, 9.81014023E-01, 2.45016802E-01, 6.12410550E-02, 1.53098508E-02, 3.82874286E-03, 9.57928022E-04, 2.38187413E-04, 6.02441899E-05, 1.44482460E-05, 3.61183274E-06, 9.02953146E-07, 2.25740056E-07, 5.64369144E-08, 1.41134897E-08]
  255.     # l4 = [3.14815323E+02, 9.61416445E+02, 3.60816867E+02, 6.84247258E+01, 1.60066843E+01, 3.93913479E+00, 9.81308365E-01, 2.45172731E-01, 6.12491721E-02, 1.54368763E-02, 3.90645617E-03, 9.24641114E-04, 2.31157243E-04, 5.77892692E-05, 1.44472601E-05, 3.61180901E-06, 9.02952284E-07, 2.25740057E-07, 5.64368925E-08, 1.41134897E-08]
  256.     # l2 = [3.12414690E+02, 9.54277666E+02, 3.58308448E+02, 6.84714427E+01, 1.60458342E+01, 3.93507879E+00, 9.65446205E-01, 2.26064204E-01, 5.91838519E-02, 1.47948051E-02, 3.69854349E-03, 9.24627030E-04, 2.31156215E-04, 5.77890156E-05, 1.44472457E-05, 3.61180909E-06, 9.02952284E-07, 2.25740057E-07, 5.64368925E-08, 1.41134897E-08]
  257.  
  258.     # x = [x for x in range(len(l10))]
  259.  
  260.     # plot.plot(x, l10)
  261.     # plot.plot(x, l6)
  262.     # plot.plot(x, l4)
  263.     # plot.plot(x, l2)
  264.     # plot.show()
  265.  
  266. main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!