1
2
3
4
5
6
7 clear all;
8 echo off;
9 format long;
10 clc;
11
12
13
14
15 function three_plots()
16 x=linspace(0,1,200);
17 y=sin(x);
18 plot(x,y)
19 area (x, y, "FaceColor", "blue")
20 title ("Area under sin(x) from x=0 to x=1");
21 set(gca, "fontsize", 8)
22 print('first.png', '-dpng', '-S600,500');
23 figure
24 x=linspace(1,3,200);
25 y=x.^5;
26 plot(x,y)
27 area (x, y, "FaceColor", "red")
28 title ("Area under x^5 from x=1 to x=3");
29 set(gca, "fontsize", 8)
30 print('second.png', '-dpng', '-S600,500');
31 x=linspace(0,50,5000);
32 y=sin(x);
33 plot(x,y)
34 area (x, y, "FaceColor", "green")
35 title ("Integral of sin(x) from x=0 to x=50");
36 set(gca, "fontsize", 8)
37 print('third.png', '-dpng', '-S600,500');
38 endfunction
39
40
41 function z=simpson(f,N,a,b)
42
43 h=(b-a)/(2*N);
44 sum=0;
45 for k=[1:1:N]
46 sum=sum+f(a+h*(2*k-2))+4*f(a+h*(2*k-1))+f(a+h*2*k);
47 endfor
48 z=h/3*sum;
49 return;
50 endfunction
51
52
53 function x=initialroots(n)
54
55 term2=1/(8*n^2);
56 term3=1/(8*n^3);
57
58 k=linspace(1,n,n);
59 cosarg=pi*(4*k-1)/(4*n+2);
60
61 x=(1-term2+term3)*cos(cosarg);
62 endfunction
63
64
65 function x=newtonmethod(x,n)
66
67 p_1=calculatepoly(n-1);
68 p_1=fliplr(p_1);
69 p_0=calculatepoly(n);
70 p_0=fliplr(p_0);
71 for i=1:n
72 diff0=1;
73 count=0;
74 maxdiff=1e-16;
75 x_initial=x(i);
76 while diff0>maxdiff
77 count=count+1;
78 temp=x(i);
79 denominator=polyval(p_1,x(i))-(x(i)*polyval(p_0,x(i)));
80 x(i)=x(i)-((1-x(i)^2)/n)*(polyval(p_0,x(i)))/denominator;
81 diff0=abs(x(i)-temp);
82 if (count > 7)
83 maxdiff=maxdiff*2;
84 x(i)=x_initial;
85 endif
86 endwhile
87 count;
88 endfor
89 x=fliplr(x);
90 endfunction
91
92
93 function bc=binomialcoefficient(n,k)
94
95 numerator=1;
96 for i=0:(k-1)
97 numerator=numerator*(n-i);
98 endfor
99 bc=numerator/factorial(k);
100 endfunction
101
102
103 function poly0=calculatepoly(n)
104
105 poly0=zeros(1,n+1);
106 for k=0:n
107 factor1=2^n;
108 factor2=binomialcoefficient(n,k);
109 factor3=binomialcoefficient((n+k-1)/2,n);
110 poly0(k+1)=factor1*factor2*factor3;
111 endfor
112 endfunction
113
114
115 function weights=calculateweights(n,x)
116
117 P=calculatepoly(n+1);
118
119 P=fliplr(P);
120 weights=zeros(1,n);
121 for k=1:n
122 k;
123 numerator=2*(1-(x(k))^2);
124 denominator=(n+1)^2*(polyval(P,x(k)))^2;
125 weights(k)=numerator/denominator;
126 sumweights=sum(weights);
127 endfor
128 endfunction
129
130
131 function u=gaussquad(f,N,a,b,x,w)
132
133 c=(b-a)/2;
134 d=(b+a)/2;
135 sum=0;
136 for k=[1:1:N]
137 w(k)*f(c*x(k)+d);
138 sum=sum+w(k)*f(c*x(k)+d);
139 endfor
140 u=c*sum;
141 return;
142 endfunction
143
144
145 function gausslegendre_value=gausslegendre(f,n,a,b)
146 x=initialroots(n);
147 x=newtonmethod(x,n);
148 w=calculateweights(n,x);
149 gausslegendre_value=gaussquad(f,n,a,b,x,w);
150 endfunction
151
152
153 function compositeGL=CompositeLG(a,b,f,n,interval)
154 start=a;
155 total=0;
156 while start<b
157 gausslegendre_value=gausslegendre(f,n,start,start+interval);
158 total=total+gausslegendre_value;
159 start=start+interval;
160 endwhile
161 compositeGL=total;
162 endfunction
163
164
165 function [interval_result,error0] = adaptive_simpson(f,a,b,abstol)
166
167
168
169
170
171 h1=(b-a)/2;
172 h2 = h1/2;
173 s1 = h1/3*(f(a)+4*f(a+h1)+f(b));;
174 s2 = h2/3*(f(a)+4*f(a+h2)+2*f(a+2*h2)+4*f(a+3*h2)+f(b));
175 error0 = abs(s1-s2)/15;
176 if (error0 < abstol)
177 interval_result = s2;
178 else
179 [ltotal, error_ltotal] = adaptive_simpson(f,a,(a+b)/2,abstol/2);
180 [rtotal, error_rtotal] = adaptive_simpson(f,(a+b)/2,b,abstol/2);
181 interval_result= ltotal+rtotal;
182 error0 = error_ltotal+error_rtotal;
183 endif
184 endfunction
185
186
187 function test()
188 f=@(x) sin(x);
189 a=0;
190 b=1;
191 n=6;
192 gausslegendre_value=gausslegendre(f,n,a,b);
193 exact_value=-cos(1)+cos(0);
194 l10difference=abs(gausslegendre_value-exact_value);
195 simpson_value=simpson(f,100*n,a,b);
196 s1difference=abs(simpson_value-exact_value);
197 fprintf('\n sin(x), 0:1: exact result: %1.15e \n',exact_value);
198 fprintf('\n sin(x), 0:1: GL result: %1.15e error %1.15e\n',gausslegendre_value,l10difference);
199 fprintf('\n sin(x), 0:1 : CSR result %1.15e error %1.15e\n',simpson_value,s1difference);
200 f=@(x) x^5;
201 a=1;
202 b=3;
203 exact_value=(1/6)*(3^6-1^6);
204 gausslegendre_value=gausslegendre(f,n,a,b);
205 l2difference=abs(gausslegendre_value-exact_value);
206 simpson_value=simpson(f,100*n,a,b);
207 s2difference=abs(simpson_value-exact_value);
208 fprintf('\n x^5, 1:3: exact result: %1.15e \n',exact_value);
209 fprintf('\n x^5, 1:3: GL result: %1.15e error %1.15e\n',gausslegendre_value,l2difference);
210 fprintf('\n x^5, 1:3 : CSR result %1.15e error %1.15e\n',simpson_value,s2difference);
211 f=@(x) sin(x);
212 a=0;
213 b=50;
214 n=6;
215 interval=1;
216 compositeGL=CompositeLG(a,b,f,n,interval);
217 exact_value=-cos(b)+cos(a);
218 l3difference=abs(compositeGL-exact_value);
219 n=12500;
220 simpson_value=simpson(f,n,a,b);
221 s3difference=abs(simpson_value-exact_value);
222 fprintf('\n sin(x), 0:50: exact result: %1.15e \n',exact_value);
223 fprintf('\n sin(x), 0:50: composite GL result: %1.15e error %1.15e\n',compositeGL,l3difference);
224 fprintf('\n sin(x), 0:50 : CSR result %1.15e error %1.15e\n',simpson_value,s3difference);
225 abstol = 10^(-15);
226 [result,error0] = adaptive_simpson(f,a,b,abstol);
227 s4difference=abs(result-exact_value);
228 fprintf("\nsin(x), 0:50 : Adaptive Simpson result %1.15e error %1.15e\n\n",result,s4difference);
229 endfunction
230
231 test()