Advertisement
Leedwon

Untitled

Jun 11th, 2018
67
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.26 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. rown1=@(x,y) exp(1).^(-x)-2*x;
  34. tspan=0:0.01:1;
  35. [x,y]=ode45(rown1, tspan, -1);
  36. plot(x,y,'r+')
  37. %tu gdzies jakies odejmowanie z tymi bledami
  38.  
  39. hold on
  40.  
  41.  
  42.  
  43. x2(1) = 0;
  44. y2(1) = y0;
  45. n2 = 1;
  46. h2 = 0.001;
  47. ilosc_dzialan_dla_h2 = 0;
  48. tic
  49. while x2<=1
  50. %tic;
  51. x2(n2+1) = x2(n2) + h2;
  52. k12 = h2*f(y2(n2),x2(n2));
  53. k22= h2*f(y2(n2)+(k12/2),x2(n2)+(h2/2));
  54. k32 = h2*f(y2(n2)+(k22/2),x2(n2)+(h2/2));
  55. k42 = h2*f(y2(n2)+k32,x2(n2)+h2);
  56. phi2 = (1/6)*(k12+2*k22+2*k32+k42);
  57. y2(n2+1) = y2(n2)+phi2;
  58. n2 = n2+1;
  59. ilosc_dzialan_dla_h2 = ilosc_dzialan_dla_h2 + 1;
  60. %toc;
  61. end
  62.  
  63. ilosc_dzialan_dla_h2
  64. y2
  65. Time2 = toc
  66. plot(x2,y2,'b-');
  67.  
  68. rown2=@(x,y) exp(1).^(-x)-2*x;
  69. tspan2=0:0.001:1;
  70. [x,y]=ode45(rown2, tspan2, -1);
  71. plot(x,y,'y-o')
  72. %tu gdzies jakies odejmowanie z tymi bledami
  73.  
  74. hold on
  75.  
  76.  
  77. x3(1) = 0;
  78. y3(1) = y0;
  79. n3 = 1;
  80. h3 = 0.0001;
  81. ilosc_dzialan_dla_h3 = 0;
  82. tic
  83. while x3=1
  84. %tic;
  85. x3(n3+1) = x3(n3) + h3;
  86. k13 = h3*f(y3(n3),x3(n3));
  87. k23 = h3*f(y3(n3)+(k13/2),x3(n3)+(h3/2));
  88. k33 = h3*f(y3(n3)+(k23/2),x3(n3)+(h3/2));
  89. k43 = h3*f(y3(n3)+k33,x3(n3)+h3);
  90. phi3 = (1/6)*(k13+2*k23+2*k33+k43);
  91. y3(n3+1) = y3(n3)+phi3;
  92. n3 = n3+1;
  93. ilosc_dzialan_dla_h3 = ilosc_dzialan_dla_h3 + 1;
  94. %toc;
  95. end
  96.  
  97. ilosc_dzialan_dla_h3
  98. y3
  99. Time3 = toc
  100. plot(x3,y3,'g-');
  101. hold on
  102.  
  103. rown1=@(x,y) exp(1).^(-x)-2*x;
  104. tspan3=0:0.0001:1;
  105. [x,y]=ode45(rown1, tspan3, -1);
  106. plot(x,y,'-o')
  107. %tu gdzies jakies odejmowanie z tymi bledami
  108.  
  109.  
  110. %legend('Wykres dla kroku : h1','Wykres dla kroku : h2','Wykres dla kroku : h3','Dokładny wykres ode23');
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement