fork download
  1.  
  2. %euler loop
  3. f = @(x, y) 2*x; % f(x, y) = 2*x
  4.  
  5.  
  6. h = 0.4;
  7. n = (2 - 0)/h;
  8.  
  9.  
  10. x(1) = 0;
  11. y(1) = 1;
  12.  
  13.  
  14. for i = 1:n
  15. x(i+1) = x(i) + h;
  16. y(i+1) = y(i) + h * f(x(i), y(i));
  17. end
  18.  
  19. %modified euler method
  20.  
  21. %write the function
  22.  
  23. f = @(x, y) exp(x^2);
  24.  
  25. h = 0.1;
  26. n = (2-0)/h;
  27.  
  28. x(1) = 0;
  29. y(1) = 1;
  30.  
  31. for i = 1:n
  32. x(i+1) = x(i) + h; %step size
  33.  
  34. y(i+1) = y(i) + h*f(x(i), y(i));
  35.  
  36. y(i+1) = y(i) + h * (f(x(i), y(i)) + f(x(i+1), y(i+1))) / 2;
  37.  
  38. end
  39. disp([x' y'])
  40.  
  41.  
  42. % simpson 1/3
  43.  
  44. f=@(x) exp(x.^2);
  45.  
  46. a=0;
  47. b=1;
  48.  
  49. n=12;
  50.  
  51. h=(b-a)/n;
  52.  
  53. sum=f(a)+f(b);
  54.  
  55.  
  56. for i=1:n-1
  57. a=a+h;
  58. if mod(i,2)==1
  59.  
  60. sum=sum+4*f(a);
  61. else
  62.  
  63. sum=sum+2*f(a);
  64. end
  65. end
  66. Sum=h/3*sum
  67.  
  68. fprintf('Approximate integral = %.10f\n',Sum);
  69.  
  70.  
  71. %simpson 3/8
  72.  
  73. f=@(x) exp(x.^2);
  74.  
  75.  
  76. a=0;
  77. b=1;
  78.  
  79. n=12;
  80.  
  81.  
  82. h=(a-b)/n;
  83.  
  84. sum= f(a)+f(b);
  85.  
  86.  
  87. for i=1:n-1
  88. a=a+h;
  89.  
  90. if mod(i,3)==0
  91. sum=sum +2*f(a);
  92.  
  93. else
  94. sum=sum+3*f(a);
  95. end
  96.  
  97. end
  98.  
  99.  
  100. sum=3*h*sum/8;
  101. fprintf('Approximate integral = %.10f\n', sum);
  102.  
  103.  
  104. %TROPZOIDAL e^x^2 from tropzoidal n=12
  105.  
  106.  
  107. f=@(x) exp(x.^2);
  108.  
  109.  
  110. a=0;
  111. b=1;
  112.  
  113. n=12
  114.  
  115. h=(b-a)/n;
  116.  
  117. sum=f(a)+f(b);
  118.  
  119. for i=1:n-1
  120.  
  121. a=a+h;
  122. sum=sum+2*f(a);
  123. end
  124. sum=h/2*sum
  125. fprintf('Approximate intregal: %.10f\n', sum')
  126.  
  127.  
  128. f= @(x) exp(x.^2);
  129. result=quad(f,0,1)
  130.  
  131.  
  132. %midpoint
  133.  
  134. f=@(x,y) x+y;
  135.  
  136. x(1)=0;
  137. y(1)=1;
  138.  
  139. xn=0.2;
  140. h=0.1;
  141.  
  142. n=(xn-x(1))/h;
  143.  
  144. for i=1:n
  145.  
  146.  
  147. k1=f(x(i),y(i));
  148. k2=f(x(i)+h/2,y(i)+h/2*k1);
  149. y(i+1)=y(i)+h*k2;
  150. x(i+1)=x(i)+h;
  151.  
  152. end
  153. [x' y']
  154.  
  155.  
  156. %ratson method
  157.  
  158. f=@(x,y) x+y;
  159.  
  160. x(1)=0;
  161. y(1)=1;
  162.  
  163. xn=0.2;
  164. h=0.1;
  165.  
  166. n=(xn-x(1))/h;
  167.  
  168.  
  169. for i=1:n
  170.  
  171.  
  172. k1=f(x(i),y(i));
  173. k2=f(x(i)+3*h/4,y(i)+3*h*k1/4);
  174. y(i+1)=y(i)+h*(k1+2*k2)/3;
  175. x(i+1)=x(i)+h;
  176.  
  177. end
  178. [x' y']
  179.  
  180.  
  181. %rk 4th method
  182.  
  183.  
  184. f=@(x,y) x+y;
  185.  
  186. x(1)=0;
  187. y(1)=1;
  188.  
  189. xn=0.2;
  190. h=0.1;
  191.  
  192. n=(xn-x(1))/h;
  193.  
  194. for i=1:n
  195.  
  196.  
  197. k1=f(x(i),y(i));
  198. k2=f(x(i)+h/2,y(i)+h*k1/2);
  199. k3=f(x(i)+h/2,y(i)+h*k2/2);
  200. k4=f(x(i)+h,y(i)+h*k3);
  201. y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;
  202. x(i+1)=x(i)+h;
  203.  
  204. end
  205. [x' y']
  206.  
  207.  
  208. %heuns
  209.  
  210.  
  211. f=@(x,y) x+y;
  212.  
  213. x(1)=0;
  214. y(1)=1;
  215.  
  216. xn=0.2;
  217. h=0.1;
  218.  
  219. n=(xn-x(1))/h;
  220.  
  221.  
  222. for i=1:n
  223.  
  224.  
  225. k1=f(x(i),y(i));
  226. k2=f(x(i)+h,y(i)+h*k1);
  227. y(i+1)=y(i)+h*(k1+k2)/2;
  228. x(i+1)=x(i)+h;
  229.  
  230. end
  231. [x' y']
  232.  
  233.  
  234. %graph of euler loop
  235.  
  236. f = @(x, y) exp(x^2);
  237.  
  238. h = 0.001;
  239. xn=2;
  240. n = (2 - 0) / h;
  241.  
  242. x(1) = 0;
  243. y(1) = 1;
  244.  
  245. for i = 1:n
  246. x(i+1) = x(i) + h;
  247. y(i+1) = y(i) + h * f(x(i), y(i));
  248. y(i+1)=y(i)+h*(f(x(i),y(i))+f(x(i+1),y(i+1)))/2; %modi-Euler
  249. end
  250. n
  251.  
  252. y_exact = x.^2 + 1;
  253.  
  254. plot(x, y, 'ro--', x, y_exact, 'b.--');
  255. legend('Euler Approx', 'Exact Solution');
  256.  
  257.  
  258. % weddles n=12 (n shold be multiple of 6)
  259.  
  260. f = @(x) exp(x.^2);
  261.  
  262. a = 0;
  263. b = 1;
  264. n = 12;
  265.  
  266. h = (b - a) / n;
  267.  
  268.  
  269. x = a:h:b;
  270. y = f(x);
  271.  
  272. sum = 0;
  273.  
  274. for k = 1:6:n
  275. sum = sum + (3*h/10) * ( ...
  276. y(k) + 5*y(k+1) + y(k+2) + ...
  277. 6*y(k+3) + y(k+4) + 5*y(k+5) + y(k+6) );
  278. end
  279.  
  280. sum
Success #stdin #stdout 0.65s 59656KB
stdin
Standard input is empty
stdout
    0.00000    1.00000
    0.10000    1.10050
    0.20000    1.20305
    0.30000    1.30979
    0.40000    1.42318
    0.50000    1.54606
    0.60000    1.68192
    0.70000    1.83521
    0.80000    2.01165
    0.90000    2.21887
    1.00000    2.46717
    1.10000    2.77076
    1.20000    3.14947
    1.30000    3.63148
    1.40000    4.25742
    1.50000    5.08677
    1.60000    6.20795
    1.70000    7.75441
    1.80000    9.93076
    1.90000   13.05575
    2.00000   17.63396
Sum =  1.4627
Approximate integral = 1.4626661264
Approximate integral = -1.4626836999
n =  12
sum =  1.4658
Approximate intregal: 1.4657942734
result =  1.4627
ans =

    0.00000    1.00000
    0.10000    1.11000
    0.20000    1.24205
    0.30000    1.30979
    0.40000    1.42318
    0.50000    1.54606
    0.60000    1.68192
    0.70000    1.83521
    0.80000    2.01165
    0.90000    2.21887
    1.00000    2.46717
    1.10000    2.77076
    1.20000    3.14947
    1.30000    3.63148
    1.40000    4.25742
    1.50000    5.08677
    1.60000    6.20795
    1.70000    7.75441
    1.80000    9.93076
    1.90000   13.05575
    2.00000   17.63396

ans =

    0.00000    1.00000
    0.10000    1.11000
    0.20000    1.24205
    0.30000    1.30979
    0.40000    1.42318
    0.50000    1.54606
    0.60000    1.68192
    0.70000    1.83521
    0.80000    2.01165
    0.90000    2.21887
    1.00000    2.46717
    1.10000    2.77076
    1.20000    3.14947
    1.30000    3.63148
    1.40000    4.25742
    1.50000    5.08677
    1.60000    6.20795
    1.70000    7.75441
    1.80000    9.93076
    1.90000   13.05575
    2.00000   17.63396

ans =

    0.00000    1.00000
    0.10000    1.11034
    0.20000    1.24281
    0.30000    1.30979
    0.40000    1.42318
    0.50000    1.54606
    0.60000    1.68192
    0.70000    1.83521
    0.80000    2.01165
    0.90000    2.21887
    1.00000    2.46717
    1.10000    2.77076
    1.20000    3.14947
    1.30000    3.63148
    1.40000    4.25742
    1.50000    5.08677
    1.60000    6.20795
    1.70000    7.75441
    1.80000    9.93076
    1.90000   13.05575
    2.00000   17.63396

ans =

    0.00000    1.00000
    0.10000    1.11000
    0.20000    1.24205
    0.30000    1.30979
    0.40000    1.42318
    0.50000    1.54606
    0.60000    1.68192
    0.70000    1.83521
    0.80000    2.01165
    0.90000    2.21887
    1.00000    2.46717
    1.10000    2.77076
    1.20000    3.14947
    1.30000    3.63148
    1.40000    4.25742
    1.50000    5.08677
    1.60000    6.20795
    1.70000    7.75441
    1.80000    9.93076
    1.90000   13.05575
    2.00000   17.63396

n =  2000
sum =  1.4627