Oct 24th, 2021
540
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
2.     flag = 0
3.     N = 10
4.     help_err = 0
5.
6.     plt.xlabel('x')
7.     plt.ylabel('y')
8.
9.     while True:
10.
11.         step = (B - A) / (N - 1)
12.         u_help = np.empty(N)
13.         arr_err = np.empty(N)
14.
15.         u_help[0] = 0
16.         u_help[1] = u_help[0] + step * g(0, u_help[0])
17.         u_help[2] = u_help[1] + step * g(step * 1, u_help[1])
18.         arr_err[0] = 0
19.
20.         error = 0
21.
22.         for i in range(3, N):
23.             a_1 = g(step * (i - 3), u_help[i - 3])
24.             a_2 = 5 * g(step * (i - 2), u_help[i - 2])
25.             a_3 = 19 * g(step * (i - 1), u_help[i - 1])
26.
27.             f_non_lin = y - u_help[i - 1] - (step / 24) * (a_3 - a_2 + a_1) - (9 / 24) * step * (1 / (i * step + 2) + p * pow((y - phi_3(i * step)), 5))
28.             f_n_l = lambdify(y, f_non_lin)
29.             der_f_n_l = f_non_lin.diff(y)
30.             deriv = lambdify(y, der_f_n_l)
31.             y_0, e = phi_3(i * step), 0.001
32.
33.             while np.fabs(f_n_l(y_0)) > e:
34.                 y_0 = y_0 - f_n_l(y_0) / deriv(y_0)
35.
36.             u_help[i] = y_0
37.
38.             if error < np.fabs(np.double(u_help[i]) - np.double(phi_3(i * step))):
39.                 error = np.fabs(u_help[i] - phi_3(i * step))
40.
41.
42.         if error < eps:
43.             flag = 1
44.
45.         print(error, end = '\t')
46.         if help_err != 0:
47.             print("p = ", np.log2(help_err / error))
48.         help_err = error
49.
50.         a = np.linspace(A, B, 160)
51.         b = phi_3(a)
52.         plt.grid()
53.         plt.plot(a, b, 'r')
54.
55.         c = np.linspace(A, B, N)
56.         plt.plot(c, arr_err, 'g')
57.         plt.plot(c, u_help, 'b')
58.         plt.pause(5)
59.
60.         if flag == 1:
61.             break
62.         N = 2 * N
63.
64.
65.     plt.show()
66.
67.     return u_help
RAW Paste Data