• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Nov 14th, 2019 117 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import math as m
2.
3.
4. def myfunc(x):
5.     return m.cos(x) * x
6.
7.
8. def f(x, g):
9.     return m.cos(x) - m.sin(x) * x + k * (g - myfunc(x))
10.
11.
12. def absErr(arr):
13.     abserr = [abs(myfunc(x[i]) - arr[i]) for i in range(0, N + 1)]
14.     return abserr
15.
16.
17. def globerr(arr):
18.     return max(absErr(arr))
19.
20.
21. def k1(i, t):
22.     return f(x[i - 1], t[i - 1])
23.
24.
25. def k2(i, t):
26.     return f(x[i - 1] + h / 2, t[i - 1] + k1(i, t) * h / 2)
27.
28.
29. def k3a(i, t):
30.     return f(x[i - 1] + h, t[i - 1] - h * k1(i, t) + 2 * h * k2(i, t))
31.
32.
33. def k3b(i, t):
34.     return f(x[i - 1] + h / 2, t[i - 1] + h * k2(i, t) / 2)
35.
36.
37. def k4(i, t):
38.     return f(x[i - 1] + h, t[i - 1] + h * k3b(i, t))
39.
40.
41. def euler():
42.     euler = start.copy()
43.     for i in range(1, N + 1):
44.         euler[i] = euler[i - 1] + h * f(x[i - 1], euler[i - 1])
45.     return euler
46.
47.
48. def rka():
49.     rka = start.copy()
50.     for i in range(1, N + 1):
51.         rka[i] = rka[i - 1] + h * (k1(i, rka) + 4 * k2(i, rka) + k3a(i, rka)) / 6
52.     return rka
53.
54.
55. def rkb():
56.     rkb = start.copy()
57.     for i in range(1, N + 1):
58.         rkb[i] = rkb[i - 1] + h * (k1(i, rkb) + 2 * k2(i, rkb) + 2 * k3b(i, rkb) + k4(i, rkb)) / 6
59.     return rkb
60.
61.
64.     for i in range(4, N + 1):
65.         adams_a[i] = adams_a[i - 1] + (h / 24) * (55 * f(x[i - 1], adams_a[i - 1]) - 59 * f(x[i - 2], adams_a[i - 2]) +
66.                                                   37 * f(x[i - 3], adams_a[i - 3]) - 9 * f(x[i - 4], adams_a[i - 4]))
68.
69.
72.     for i in range(4, N + 1):
73.         adams_b[i] = adams_b[i - 1] + (h / 24) * (55 * f(x[i - 1], adams_b[i - 1]) - 59 * f(x[i - 2], adams_b[i - 2]) +
74.                                                   37 * f(x[i - 3], adams_b[i - 3]) - 9 * f(x[i - 4], adams_b[i - 4]))
76.
77.
78. def abm_a():
80.     abm_a = startrka.copy()
81.     for i in range(4, N + 1):
82.         abm_a[i] = abm_a[i - 1] + h * (9 * f(x[i], predicta[i]) + 19 * f(x[i - 1], abm_a[i - 1]) -
83.                                        5 * f(x[i - 2], abm_a[i - 2]) + f(x[i - 3], abm_a[i - 3]))
84.     return abm_a
85.
86.
87. def abm_b():
89.     abm_b = startrkb.copy()
90.     for i in range(4, N + 1):
91.         abm_b[i] = abm_b[i - 1] + h * (9 * f(x[i], predictb[i]) + 19 * f(x[i - 1], abm_b[i - 1]) -
92.                                        5 * f(x[i - 2], abm_b[i - 2]) + f(x[i - 3], abm_b[i - 3]))
93.     return abm_b
94.
95.
96. def gir_a():
97.     gir_a = startrka.copy()
98.     for i in range(4, N + 1):
99.         gir_a[i] = (48 * gir_a[i - 1] - 36 * gir_a[i - 2] + 16 * gir_a[i - 3] -
100.                     3 * gir_a[i - 4] + 12 * h * (m.cos(x[i]) - m.sin(x[i]) * x[i] -
101.                     k * myfunc(x[i]))) / (25 - k * 12 * h)
102.     return gir_a
103.
104.
105. def gir_b():
106.     gir_b = startrkb.copy()
107.     for i in range(4, N + 1):
108.         gir_b[i] = (48 * gir_b[i - 1] - 36 * gir_b[i - 2] + 16 * gir_b[i - 3] -
109.                     3 * gir_b[i - 4] + 12 * h * (m.cos(x[i]) - m.sin(x[i]) * x[i] -
110.                                                  k * myfunc(x[i]))) / (25 - k * 12 * h)
111.     return gir_b
112.
113.
114. def pred_corr():
115.     pred = start.copy()
116.     corr = [0 for i in range(0, N + 1)]
117.     corr[0] = pred[0]
118.     stop = 1e-10
119.     for i in range(1, N + 1):
120.         pred[i] = pred[i - 1] + h * f(x[i - 1], pred[i - 1])
121.         corr[i] = corr[i - 1] + (h / 2) * (f(x[i - 1], corr[i - 1]) + f(x[i], pred[i]))
122.         z = 0
123.         while abs(corr[i] - pred[i]) / abs(pred[i]) >= stop:
124.             pred[i] = corr[i]
125.             corr[i] = corr[i - 1] + (h / 2) * (f(x[i - 1], corr[i - 1]) + f(x[i], corr[i]))
126.             z += 1
127.     return corr
128.
129.
130. k = 0
131. N = 1000000
132. h = 2 / N
133. x = [-1 + i * h for i in range(N + 1)]
134. y = [0 for i in range(N + 1)]
135. y[0] = myfunc(x[0])
136. start = [0 for i in range(0, N + 1)]
137. start[0] = myfunc(x[0])
138. startrka = [0 for i in range(0, N + 1)]
139. startrkb = [0 for i in range(0, N + 1)]
140. for i in range(0, 4):
141.     startrka[i] = rka()[i]
142.     startrkb[i] = rkb()[i]
143. # plot = open('J:\Test2\ ' + str(N) + '.txt', 'w')
144. # for i in range(N+1):
145. #    plot.write(str(x[i]) + ' ' + str(myfunc(x[i])) + ' ' + str(adams_a()[i]) + '\n')
146. # plot.close()
147. test = open('H:\Test\ ' + str(N) + '.txt', 'w')
148. test.write(str(globerr(euler())) + ' ' + str(globerr(rka())) + ' ' + str(globerr(rkb())) + ' ' +
149.            str(globerr(adams_a())) + ' ' + str(globerr(adams_b())) +  ' ' + str(globerr(gir_a())) + ' ' +
150.            str(globerr(gir_b())) + ' ' + str(globerr(pred_corr())) + '\n')
151. test.close()
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.

Top