Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- clear all
- clc
- %RK IV
- f=@(x,y) exp(1).^(-1*x) -2*x;
- y0 = -1;
- x1(1) = 0;
- y1(1) = y0;
- n1 = 1;
- h1 = 0.01;
- ilosc_dzialan_dla_h1 = 0;
- tic
- while x1<=1
- %tic;
- x1(n1+1) = x1(n1) + h1;
- k11 = h1*f(y1(n1),x1(n1));
- k21 = h1*f(y1(n1)+(k11/2),x1(n1)+(h1/2));
- k31 = h1*f(y1(n1)+(k21/2),x1(n1)+(h1/2));
- k41 = h1*f(y1(n1)+k31,x1(n1)+h1);
- phi1 = (1/6)*(k11+2*k21+2*k31+k41);
- y1(n1+1) = y1(n1)+phi1;
- n1 = n1+1;
- ilosc_dzialan_dla_h1 = ilosc_dzialan_dla_h1 + 1;
- %toc;
- end
- ilosc_dzialan_dla_h1
- y1
- Time1 = toc
- plot(x1,y1,'r-');
- hold on
- x2(1) = 0;
- y2(1) = y0;
- n2 = 1;
- h2 = 0.001;
- ilosc_dzialan_dla_h2 = 0;
- tic
- while x2<=1
- %tic;
- x2(n2+1) = x2(n2) + h2;
- k12 = h2*f(y2(n2),x2(n2));
- k22= h2*f(y2(n2)+(k12/2),x2(n2)+(h2/2));
- k32 = h2*f(y2(n2)+(k22/2),x2(n2)+(h2/2));
- k42 = h2*f(y2(n2)+k32,x2(n2)+h2);
- phi2 = (1/6)*(k12+2*k22+2*k32+k42);
- y2(n2+1) = y2(n2)+phi2;
- n2 = n2+1;
- ilosc_dzialan_dla_h2 = ilosc_dzialan_dla_h2 + 1;
- %toc;
- end
- ilosc_dzialan_dla_h2
- y2
- Time2 = toc
- plot(x2,y2,'b-');
- hold on
- x3(1) = 0;
- y3(1) = y0;
- n3 = 1;
- h3 = 0.0001;
- ilosc_dzialan_dla_h3 = 0;
- tic
- while x3<=1
- %tic;
- x3(n3+1) = x3(n3) + h3;
- k13 = h3*f(y3(n3),x3(n3));
- k23 = h3*f(y3(n3)+(k13/2),x3(n3)+(h3/2));
- k33 = h3*f(y3(n3)+(k23/2),x3(n3)+(h3/2));
- k43 = h3*f(y3(n3)+k33,x3(n3)+h3);
- phi3 = (1/6)*(k13+2*k23+2*k33+k43);
- y3(n3+1) = y3(n3)+phi3;
- n3 = n3+1;
- ilosc_dzialan_dla_h3 = ilosc_dzialan_dla_h3 + 1;
- %toc;
- end
- ilosc_dzialan_dla_h3
- y3
- Time3 = toc
- plot(x3,y3,'g-');
- hold on
- rown1=@(x,y) exp(1).^(-x)-2*x;
- tspan=0:0.0001:1;
- [x,y]=ode45(rown1, tspan, -1);
- plot(x,y)
- 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