% laplace_relax generates a solution to Laplace's equation with % specified boundary conditions using a simple relaxation method. % The scalar potential phi(x,y) is represented as the matrix phi(i,j). % A new value for phi(i,j) is given by averaging the neighboring points: % phi_new = ( phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1) )/4 % Programmer: Dick Furnstahl furnstahl.1@osu.edu % Revision history: % 12-May-2005 -- original version for Physics 263 close all; clear; % this closes plots and sets variables to zero % Set up the x-y grid with "gridsize" number of points in each direction xmin = 0; xmax = 1; % set the minimum and maximum values of x ymin = 0; ymax = 1; % set the minimum and maximum values of y gridsize = 20; [X Y] = meshgrid( linspace(xmin,xmax,gridsize), linspace(ymin,ymax,gridsize) ); dx = (xmax-xmin)/(gridsize-1); % spacing between X points dy = (ymax-ymin)/(gridsize-1); % spacing between Y points % Set up boundary conditions and initial values phi = zeros(gridsize,gridsize); % set phi to all zeros phi(1,1:gridsize) = 100.; % set the nonzero boundary condition relax = 1; % parameter controlling mix of old and new phi values iterMax = 50; % how many iterations to do % Do the iterations for iter = 1:iterMax % step through the interior of the lattice, with each point labeled i,j for i = 2:gridsize-1 for j = 2:gridsize-1 phi_new = 0.25*( phi(i+1,j) + phi(i-1,j) + ... phi(i,j+1) + phi(i,j-1) ); % new phi via average phi(i,j) = relax*phi_new + (1-relax)*phi(i,j); % set current phi end end figure(1); % plot the electric potential phi at each iteration pcolor(X,Y,phi); colorbar; title(['iteraction #', int2str(iter)]); % put the iteration # in the title xlabel('x-axis'); ylabel('y-axis'); % add labels pause(.1); % pause for this many seconds end % Now calculate the electric field from minus the gradient of phi [Ex Ey] = gradient(-phi,dx,dy); pause; % wait for a return to show the other plots figure(2); % figure 2 is a surface plot of the electric potential phi surf(X,Y,phi); colorbar; xlabel('x-axis'); ylabel('y-axis'); title('Electric Potential Phi'); figure(3); % figure 3 shows the electric potential AND electric field num_contours = 6; contourf(X,Y,phi,num_contours); colorbar; % make the plot axis equal; % make the axes equal lengths hold on; % hold the current plot scale = 5; quiver(X,Y,Ex,Ey,scale); % now add the electric field xlabel('x-axis'); ylabel('y-axis'); title('Electric Field'); hold off; % release the hold