Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford

This notebook demonstrates an application of discrete partial differential equations. 

Here we take a rectangular patch with polygonal  boundaries and minimize Laplace's equation, f_xx + f_yy = 0, to find the interior points. The number of iterations done to minimize Laplace's equation determines the quality of the solution.

Elements to play with:
1)  The number of points on the boundary -- the discretization of the polynomials that define the boundaries. (This is "n" below.)
2) The number of iterations. In the first cell, you can manually set "numiter"; in the manipulate cell below, Mathematica will animate this variable.
3) Try different boundary curves -- fo, f1, g0, g1

In the first cell below, you change the values manually.
In the second cell below, you can manipulate the "numiter" variable.

Figure 11.11 was created with this notebook.

See Chapter 11: Computing with Multivariate Functionsfor details.

Updated May 2008

(* Set the output folder to be the one this file lives in *)

(SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];) 

(* Define the four boundary curves *)
(* f0 and f1 are opposing; g0 and g1 are opposing *)

(f0[x_]:=8.` x (1-x);)

 (f1[x_]:=8.` x (1-x);) 

(g0[y_]:=8.` y (1-y);)

 (g1[y_]:=8.` y (1-y);) 

(* Define the discretization of the boundary curves *)
(n=10;) 

(* *********** *)
(* How many iterations will we do to satisfy the equation *)
(numiter=50;) 
(* *********** *)

(* Initialize the nxnx table with zeros *)
(F=Table[0.`,{i,1,n},{j,1,n}];) 

(* Set up a simple name for n-1 *)
(n1=n-1;) 

(* Discretize the boundary curves and load table F with values at n places *)
Do[F[[1,i]]=f0[((i-1)/n1)],{i,1,n}] ;

(Do[F[[n,i]]=f1[((i-1)/n1)],{i,1,n}];) 

(Do[F[[i,1]]=g0[((i-1)/n1)],{i,1,n}];)

 (Do[F[[i,n]]=g1[((i-1)/n1)],{i,1,n}];) 

G = F;

(* F[[1,2]] *)


(* Update the mask numiter times *)
(* The mask is minimizing Laplace's equation fxx + fyy = 0 *)

(Do[Do[F[[i,j]]=(1/(4.`))(2 F[[i+1,j]]+2 F[[i-1,j]]+2 F[[i,j-1]]+2 F[[i,j+1]]-F[[i-1,j-1]]-F[[i+1,j-1]]-F[[i+1,j+1]]-F[[i-1,j+1]]),{i,2,n-1},{j,2,n-1}],{k,1,numiter}];)



 (* F[[1,2]]  *)

(* Create a plot of the table -- function *)

(* make the surface a bit opaque *) 

picfct=ListPlot3D[F,Mesh->False,Axes->False];
 pic0= ListPlot3D[G, Mesh->True,Axes->False];
Show[picfct, pic0]



(* Output the postscript file *)
(* THis takes more time for the nb to run. Uncomment when wanted *)
(* (Export["Coons.eps",pic];) *)

(* ******
Now let's use the manipuate command to change the number of iterations -- numiter 
****** *)

createF[numI_] := (F=Table[0.`,{i,1,n},{j,1,n}];
Do[F[[1,i]]=f0[((i-1)/n1)],{i,1,n}] ;
Do[F[[n,i]]=f1[((i-1)/n1)],{i,1,n}];
Do[F[[i,1]]=g0[((i-1)/n1)],{i,1,n}];
 Do[F[[i,n]]=g1[((i-1)/n1)],{i,1,n}];
Do[Do[F[[i,j]]=(1/(4.`))(2 F[[i+1,j]]+2 F[[i-1,j]]+2 F[[i,j-1]]+2 F[[i,j+1]]-F[[i-1,j-1]]-F[[i+1,j-1]]-F[[i+1,j+1]]-F[[i-1,j+1]]),{i,2,n-1},{j,2,n-1}],{k,1,numI}]; Return[ListPlot3D[F,Mesh->False,Axes->False, PlotRange->All]];)


Manipulate[createF[numI],{{numI,0}, 0,400, 5}]