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)