Advertisement
legrahamcracker

Composite Simpson's Rule

May 22nd, 2021
2,101
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Octave 1.31 KB | None | 0 0
  1. function [outI, outX] = composite_two_simpson (func, a, b, h)
  2. ## because of how matlab indexes vectors (1 to "k") instead of (0 to "our n")
  3. ## have to change the way we count the starting and ending points of intervals
  4. ## h represents distance between mesh points
  5. ## can input (b-a)/n in for h if you know how many subintervals you want
  6.  
  7. if h < (b-a)/2 &&  h > 0
  8.  
  9. ## number of intervals labeled k and initalize other sums
  10.   k = 1+(b-a)/h;
  11.   start_pts = 0;
  12.   end_pts = 0;
  13.  
  14. ## construct outX vector which contains all mesh points
  15.   outX = a:h:b;
  16.  
  17. ## 4f(x(2*i-1)) from 1 to n/2
  18. ## instead go from 1 to (k/2)-1
  19.   for i = 1:floor(k/2-1)
  20.     start_pts = start_pts + func(outX(2*i));
  21.     i = i + 1;
  22.  
  23.   endfor
  24.  
  25. ## 2f(x(2*i)) from 1 to (n/2)-1
  26. ## instead go from 2 to (k/2)
  27.  
  28.   for j = 2:(floor(k/2))
  29.     end_pts = end_pts + func(outX(2*j-1));
  30.     j = j + 1;
  31.  
  32.   endfor
  33.   outI = h/3*(func(a) + 4*start_pts + 2*end_pts + func(b));
  34.  
  35. ## elseif for simple simpson rule with 2 subintervals
  36. ## m is the midpoint
  37. ## outI is the approximation of the integral
  38. ## outX is the vector containing unique mesh points
  39. elseif h == (b-a)/2
  40.   m = (a+b)/2;
  41.   outI = h/3*(func(a) + 4*func((a+b)/2) + func(b));
  42.   outX = [a, m, b];
  43.  
  44. else
  45.   print("distance h is too large or a nonpositive distance")
  46.  
  47. endif
  48.  
  49. endfunction
  50.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement