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
V_bem(m) = v(1); % computing one of the (identical because of symmetry) values of du/dn
V = -pressure_exact_dr(R,k,[R,0]); % the exact normal derivative
Rel_error(m) = abs((V-V_bem(m))/V); % 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 normal velocity')
grid
disp(' ')
disp(' h/lambda Relative error')
disp([h',Rel_error'])