SHARE
TWEET

Untitled

a guest Mar 26th, 2019 85 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. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top