Leedwon

Untitled

Jun 12th, 2018
40
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.84 KB | None | 0 0
  1. clear all
  2. clc
  3.  
  4. %RK IV
  5.  
  6. f=@(x,y) exp(1).^(-1*x) -2*x;
  7. y0 = -1;
  8.  
  9. x1(1) = 0;
  10. y1(1) = y0;
  11. n1 = 1;
  12. h1 = 0.01;
  13. ilosc_dzialan_dla_h1 = 0;
  14. tic
  15. while x1<=1
  16. %tic;
  17. x1(n1+1) = x1(n1) + h1;
  18. k11 = h1*f(y1(n1),x1(n1));
  19. k21 = h1*f(y1(n1)+(k11/2),x1(n1)+(h1/2));
  20. k31 = h1*f(y1(n1)+(k21/2),x1(n1)+(h1/2));
  21. k41 = h1*f(y1(n1)+k31,x1(n1)+h1);
  22. phi1 = (1/6)*(k11+2*k21+2*k31+k41);
  23. y1(n1+1) = y1(n1)+phi1;
  24. n1 = n1+1;
  25. ilosc_dzialan_dla_h1 = ilosc_dzialan_dla_h1 + 1;
  26. %toc;
  27. end
  28. ilosc_dzialan_dla_h1
  29. y1
  30. Time1 = toc
  31. plot(x1,y1,'r-');
  32.  
  33. hold on
  34.  
  35.  
  36.  
  37. x2(1) = 0;
  38. y2(1) = y0;
  39. n2 = 1;
  40. h2 = 0.001;
  41. ilosc_dzialan_dla_h2 = 0;
  42. tic
  43. while x2<=1
  44. %tic;
  45. x2(n2+1) = x2(n2) + h2;
  46. k12 = h2*f(y2(n2),x2(n2));
  47. k22= h2*f(y2(n2)+(k12/2),x2(n2)+(h2/2));
  48. k32 = h2*f(y2(n2)+(k22/2),x2(n2)+(h2/2));
  49. k42 = h2*f(y2(n2)+k32,x2(n2)+h2);
  50. phi2 = (1/6)*(k12+2*k22+2*k32+k42);
  51. y2(n2+1) = y2(n2)+phi2;
  52. n2 = n2+1;
  53. ilosc_dzialan_dla_h2 = ilosc_dzialan_dla_h2 + 1;
  54. %toc;
  55. end
  56.  
  57. ilosc_dzialan_dla_h2
  58. y2
  59. Time2 = toc
  60. plot(x2,y2,'b-');
  61.  
  62. hold on
  63.  
  64.  
  65. x3(1) = 0;
  66. y3(1) = y0;
  67. n3 = 1;
  68. h3 = 0.0001;
  69. ilosc_dzialan_dla_h3 = 0;
  70. tic
  71. while x3<=1
  72. %tic;
  73. x3(n3+1) = x3(n3) + h3;
  74. k13 = h3*f(y3(n3),x3(n3));
  75. k23 = h3*f(y3(n3)+(k13/2),x3(n3)+(h3/2));
  76. k33 = h3*f(y3(n3)+(k23/2),x3(n3)+(h3/2));
  77. k43 = h3*f(y3(n3)+k33,x3(n3)+h3);
  78. phi3 = (1/6)*(k13+2*k23+2*k33+k43);
  79. y3(n3+1) = y3(n3)+phi3;
  80. n3 = n3+1;
  81. ilosc_dzialan_dla_h3 = ilosc_dzialan_dla_h3 + 1;
  82. %toc;
  83. end
  84.  
  85. ilosc_dzialan_dla_h3
  86. y3
  87. Time3 = toc
  88. plot(x3,y3,'g-');
  89. hold on
  90.  
  91. rown1=@(x,y) exp(1).^(-x)-2*x;
  92. tspan=0:0.0001:1;
  93. [x,y]=ode45(rown1, tspan, -1);
  94. plot(x,y)
  95.  
  96.  
  97. legend('Wykres dla kroku : h1','Wykres dla kroku : h2','Wykres dla kroku : h3','Dokładny wykres ode45');
Add Comment
Please, Sign In to add comment