fork download
copy 
  1. import matplotlib.pyplot as plt
  2. import numpy as np
  3. import scipy.integrate as si
  4. import scipy.interpolate as interp
  5. from array import *
  6. import math
  7.  
  8. #f(x,y)
  9. def f(x,y):
  10.     return -y - x**2
  11.  
  12. def Adams(x, y, num, h):
  13.     for i in range(num-1):
  14.         y.insert(len(y), Euler(y,x,h))
  15.     F = [f(x[0],y[0]), f(x[1],y[1]), f(x[2],y[2]), f(x[3],y[3])]
  16.     coef = 24
  17.     for i  in range(len(x)):
  18.         if i >=  (len(y)):
  19.             y.insert(i, y[i-1] + h/coef*(55*F[3] - 59*F[2] + 37*F[1] - 9*F[0]))
  20.             for j in range(3):
  21.                 F[j] = F[j+1]
  22.             F[3] = f(x[i],y[i])
  23.     return y
  24.  
  25.  
  26. #Метод эйлера
  27. def Euler(y, x, h):
  28.   return y[len(y)-1] + h*f(x[len(y) - 1], y[len(y) - 1]) 
  29.  
  30.  
  31. def RRKM(b, h, e):
  32.      y = array('f', [10])
  33.      x = array('f', [0])
  34.      k = array('f', [0,0,0,0,0,0])
  35.      flag = 0
  36.  
  37.      while x[len(x) - 1]<=b:
  38.         while True:
  39.             k[1] = h*f(x[len(y) - 1], y[len(y) - 1])
  40.             k[2] = h*f(x[len(y) - 1] + 1/3*h, y[len(y) - 1] +1/3*k[1])
  41.             k[3] = h*f(x[len(y) - 1] + 1/3*h, y[len(y) - 1] + 1/6*k[1] + 1/6*k[2])
  42.             k[4] = h*f(x[len(y) - 1] + 1/2*h, y[len(y) - 1] + 1/8*k[1] + 3/8*k[3])
  43.             k[5] = h*f(x[len(y) - 1] + h, y[len(y) - 1] + 1/2*k[1] - 3/2*k[3] + 2*k[4])
  44.             q=(8*k[4] - 9*k[3] - k[5] + 2*k[1])/30
  45.             if abs(q)<e*abs(y[len(y) - 1]):
  46.                 break
  47.             h = h/2
  48.         x.insert(len(x), x[len(x) - 1] + h)
  49.         y.insert(len(y), y[len(y) - 1] + 1/6*k[1] + 2/3*k[4] + 1/6*k[5])
  50.         if abs(q)<=e*abs(y[len(y) - 1])/32:
  51.             h = h*2
  52.      return x,y
  53.  
  54. def RK4(x, h, b, n):
  55.     y = array('f', [10])
  56.     k = array('f', [0,0,0,0,0])
  57.     while len(y)!=len(x):
  58.         k[1] = f(x[len(y) - 1], y[len(y) - 1])
  59.         k[2] = f(x[len(y) - 1] + 1/2*h, y[len(y) - 1] + h*k[1]/2)
  60.         k[3] = f(x[len(y) - 1] + 1/2*h, y[len(y) - 1] + h*k[2]/2)
  61.         k[4] = f(x[len(y) - 1] + h, y[len(y) - 1] + h*k[3])
  62.         y.insert(len(y), y[len(y) - 1] +h*(k[1] + 2*k[2] + 2*k[3] + k[4])/6)
  63.     return y
  64.  
  65.  
  66.  
  67.  
  68. print ("Введите кол-во шагов")
  69. h = int(input())
  70. a = 0
  71.  
  72. print ("Введите конец интервала")
  73. b = int(input())
  74.  
  75. print ("Введите точность")
  76. e = float(input())
  77. y1 = array('f')
  78. x, step = np.linspace(a, b, h, True, True)
  79. y1 = -x*x + 2*x - 2 + 12*math.e**(-x)
  80. y = array('f', [10])
  81. yA = array('f', [10])
  82. yRKM = array('f', [10])
  83.  
  84.  
  85. for i in x:
  86.     if i != x[len(x)-1]:
  87.         y.insert(len(y), Euler(y, x, step))
  88.  
  89. yA = Adams(x,yA,5,step)
  90. xRKM,yRKM = RRKM(b, step, e)
  91. yRK4 = RK4(x, step, b, h)
  92.  
  93.  
  94.  
  95. fig, ax = plt.subplots()
  96. ax.set_xlim(0,5)
  97. ax.set_ylim(-20,20)
  98. ax.set_title("Графики")
  99. ax.set_xlabel("x")
  100. ax.set_ylabel("y")
  101. ax.plot(x, y, color = 'g', label = 'Euler')
  102. ax.plot(x, yRK4, color = 'y', label = 'RK4')
  103. ax.plot(x, yA, color = 'c', label = 'Adams')
  104. ax.plot(xRKM,yRKM, color = 'm', label = 'RKM')
  105. ax.scatter(xRKM,yRKM,c = 'deeppink')
  106. ax.plot(x, y1, color = 'r', label = 'Solution')
  107. ax.legend()
  108. plt.show()
  109.  
Runtime error #stdin #stdout #stderr 0.6s 71500KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Введите кол-во шагов
stderr
copy 
Traceback (most recent call last):
  File "prog.py", line 69, in <module>
EOFError: EOF when reading a line
https://ideone.com/SzTV5s
language:
Python (cpython 2.7.16)
created:
2 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!