clear
k = 2*pi; % so that the wavelength is 1
lambda = 2*pi/k; % the wavelength
R = 0.5 % so that the radius of the circle is 0.5 wavelengths
x = [R/2,R/2] % picking a position at which we will compute the pressure
x0 = [0,0] % putting the source at the centre of the circle
N = [8 16 32 64 128 256];
for m = 1:length(N)
[xx,A] = circ(N(m),R); % dividing the circle into N boundary elements
h(m) = A; % the length of each boundary element
v = bem_solve(xx,A*ones(N(m),1),x0,k,N(m)); % solving the integral equation by BEM
u_bem(m) = pressure(xx,A*ones(N(m),1),x0,k,N(m),v,x); % computing pressure at x by BEM
u = pressure_exact(R,k,x); % the exact pressure
Rel_error(m) = abs((u-u_bem(m))/u); % the relative error in the BEM pressure value at x
end
plot(h/lambda,Rel_error)
xlabel('Element length as a fraction of the wavelength')
ylabel('Relative error')
title('Error in BEM calculation of pressure')
grid
disp(' ')
disp(' h/lambda Relative error')
disp([h',Rel_error'])