Advertisement
ayurchyk1998

Untitled

Jan 30th, 2021
1,645
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Scilab 1.93 KB | None | 0 0
  1. clc;
  2. clear;
  3. xdel(winsid());
  4.  
  5. // Initializing variables
  6. h = 0.05;
  7. x = 3;
  8. y = 1;
  9. z = 2;
  10. xt = 0.2;
  11. yt = 0.1;
  12. zt = 0.05;
  13.  
  14. // Start values equals zero
  15. for a = 1:2000
  16.  x_st = 0;
  17.  y_st = 0;
  18.  z_st = 0;
  19.  xt_st = 0;
  20.  yt_st = 0;
  21.  zt_st = 0;
  22. end
  23.  
  24. // myval equals zero in all dimensions
  25. for i = 1:2000
  26.  for j = 1:3
  27.     myvalue = 0;
  28.  end;
  29. end
  30.  
  31. // Using RK4 method to solve equations for 2000 points
  32. for i = 1:2000
  33.  myvalue(i,1) = y + z;
  34.  myvalue(i,2) = -x + 0.5*y;
  35.  myvalue(i,3) = x^2 - z;
  36.  
  37.  kx1 = h*(y + z);
  38.  ky1 = h*(-x + 0.5*y);
  39.  kz1 = h*(x^2 - z);
  40.  
  41.  kx2 = h*((y+ky1)/2 + (z+kz1)/2);
  42.  ky2 = h*((-x+kx1)/2 + 0.5*((y+ky1)/2));
  43.  kz2 = h*(((x+kx1)^2)/2 - (z+kz1)/2);
  44.  
  45.  kx3 = h*((y+ky2)/2 + (z+kz2)/2);
  46.  ky3 = h*((-x+kx2)/2 + 0.5*((y+ky2)/2));
  47.  kz3 = h*(((x+kx2)^2)/2 - (z+kz2)/2);
  48.  
  49.  kx4 = h*((y+ky3) + (z+kz3));
  50.  ky4 = h*((-x+kx3) + 0.5*(y+ky3));
  51.  kz4 = h*((x+kx3)^2 - (z+kz3));
  52.  
  53.  xt = (kx1+2*kx2+2*kx3+kx4)/6 + xt;
  54.  yt = (ky1+2*ky2+2*ky3+ky4)/6 + yt;
  55.  zt = (kz1+2*kz2+2*kz3+kz4)/6 + zt;
  56.  
  57.  xt_st(i) = xt;
  58.  yt_st(i) = yt;
  59.  zt_st(i) = zt;
  60.  
  61.  x_st(i) = x;
  62.  y_st(i) = y;
  63.  z_st(i) = z;
  64.  
  65.  x = x + h;
  66.  y = y + h;
  67.  z = z + h;
  68.  
  69.  if i < 51
  70.      mprintf("x(t) = %f, y(t) = %f, z(t) = %f\n", xt_st(i), yt_st(i), zt_st(i));
  71.  end;
  72.  
  73. end;
  74.  
  75. plot(x_st, xt_st, 'g*-');
  76. plot(y_st, yt_st, 'b*-');
  77. plot(z_st, zt_st, 'r*-');
  78. h1 = legend(['X(t)'; 'Y(t)'; 'Z(t)']);
  79. ylabel("t");
  80.  
  81. figure();
  82. param3d(xt_st, yt_st, zt_st);
  83. h1=legend(['X(t) = f(Y(t), Z(t))']);
  84.  
  85. //graphs after compression
  86. minValue = min(xt_st);
  87. maxValue = max(yt_st);
  88. for i = 1:2000
  89.  xt_st(i) = ((xt_st(i)-minValue) * 2)/(maxValue-minValue) - 1;
  90.  yt_st(i) = ((yt_st(i)-minValue) * 2)/(maxValue-minValue) - 1;
  91.  zt_st(i) = ((zt_st(i)-minValue) * 2)/(maxValue-minValue) - 1;
  92. end
  93. figure();
  94.  
  95. plot(x_st, xt_st,'g*-');
  96. plot(y_st, yt_st,'b*-');
  97. plot(z_st, zt_st,'r*-');
  98. h1=legend(['X(t)'; 'Y(t)'; 'Z(t)']);
  99. ylabel("t");
  100.  
  101. figure();
  102. param3d(xt_st, yt_st, zt_st);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement