1 %{
  2    Composite Simpson and Gauss-Legendre Integration
  3    Octave 5.1.0
  4    John Bryan 2019
  5 %}
  6
  7 clear all;
  8 echo off;
  9 format long;
 10 % :set ft=octave
 11
 12
 13 function z=simpson(f,N,a,b)
 14    % Equation 1 implementation
 15    h=(b-a)/(2*N);
 16    sum=0;
 17    for k=[1:1:N]
 18        sum=sum+f(a+h*(2*k-2))+4*f(a+h*(2*k-1))+f(a+h*2*k);
 19    endfor
 20    z=h/3*sum;
 21    return;
 22 endfunction
 23
 24
 25 function bc=binomialcoefficient(n,k)
 26      % Returns binomial coefficient.
 27      numerator=1;
 28      for i=0:(k-1)
 29         numerator=numerator*(n-i);
 30      endfor
 31      bc=numerator/factorial(k);
 32 endfunction
 33
 34
 35 function poly0=calculatepoly(n)
 36    % Equation 7 implementation.
 37    poly0=zeros(1,n+1);
 38    for k=0:n
 39        factor1=2^n;
 40        factor2=binomialcoefficient(n,k);
 41        factor3=binomialcoefficient((n+k-1)/2,n);
 42        poly0(k+1)=factor1*factor2*factor3;
 43    endfor
 44 endfunction
 45
 46
 47 function u=gaussquad(f,N,a,b,x,w)
 48     % Equation 2 implementation
 49     c=(b-a)/2;
 50     d=(b+a)/2;
 51     sum=0;
 52     for k=[1:1:N]
 53         sum=sum+w(k)*f(c*x(k)+d);
 54     endfor
 55     u=c*sum;
 56     return;
 57 endfunction
 58
 59
 60 function y=fo(x)
 61     y=exp(x)*cos(x);
 62     return;
 63 endfunction
 64
 65
 66 function z=go(x)
 67     z=x^8-x^7;
 68     return;
 69 endfunction
 70
 71
 72 function x=initialroots(n)
 73    % Equation 3 implementation
 74    term2=1/(8*n^2);
 75    term3=1/(8*n^3);
 76    % linspace(start,end,number_of_elements)
 77    k=linspace(1,n,n)
 78    cosarg=pi*(4*k-1)/(4*n+2);
 79    % roots in descending order
 80    x=(1-term2+term3)*cos(cosarg);
 81 endfunction
 82
 83
 84 function x=newtonmethod(x,n)
 85    % Equation 4 implementation
 86    p_1=calculatepoly(n-1);
 87    p_1=fliplr(p_1);
 88    p_0=calculatepoly(n);
 89    p_0=fliplr(p_0);
 90    for i=1:n
 91       for k=1:10
 92          denominator=polyval(p_1,x(i))-(x(i)*polyval(p_0,x(i)));
 93          x(i)=x(i)-((1-x(i)^2)/n)*(polyval(p_0,x(i)))/denominator;
 94       endfor
 95    endfor
 96 endfunction
 97
 98
 99 function weights=calculateweights(n,x)
100    % Equation 5 implementation
101    P=calculatepoly(n+1);
102    % fliplr(p) to use with polyval()
103    P=fliplr(P);
104    weights=zeros(1,n);
105    for k=1:n
106       k;
107       numerator=2*(1-(x(k))^2);
108       denominator=(n+1)^2*(polyval(P,x(k)))^2;
109       weights(k)=numerator/denominator;
110    endfor
111 endfunction
112
113
114 function gausslegendre_value=gausslegendre(f,n,a,b)
115    x=initialroots(n);
116    x=newtonmethod(x,n);
117    w=calculateweights(n,x);
118    gausslegendre_value=gaussquad(f,n,a,b,x,w);  
119 endfunction
120
121
122 f=@fo;
123 N=5;
124 a=1;
125 b=4;
126 gausslegendre_value=gausslegendre(f,N,a,b)  
127 simpson_value=simpson(f,N,a,b)
128 simpson_value_N_doubled=simpson(f,2*N,a,b)
129 exact_value=(1/2*exp(4))*(cos(4)+sin(4))-(1/2*exp(1))*(cos(1)+sin(1))
130 g=@go;
131 N=6;
132 a=3;
133 b=5;
134 gausslegendre_value2=gausslegendre(g,N,a,b)  
135 simpson_value=simpson(g,N,a,b)
136 simpson_value_N_doubled=simpson(g,2*N,a,b)
137 exact_value2=((1/9)*5^9-(1/8)*5^8)-((1/9)*3^9-(1/8)*3^8)