SHARE
TWEET

Untitled

a guest Nov 14th, 2019 124 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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.                 if l is None:
  53.                     y.append(yp)
  54.                 else:
  55.                     while True:
  56.                         yc = y[i - 1] + (self.f(x[i - 1], y[i - 1]) + self.f(x[i], yp)) * h / 2
  57.                         if abs(yc - yp) / abs(yp) <= l:
  58.                             y.append(yc)
  59.                             break
  60.                         yp = yc
  61.  
  62.         return y
  63.  
  64.     def y_abm(self, n):
  65.         y = self.y_runge_cutte_a(n)[:4]
  66.         x, n, h = self.get_setka(n)
  67.         for i in range(4, n):
  68.             f1 = self.f(x[i - 1], y[i - 1])
  69.             f2 = self.f(x[i - 2], y[i - 2])
  70.             f3 = self.f(x[i - 3], y[i - 3])
  71.             f4 = self.f(x[i - 4], y[i - 4])
  72.             yc = y[i - 1] + h / 24 * (55 * f1 - 59 * f2 + 37 * f3 - 9 * f4)
  73.             y.append(y[i - 1] + h  / 24 * (9 * self.f(x[i], yc) + 19 * f1 - 5 * f2 + f3))
  74.         return y
  75.  
  76.     def y_runge_cutte_a(self, n):
  77.         x, n, h = self.get_setka(n)
  78.         y = []
  79.         for i in range(n):
  80.             if i == 0:
  81.                 y.append(self.phi(x[i]))
  82.             else:
  83.                 k1 = self.f(x[i - 1], y[i - 1])
  84.                 k2 = self.f(x[i - 1] + h / 2, y[i - 1] + k1 * h / 2)
  85.                 k3 = self.f(x[i - 1] + h, y[i - 1] - k1 * h + 2 * k2 * h)
  86.                 y.append(y[i - 1] + h * (k1 + 4 * k2 + k3) / 6)
  87.         return y
  88.  
  89.     def y_runge_cutte_b(self, n):
  90.         x, n, h = self.get_setka(n)
  91.         y = []
  92.         for i in range(n):
  93.             if i == 0:
  94.                 y.append(self.phi(x[i]))
  95.             else:
  96.                 k1 = self.f(x[i - 1], y[i - 1])
  97.                 k2 = self.f(x[i - 1] + h / 2, y[i - 1] + k1 * h / 2)
  98.                 k3 = self.f(x[i - 1] + h / 2, y[i - 1] +  k2 * h / 2)
  99.                 k4 = self.f(x[i - 1] + h, y[i - 1] + k3 * h)
  100.                 y.append(y[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6)
  101.         return y
  102.  
  103.     def y_gira(self, n, eps=1e-3):
  104.         y = self.y_runge_cutte_b(n)[:4]
  105.         x, n, h = self.get_setka(n)
  106.         for i in range(4, n):
  107.             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]
  108.             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)
  109.         return y
  110.  
  111.     def local_regarding_eps(self, y):
  112.         x, n, h = self.get_setka(len(y) - 1)
  113.         eps = []
  114.         for i in range(n):
  115.             phi_i = self.phi(x[i])
  116.             eps.append(abs(phi_i - y[i]) / max(abs(phi_i), abs(y[i])))
  117.         return eps
  118.  
  119.     def absolute_eps(self, y):
  120.         x, n, h = self.get_setka(len(y) - 1)
  121.         eps = []
  122.         for i in range(n):
  123.             eps.append(abs(self.phi(x[i]) - y[i]))
  124.         return eps
  125.  
  126.     def valuation_runge(self, yh1, yh2, p):
  127.         diff = []
  128.         for i in range(len(yh1)):
  129.             diff.append((2 ** p / (2 ** p - 1)) * abs(- yh2[2 * i] + yh1[i]))
  130.         return diff
  131.  
  132.     def plot(self, y):
  133.         x, n, h = self.get_setka(len(y) - 1)
  134.         phi_x = [self.x0 + i * (self.xn - self.x0) / 1000 for i in range(1000)]
  135.         phi_y = [self.phi(xi) for xi in phi_x]
  136.         plot.plot(phi_x, phi_y, color='red', label='phi')
  137.         plot.plot(x, y, color='blue', label='y')
  138.         plot.show()
  139.  
  140. def plot_emax(f, error, p, **kwargs):
  141.     a = 100
  142.     b = 1000
  143.     step = 1
  144.  
  145.     label = kwargs['label']
  146.     del kwargs['label']
  147.  
  148.     c = 20
  149.     x = []
  150.     e_max = []
  151.     sf = []
  152.     for n in range(a, b, step):
  153.         x.append(n)
  154.         e_max.append(max(error(f(n=n, **kwargs))))
  155.         # sf.append(c * n ** (-p))
  156.  
  157.     plot.xscale('log')
  158.     plot.yscale('log')
  159.     plot.plot(x, e_max, label = label)
  160.     # plot.plot(x, sf, color='red')
  161.     # plot.show()
  162.  
  163.  
  164. def phi(x):
  165.     return math.cos(x ** 3)
  166.  
  167. def phi_derivative(x):
  168.     return -math.sin(x ** 3) * 3 * x ** 2
  169.  
  170. def main():
  171.     n = 10
  172.     p = 2
  173.    
  174.     solver = CauchySolver(phi=phi, derivative=phi_derivative, x0=-math.pi, xn=math.pi, k=0.8)
  175.  
  176.     # y = solver.y_runge_cutte_b(n=n)
  177.     # y = solver.y_adams(n=n)
  178.     # y2= solver.y_adams(n=2*n)
  179.     # y = solver.y_pred_corr(l = 0.01, n=n)
  180.     # y2 = solver.y_pred_corr(l = 0.01, n=n)
  181.     # y = solver.y_abm(n=n)
  182.     # y2 = solver.y_abm(n=2* n)
  183.     # y = solver.y_euler(n=n)
  184.     # y2 = solver.y_euler(n = 2*n)
  185.     # y = solver.y_runge_cutte_a(n=n)
  186.     # y = solver.y_gira(n=n)
  187.     # y2 = solver.y_gira(n=2*n)
  188.     # y2 = solver.y_runge_cutte_a(n=2 * n)
  189.     # pz = math.log( abs(max(solver.local_regarding_eps(y)) / max(solver.local_regarding_eps(y2))), 2)
  190.     # pz2 = math.log( abs(max(solver.absolute_eps(y)) / max(solver.absolute_eps(y2))), 2)
  191.  
  192.  
  193.     # n = 10
  194.     # e_max = None
  195.  
  196.     # while True:
  197.         # y = solver.y_euler(n=n)
  198.         # e = max(solver.absolute_eps(y))
  199.         # if e_max is not None:
  200.             # if abs(e - e_max) < 10e-3:
  201.                 # print(n, (solver.xn - solver.x0) / n)
  202.                 # break
  203.        
  204.         # n = n * 2
  205.         # e_max = e
  206. # графики е_max
  207.     # 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, 'Гира)')]
  208.     # for f, p, l in aa:
  209.         # plot_emax(f=f, error=solver.absolute_eps, p=p, label=l)
  210.     # plot.legend()
  211.     # plot.show()
  212. # итерац.метод
  213.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=None, label='10000')
  214.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-10, label='-10')
  215.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-8, label='-8')
  216.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-6, label='-6')
  217.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-4, label='-4')
  218.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-2, label='-2')
  219.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-1, label='-1')
  220.     plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-0, label='-0')
  221.     plot.legend()
  222.     plot.show()
  223.  
  224.     # for i in range(20):
  225.         # y = solver.y_abm(n=n)
  226.         # y2 = solver.y_abm(n=2*n)
  227.         # y = solver.y_gira(n=n)
  228.         # y2 = solver.y_gira(n=2*n)
  229.         # y = solver.y_runge_cutte_b(n=n)
  230.         # y2 = solver.y_runge_cutte_b(n=2*n)
  231.         # y = solver.y_runge_cutte_a(n=n)
  232.         # y2 = solver.y_runge_cutte_a(n=2*n)
  233.         # y = solver.y_euler(n=n)
  234.         # y2 = solver.y_euler(n= 2 * n)
  235.         # y = solver.y_adams(n=n)
  236.         # y2 = solver.y_adams(n=2*n)
  237.         # y = solver.y_pred_corr(l = 10 ** (-1), n=n)
  238.         # y2 = solver.y_pred_corr(l = 10 ** (-1), n=2 * n)
  239.         # print(solver.local_regarding_eps(y))
  240.         # pz = math.log( abs(max(solver.local_regarding_eps(y)) / max(solver.local_regarding_eps(y2))), 2)
  241.         # pz2 = math.log( abs(max(solver.absolute_eps(y)) / max(solver.absolute_eps(y2))), 2)
  242.         # print('p* (loc) =  %.2E, p* (abs) = %.2E' % (pz, pz2))
  243.         # print('max локальная относительная %.2E, max абсолютная %.6E' % (max(solver.local_regarding_eps(y)), max(solver.absolute_eps(y))))
  244.         # print('%d оценка по правилу Рунге %.2E, %.2E' % (n, max(solver.valuation_runge(y, y2, p)), max(solver.absolute_eps(y))))
  245.         # plot_emax(f=solver.y_gira, error=solver.absolute_eps, p=p)
  246.         # n = n * 2
  247.     # print('оценка по правилу Рунге %.6E' % (max(solver.valuation_runge(y, y2, p))))
  248.     # print('p* (loc) =  %.6f, p* (abs) = %.6f' % (pz, pz2))
  249.     # print('max локальная относительная %.2E, max абсолютная %.2E' % (max(solver.local_regarding_eps(y)), max(solver.absolute_eps(y))))
  250.  
  251.     # solver.plot(y)
  252.     # solver.plot(solver.absolute_eps)
  253.     # plot_emax(f=solver.y_gira, error=solver.absolute_eps, p=p)
  254.  
  255.     # 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]
  256.     # 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]
  257.     # 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]
  258.     # 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]
  259.  
  260.     # x = [x for x in range(len(l10))]
  261.  
  262.     # plot.plot(x, l10)
  263.     # plot.plot(x, l6)
  264.     # plot.plot(x, l4)
  265.     # plot.plot(x, l2)
  266.     # plot.show()
  267.  
  268. main()
RAW Paste Data
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!
 
Top