Problem 8.1

Nasser Abbasi

OUTPUT

» nma_problem_8_1

Solves problem 8.1 in book

Nasser Abbasi

Enter number of grid points N. (Suggest 50) :50

Enter Lx, the length in x direction (Suggest 1):1

Enter Ly, the length in y direction (Suggest 1):1

Enter PHI_0, the constant term (Suggest 1):1

average PHI= 0.249125, maxTerm=11

average PHI= 0.249557, maxTerm=21

average PHI= 0.249697, maxTerm=51

Number of terms needed to get 1 percent accuracy in solution = 174

CODE

function nma_problem_8_1()

% Solves problem 8.1 in book

% Nasser Abbasi

clear all; help nma_problem_8_1;

N = input('Enter number of grid points N. (Suggest 50) :');

Lx = input('Enter Lx, the length in x direction (Suggest 1):');

Ly = input('Enter Ly, the length in y direction (Suggest 1):');

%maxTerm = input('Enter max n to terminate sum at (Suggest 11 or 21 or 51):');

PHI_0 = input('Enter PHI_0, the constant term (Suggest 1):');

hx = Lx/(N-1);

hy = Ly/(N-1);

x = (0:N-1)*hx;

y = (0:N-1)*hy;

maxTerms = [11 21 51]; % run the sum up to these values.

for(k=1:length(maxTerms))

PHI = zeros(N);

maxTerm = maxTerms(k);

for(i = 1:N)

for(j = 1:N)

PHI(i,j) = PHI_0 * calcPHI(x(i),y(j),Lx,Ly,maxTerm);

end

end

plotit(x,y,PHI,maxTerm,N);

V(k) = sum(sum(PHI))/(N*N);

fprintf('average PHI= %f, maxTerm=%d\n',V(k),maxTerm);

end

phiChanged = V(end) - V(1);

termsChanged = maxTerms(end)-maxTerms(1);

changePerTerm = phiChanged / termsChanged;

phiChangedPercent = phiChanged*100/V(1);

phiChangedPerOnePercent = phiChanged/phiChangedPercent;

numberOfTermsNeeded = phiChangedPerOnePercent/changePerTerm;

fprintf('Number of terms needed to get 1 percent accuracy in solution = %d\n',...

round(numberOfTermsNeeded));

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotit(x,y,PHI,maxTerm,N)

PHI=PHI'; % flip PHI to make it match physical layout of x,y

figure;

mesh(x,y,PHI);

xlabel('x');

ylabel('y');

zlabel('\rho(x,y)');

title(sprintf('Laplace PDE solution, direct sum. n=%d,N=%d',maxTerm,N));

figure;

cs= contour(x,y,PHI,(0:(0.1):1));

clabel(cs);

title(sprintf('Laplace PDE solution, direct sum. n=%d,N=%d',maxTerm,N));

xlabel('x');

ylabel('y');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function PHI=calcPHI(x,y,Lx,Ly,maxTerm)

PHI=0;

for(n=1:2:maxTerm)

pi_n = pi*n;

PHI = PHI + (4/pi_n) * sin( pi_n * x / Lx) * sinh(pi_n * y/Lx) / sinh(pi_n * Ly/Lx);

end