• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Mar 26th, 2019 87 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. restart;
2. #
3. a := -1.;
4. #
5. b := 0.;
6. #
7. H := 0.786;
8. t := 1/2;
9. #
10. yy:=(x)->sin(x^3+1);
11. #
12. y0:=yy(a);
13. #
14. F:=(x,y)->3*x^2*cos(x^3+1)-sin(x^3+1)+y;
15. #
16. eps:=0.0001;
17. #
18. RGKSTEP := proc(XOLD, YOLD, H, t)
19. local K1, K2, K3, K4, xnew, ynew;
20. K1 := H * F(XOLD,YOLD);
21. K2 := H * F(XOLD + H/2, YOLD + K1/2);
22. K3 := H * F(XOLD + H, YOLD + 2 * K2 - K1);
23. ynew := evalf(YOLD + (1/6) * (K1 + 4 * K2 + K3));
24. xnew := XOLD + H;
25. return (xnew, ynew);
26. end proc:
27. #
28. RGKSTEP(a, y0, H, t);
29. yy(a+H);
30. #
31. RGKH := proc(X0, Y0, B, n, t)
32. local X, Y, xx, yy, i, H;
33. H:=(B-X0)/n;
34. X := array(0..n):
35. Y := array(0..n):
36. X[0] := X0;
37. Y[0] := Y0;
38. xx := X0;
39. yy := Y0;
40. for i from 1 to n do
41.    X[i], Y[i] := RGKSTEP(xx, yy, H, t);
42.    xx := X[i];
43.    yy := Y[i];
44. end do:
45. return X, Y;
46. end proc:
47. #
48. N:=5;
49. #
50. XX, YY := RGKH(a, y0, b, N, t);
51. #
52. for i from 0 to N do
53.    print("x=", XX[i], " yrk=",YY[i], " yt=",yy(XX[i])," pogr=",abs(YY[i]-yy(XX[i])));
54. end do:
55. #
56. RGKAUTOH := proc(X0, Y0, H, B, EPS, t)
57. local h, X, Y, xx, yy, y1, pr, y2, xi, k, i, cd, xnew;
58. X[0] := X0;
59. Y[0] := Y0;
60. xx := X0;
61. yy := Y0;
62. h := H;
63. i := 1;
64. cd:=0;
65. while (xx + h <= B) and (cd =0) do
66.   xnew, y1 := RGKSTEP(xx, yy, h, t);
67.   xnew, pr := RGKSTEP(xx, yy, h/2, t);
68.   xnew, y2 := RGKSTEP(xx + h/2, pr, h/2, t);
69.   k := 0;
70.   while (k < 5) and (abs ((y1 - y2)/7) >= EPS) do
71.     k := k + 1;
72.     h := h/2;
73.     y1:= pr;
74.     xnew, pr := RGKSTEP(xx, yy, h/2, t);
75.     xnew, y2 := RGKSTEP(xx + h/2, pr, h/2, t);
76.   end do:
77.   if k < 5 then
78.     X[i] := xx + h;
79.     Y[i] := y2;
80.     xx := X[i];
81.     yy := Y[i];
82.     i := i + 1;
83.     if (abs((y1 - y2)/(7)) < EPS/16) then
84.       h := h*2;
85.     end if:
86.     if (xx + h > B) then
87.       h := B - xx;
88.       if (h = 0) then
89.         cd := 99;
90.       fi:
91.     fi:
92.   else
93.     cd := 99;
94.   end if:
95. end do:
96.
97. return X, Y, i-1 ;
98. end proc:
99. #
100. h0:=0.01;
101. a; b; yy(a);
102. #
103. Xh, Yh, ih := RGKAUTOH(a, yy(a), h0, b, 0.000001, t):
104. print ("Значение X: ",Xh[0],
105.        " Посчитанное значение в точке X: ", Yh[0],
106.        " Точное значение в X: ", evalf(yy(Xh[0])),
107.        " Шаг: ", evalf(h0),
108.        " Погрешность: ", abs(Yh[0]- evalf(yy(Xh[0])))
109.       );
110. for i from 1 to ih  do
111.    print ("Значение X: ", Xh[i],
112.           " Посчитанное значение в точке X: ", Yh[i],
113.           " Точное значение в X: ", evalf(yy(Xh[i])),
114.           " Шаг: ", evalf(Xh[i]-Xh[i-1]),
115.           " Погрешность: ", abs(Yh[i]- evalf(yy(Xh[i])))
116.          );
117. end do:
118. #
119. with(plots):
120. gr1:=pointplot({seq([Xh[i], Yh[i]], i=0..ih)}, color = red):
121. gr2:=plot(yy(x), x=a..b, color=blue, thickness=2):
122. display(gr1, gr2);
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