Partial Differential Equations

One Dimension - Implementation 4.1

The function below makes the one-dimensional finite difference matrix (retD) and right-hand side (the combination of source minus boundaries) retsmb.  You provide p, q, and s, the function appearing in the general form (4.38): ((d^2)/(dx^2))f(x) + p(x) (d/dx)f(x) + q(x) f(x) = s(x).  The left endpoint is at 0, the right at xf, and we make a grid with Ns internal points - the values of f at 0 and xf must be provided (called f0 and fNp1).

MakeDifference[p_, q_, s_, f0_, fNp1_, xf_, Ns_] := Module[{retD, retsmb,index,dx,x},
retD = Table[0.0,{i,1,Ns},{j,1,Ns}];
retsmb = Table[0.0,{i,1,Ns}];
dx = xf/(Ns + 1);
x = dx;
retD[[1,1]] = (q[x] - 2/dx^2);
retD[[1,2]] = (p[x]/(2 dx) + 1/dx^2);
retsmb[[1]] = s[x] - (-p[x]/(2 dx) + 1/dx^2) f0;

For[index = 2, index < Ns, index=index+1,
x = x + dx;
retD[[index,index-1]] = (-p[x]/(2 dx) + 1/dx^2);
retD[[index,index]] = (q[x] - 2/dx^2);
retD[[index,index+1]] = (p[x]/(2 dx) + 1/dx^2);
retsmb[[index]] = s[x];
];
x = x + dx;
retD[[Ns,Ns-1]] = (-p[x]/(2 dx) + 1/dx^2);
retD[[Ns,Ns]] = (q[x]- 2/dx^2);
retsmb[[Ns]] = s[x] - (p[x]/(2 dx) + 1/dx^2)fNp1;
Return[N[{retD,retsmb}]];
]

Example 1

Choose the coefficients appropriate for the ODE: ((d^2)/(dx^2))f(x) + (((3 π))^2) f(x) = 0:

p[x_] := 0

q[x_] := (3 Pi)^2

s[x_] := 0

Make the difference matrix and rhs vector for this problem

{DD,bb} = MakeDifference[p,q,s,1,-1,1,100];

invert (using the built-in function Inverse, for now)

f = Inverse[DD].bb;

G0 = ListPlot[f]

We can make the exact solution and plot the residual:

dx = N[1/101];
ef = Table[Cos[3 Pi x],{x,dx,1-dx,dx}];
ListPlot[ef-f]

Example 2 - Figure 4.1

Set up ODE coefficients for the damped driven case (4.48)

p[x_] := 1

q[x_] := (2 Pi)^2

s[x_] :=  Cos[2 Pi x] Sin[10 Pi x]

Generate matrix/vector pair and invert (using the built-in LinearSolve function, this time, just to mix things up) to get the numerical solution on a grid with 300 points:

Ns = 200;
{DD,bb} = MakeDifference[p,q,s,0,0.,1,Ns];
f= LinearSolve[DD,bb];
dx = 1/(Ns+1);
fplot = Table[{j dx, f[[j]]},{j,1,Length[f]}];
G1 = ListPlot[fplot,PlotStyle->GrayLevel[0.75]]

In this case, we know the solution - it’s not pretty:

exact[s_] := ((E^(-s/2)) (-(E^(s/2)) (2 (9+1225 (π^2)) Cos[8 π s]+3 (4+225 (π^2)) Cos[12 π s]+5 π (3 (9+1225 (π^2)) Sin[8 π s]+7 (4+225 (π^2)) Sin[12 π s]))-5 (6+625 (π^2)) Csc[(1/2) (Sqrt[-1+16 (π^2)])] (Sin[(1/2) (Sqrt[-1+16 (π^2)]) (-1+s)]-(Sqrt[E]) Sin[(1/2) (Sqrt[-1+16 (π^2)]) s])))/(8 π (4+225 (π^2)) (9+1225 (π^2)))

G0 =Plot[exact[s],{s,0,1}, PlotStyle->GrayLevel[0]]

We can put the points on top of the continuous curve:

Show[G0,G1]

Example 3 - Figure 4.2

Set up functions for the third example

p[x_] := 1

q[x_] := (3 Pi)^2

s[x_] :=  Cos[Pi x/2]

There’s a “magic” value for sf,

Ns = 1000;
sf = 25.01167873217457;

{DD,bb} = MakeDifference[p,q,s,-1/2,0.,sf,Ns];
f= LinearSolve[DD,bb];
dx = sf/(Ns+1);
fplot=Table[{j dx, f[[j]]},{j,1,Length[f]}];
G1 = ListPlot[fplot, PlotRange-> All,PlotStyle->GrayLevel[.75]]

exact[s_] := 0.00002631907568456862` (E^(-s/2)) (-19437.45062349339` Cos[(1/2) (Sqrt[-1+36 (π^2)]) s]+439.822971502571` (E^(s/2)) Cos[((π s)/2)] ((Cos[(1/2) (Sqrt[-1+36 (π^2)]) s])^2)+8.` (E^(s/2)) ((Cos[(1/2) (Sqrt[-1+36 (π^2)]) s])^2) Sin[((π s)/2)]-1033.9781949542694` Sin[(1/2) (Sqrt[-1+36 (π^2)]) s]+439.822971502571` (E^(s/2)) Cos[((π s)/2)] ((Sin[(1/2) (Sqrt[-1+36 (π^2)]) s])^2)+8.` (E^(s/2)) Sin[((π s)/2)] ((Sin[(1/2) (Sqrt[-1+36 (π^2)]) s])^2))

G0 =Plot[exact[s] ,{s,0,sf},PlotRange->All,PlotStyle->GrayLevel[0]]

Show[G0,G1]

Difference Matrices in Two Dimensions (Cartesian)

For a uniform Cartesian grid, we have the following difference matrices (from Section 4.3):

DXmat[Nx_, Ny_, dx_, dy_, left_, right_] := Module[{g,retmat,bside, index, indey,i,j,k},
g[n_,m_] := (m - 1)*Nx + n;

retmat = Table[0.0,{i,1,Nx Ny},{j,1,Nx Ny}];
bside = Table[0.0,{i,1,Nx Ny}];

For[indey = 1, indey ≤ Ny, indey = indey+1,
j = g[1,indey];
k = g[2,indey];

retmat[[j,j ]] += -2.0/dx^2;
retmat[[j,k ]] += 1.0/dx^2;
bside[[ j ]] += -left[indey dy]/dx^2;

For[index = 2, index ≤ Nx-1, index = index+1,
i = g[index-1,indey];
j = g[index,indey];
k = g[index + 1,indey];
retmat[[j,i]] += 1.0/dx^2;
retmat[[j,j]] += -2.0/dx^2;
retmat[[j,k]] += 1.0/dx^2;
];
i = g[Nx-1,indey];
j = g[Nx,indey];
retmat[[j,i ]] += 1.0/dx^2;
retmat[[ j,j]] += -2.0/dx^2;
bside[[ j]] += -right[indey dy]/dx^2;
];
Return[{retmat, bside}];
]

DYmat[Nx_, Ny_, dx_, dy_, top_, bottom_] := Module[{g,retmat,bside, index, indey,i,j,k},
g[n_,m_] := (m - 1)*Nx + n;

retmat = Table[0.0,{i,1,Nx Ny},{j,1,Nx Ny}];
bside = Table[0.0,{i,1,Nx Ny}];

For[index = 1, index ≤ Nx, index = index+1,
j = g[index,1];
k = g[index, 2];

retmat[[j,j ]] += -2.0/dy^2;
retmat[[j,k ]] += 1.0/dy^2;
bside[[ j ]] += -bottom[index dx]/dy^2;

For[indey = 2, indey ≤ Ny-1, indey = indey+1,
i = g[index,indey-1];
j = g[index,indey];
k = g[index ,indey+1];
retmat[[j,i]] += 1.0/dy^2;
retmat[[j,j]] += -2.0/dy^2;
retmat[[j,k]] += 1.0/dy^2;
];
i = g[index, Ny-1];
j = g[index, Ny];
retmat[[j,i ]] += 1.0/dy^2;
retmat[[ j,j]] += -2.0/dy^2;
bside[[ j]] += -top[index dx]/dy^2;
];
Return[{retmat, bside}];
]

Notice that these functions return the appropriate x and y difference matrices, together with the boundary vectors for the x and y directions (hence the need to provide the left and right boundary functions for DXmat, top and bottom ones for DYmat).  Finally, we need a function that makes the source portion of the right-hand side - this is the explicit source, not the effective source that comes from the boundary conditions (those are constructed inside DXmat and DYmat).

MakeSource[ Nx_, Ny_, dx_, dy_, source_] := Module[{svec, g,j, index, indey},
g[n_,m_] := (m - 1)*Nx + n;

svec = Table[0.0,{i,1,Nx Ny}];
For[index = 1, index ≤ Nx, index=index+1,
For[indey = 1, indey ≤ Ny, indey = indey + 1,
j = g[index,indey];
svec[[j ]] = source[ index dx, indey dy];
];
];
Return[svec];
]

Example Nx = Ny = 3

This three by three case is included to show the structure of the two matrices, DX and DY, and the form of the rhs vector (with both sources and boundary conditions included).

RHS[x_,y_] := cv[x,y]

top[x_] := T[x]

bottom[x_] := B[x]

left[y_] := L[y]

right[y_] := R[y]

Nx = 3;
Ny = 3;
dx = 1/(Nx + 1);
dy = 1/(Ny+1);

{DX, bx} = DXmat[Nx, Ny, dx, dy, left, right] ;
{DY, by} = DYmat[Nx, Ny, dx, dy, top, bottom];
bs = MakeSource[ Nx, Ny, dx, dy,RHS] ;
rhsvec = bx + by + bs;

MatrixForm[dx^2 DX]

MatrixForm[dy^2 DY]

MatrixForm[Chop[rhsvec]]

Laplace Problem - Figure 4.5

Setup

For the Laplace problem, there are no physical sources in our domain, so we set:

RHS[x_,y_] := 0

We can set boundary conditions around the edges - let's ground three of the edges and put a sinusoidal potential on the bottom:

top[x_] := 0

bottom[x_] := Sin[Pi x]

left[y_] := 0

right[y_] := 0

Below, we take a 25 by 25 grid with the boundaries on 0 and 1 (so the grid indices go from 1 to NX corresponding to the points dx to 1-dx).

Nx= 25;
Ny= 25;
dx = 1/(Nx + 1);
dy = 1/(Ny+1);

We are ready to make the source vector and operator matrix:

{DX, bx} = DXmat[Nx, Ny, dx, dy, left, right] ;
{DY, by} = DYmat[Nx, Ny, dx, dy, top, bottom];
bs = MakeSource[ Nx, Ny, dx, dy,RHS] ;
rhsvec = bx + by + bs;
D2mat = DX + DY;

Again - we can use whatever matrix solving routine we want, we'll use the built-in LinearSolve command (we’ll come back to this matrix inversion problem in Chapter 10):

fxy = LinearSolve[D2mat, rhsvec];

Visualization

The contour plot is probably easiest for looking at solutions.  To get the appropriate value out of our solution vector, we’ll need to define the function g:

g[n_,m_] := (m - 1)*Nx + n

Datatable= Table[{ii dx, jj dy, fxy[[ g[ii,jj] ]]},{ii,1, Nx}, {jj,1,Ny}];
Datatable = Partition[Flatten[Datatable],3];

Ga = ListContourPlot[Datatable,Contours-> 20,ColorFunction-> "GrayTones"]

Now for the exact solution:

Vsol[x_,y_] := (1/4) Sin[Pi x]((1 - Coth[Pi]) Exp[Pi y] + (1 + Coth[Pi]) Exp[-Pi y])

Datatable= Table[{ii dx, jj dy, Vsol[ii dx, jj dy]},{ii,1, Nx}, {jj,1,Ny}];
Datatable = Partition[Flatten[Datatable],3];

Gb = ListContourPlot[Datatable,Contours-> 20, ColorFunction-> "GrayTones"]

Grout= Show[GraphicsRow[{Ga,Gb}]]

Pretty good approximation (at least by eye).

Poisson Problems

Parallel Plates - Figure 4.6

With sources in the domain, the “RHS” of Poisson’s problem is non-trivial.  Below, we define parallel plates:

RHS[x_,y_] := Module[{retval, rho},
rho = 1;
retval = 0;
If[y == Round[Ny/4] dy && (x ≤(Nx -  Round[Nx/5]+1) dx && x ≥ (Round[Nx/5] dx)) , 
retval = -rho;
];
If[y == (Ny - Round[Ny/4]+1) dy && (x ≤(Nx -  Round[Nx/5]+1) dx && x ≥ (Round[Nx/5] dx)) , 
retval = rho;
];
Return[-retval];
]

We can ground the edges of our box:

top[x_] := 0

bottom[x_] := 0

left[y_] := 0

right[y_] := 0

Grid definition:

Nx =40;
Ny = 40;
dx = 1/(Nx+ 1);
dy = 1/(Ny+1);

{DX, bx} = DXmat[Nx, Ny, dx, dy, left, right] ;
{DY, by} = DYmat[Nx, Ny, dx, dy, top, bottom];
bs = MakeSource[ Nx, Ny, dx, dy,RHS] ;
rhsvec = bx + by + bs;
D2mat = DX + DY;

Solve:

fxy = LinearSolve[D2mat, rhsvec];

and make the contour plot - we can define a function that generates these:

PlotPoisson[invec_, Nx_,dx_, Ny_, dy_,ncnts_:20] := Module[{g, Datatable,retplot},
g[n_,m_] := (m - 1)*Nx + n;
Datatable= Table[{ii dx, jj dy, invec[[ g[ii,jj] ]]},{ii,1, Nx}, {jj,1,Ny}];
Datatable = Partition[Flatten[Datatable],3];
 retplot = ListContourPlot[Datatable,Contours-> ncnts, ColorFunction->"GrayTones"];
Return[retplot];
]

Ga = PlotPoisson[fxy, Nx, dx, Ny, dy]

We can also visualize the potential by plotting it in three dimensions:

PlotPoisson3D[invec_, Nx_,dx_, Ny_, dy_] := Module[{g, Datatable,retplot},
g[n_,m_] := (m - 1)*Nx + n;
Datatable= Table[{ii dx, jj dy, invec[[ g[ii,jj] ]]},{ii,1, Nx}, {jj,1,Ny}];
Datatable = Partition[Flatten[Datatable],3];
 retplot = ListPlot3D[Datatable, ColorFunction->"GrayTones"];
Return[retplot];
]

PlotPoisson3D[fxy,Nx,dx,Ny,dy]

Electric Field Figure 4.6

We can get the electric field from the potential - the function “GradV” below generates finite difference approximations to the gradient at all grid locations (to avoid special-cases at the boundaries, we go from 2 -> N-1).  Since the electric field is the negative gradient, “sign” is included to multiply by -1 if desired.

GradV[Vxy_,Nx_, Ny_, dx_,dy_,sign_] := Module[{g, gradout, i, j,gv,ij,hp1,hm1,vp1,vm1},
g[n_,m_] := (m - 1)*Nx + n;

gradout = {};
For[i = 2, i ≤ Nx - 1, i=i+1,
For[j = 2, j≤ Ny - 1, j = j + 1,
hp1 = g[i+1,j];
hm1 = g[i-1,j];
vp1 = g[i,j+1];
vm1 = g[i,j-1];
gv ={ (Vxy[[hp1]] - Vxy[[hm1]])/(2 dx),(Vxy[[vp1]] - Vxy[[vm1]])/(2 dy)};
gradout = Join[gradout,{{{i dx, j dy},sign gv}}];
];
];
Return[gradout];
]

Efield = N[GradV[fxy,Nx, Ny,dx,dy,-1]];

GE = ListVectorPlot[Efield, VectorStyle->GrayLevel[0.5]]

Grout = GraphicsRow[{Ga,GE}]

Box - Figure 4.10

Ground the edges of our domain:

top[x_] := 0

bottom[x_] := 0

left[y_] := 0

right[y_] := 0

Grid definition:

Nx =50;
Ny = 50;
dx = 1/(Nx + 1);
dy = 1/(Ny+1);

Below, we define a box that goes around the domain, five gridpoints in:

RHS[x_,y_] := Module[{retval,rho},
rho = -1000.0;
retval = 0;
If[y == 5 dy &&  (x ≥ 5 dx && x ≤ 46 dx) ,  retval = rho;];
If[y == 46 dy &&  (x ≥ 5 dx && x ≤ 46 dx) ,  retval = rho;];
If[x == 5 dx &&  (y ≥ 5 dy && y ≤ 46 dy) ,  retval = rho;];
If[x == 46 dx &&  (y ≥ 5 dy && y ≤ 46 dy) ,  retval = rho;];
Return[retval];
]

{DX, bx} = DXmat[Nx, Ny, dx, dy, left, right] ;
{DY, by} = DYmat[Nx, Ny, dx, dy, top, bottom];
bs = MakeSource[ Nx, Ny, dx, dy,RHS] ;
rhsvec = bx + by + bs;
D2mat = DX + DY;

fxy = LinearSolve[D2mat,rhsvec];

g[n_,m_] := (m - 1)*Nx + n

We have to stay away from the edges here to avoid washing out the solution:

Datatable= Table[{ii dx, jj dy, fxy[[ g[ii,jj] ]]},{ii,5, Nx-5}, {jj,5,Ny-5}];
Datatable = Partition[Flatten[Datatable],3];

Ga = ListContourPlot[Datatable,ColorFunction->"GrayTones"]

Efield = N[GradV[fxy,Nx,Ny,dx,dy,-1]];

Efr = Select[Efield,#[[1,1]] < .8 && #[[1,1]] >.2 && #[[1,2]]<.8 &&#[[1,2]] >.2&];

GE = ListVectorPlot[Efr, PlotRange-> {{.15,.85},{.15,.85}},VectorStyle->GrayLevel[0.5]]

Grout = Show[GraphicsGrid[{{Ga,GE}}]]

Exact Solution - Figure 4.9

V[x_,y_] := 2 + y ArcTan[(x-1)/y] - y ArcTan[(x+1)/y] + (x-1) Log[Sqrt[(x-1)^2 + y^2]] - (x+1) Log[Sqrt[(x+1)^2 + y^2]]

Vb[x_,y_] := V[x,y+1] + V[y,1-x] + V[-x,1-y] + V[-y,x+1]

Ga = ContourPlot[Vb[x,y],{x,-1,1},{y,-1,1}, ColorFunction->"GrayTones"]

Arbitrary Grounded Boundary - Figure 4.11

This function generates the information described in the “Arbitrary Grounded Boundary” subsection - we’ll test it and generate the content of (4.55)

HashMat[inmat_] := Module[{hvec, outmat,i,j,left,right,top,bottom,index,count,NN},
NN = Length[inmat];
outmat = Table[0.0,{i,1,NN},{j,1,NN}];
count = Length[Select[Flatten[inmat],#== 1&]];
hvec = Table[0,{j,1,count}];
index = 1;
Off[Part::partd];
Off[Part::partw];
For[i = NN, i ≥  1, i = i-1,
For[j = 1, j ≤ NN, j = j+1,
If[inmat[[i,j]] == 1,
If[inmat[[i+1,j]] ==0 || i == NN,
bottom = 0.0;
,
bottom = {j,NN-(i+1)+1};
];
If[inmat[[i-1,j]] == 0 || i - 1 == 0,
top = 0.0;
,
top = {j, NN-(i-1) + 1};
];
If[inmat[[i,j+1]] == 0 || j == NN,
right = 0.0;
,
right = {j+1, NN-i+1};
];
If[inmat[[i,j-1]] == 0 || j -1 == 0,
left = 0.0;
,
left = {j-1,NN-i+1};
];

hvec[[index]] = {{j, NN-i+1},left,right,top,bottom};
outmat[[i,j]] = index;
index = index+1;
];
];
];
On[Part::partd];
On[Part::partw];
Return[{hvec, outmat}];
]

Now for the test case, compare with (4.55):

test = {{0,0,0,0,0},{0,0,1,0,0},{0,1,1,1,0},{0,0,1,1,0},{0,0,0,0,0}};

{hv,hm} = HashMat[test]

hv[[5]]

MatrixForm[hm]

This generates the finite difference matrix given the ordering defined by hmat and the associated vector.

cMakeD2[vec_,hmat_, dx_, dy_] := Module[{g,retmat,index,NN,n,m,left,right,top,bottom},
left = 2;
right= 3;
top = 4;
bottom = 5;
NN = Length[vec];

g[n_,m_] := hmat[[Length[hmat]-(m-1),n]];

retmat = Table[0.0,{i,1,NN},{j,1,NN}];
For[index = 1, index ≤ NN, index=index+1,
retmat[[index,index]] = -2.0/dx^2 -2.0/dy^2;
If[Length[vec[[index,left]]] ==2,
{n,m} = vec[[index,left]];
retmat[[ index,g[n,m] ]] = 1.0/dx^2;
];
If[Length[vec[[index,right]]] ==2,
{n,m} = vec[[index,right]];
retmat[[ index,g[n,m] ]] = 1.0/dx^2;
];
If[Length[vec[[index,top]]] ==2,
{n,m} = vec[[index,top]];
retmat[[ index,g[n,m] ]] = 1.0/dy^2;
];
If[Length[vec[[index,bottom]]] ==2,
{n,m} = vec[[index,bottom]];
retmat[[ index,g[n,m] ]] = 1.0/dy^2;
];
];
Return[retmat];
]

We need to modify the source function to refer to the relevant locations.

cMakeSource[ hvec_, dx_, dy_, source_] := Module[{svec,index,NN,n,m},
NN = Length[hvec];
svec = Table[0.0,{i,1,NN}];
For[index = 1, index ≤ NN, index=index+1,
{n,m} = hvec[[index,1]];
svec[[index]] = source[n dx, m dy];
];
Return[svec];
]

The function below generates the matrix appropriate for the conducting plate with a hole in the center:

MakeMat[NN_] := Module[{retmat,ix,jx},
retmat = Table[1,{i,1,NN},{j,1,NN}];
ix = NN/2;
For[jx = 1, jx ≤ NN/2-5, jx = jx+1,
retmat[[ix,jx]] = 0;
];
For[jx = NN/2+5, jx ≤ NN, jx = jx+1,
retmat[[ix,jx]] = 0;
];
Return[retmat];
]

MatrixPlot[MakeMat[60]]

Set up the grid:

Nx = 60;
Ny = 60;
dx = 1.0/(Nx+1);
dy = 1.0/(Ny+1);
Mat = MakeMat[Nx];
{hvec,hmat} = HashMat[Mat];

define the source:

source[x_,y_] := If[x/dx == 30&&  y/dy==10 , 100,0]

make the matrices and solve.

D2 = cMakeD2[hvec,hmat, dx,dy];
svec = cMakeSource[hvec,dx,dy,source];
fxy = LinearSolve[D2,svec];

Datatable= Table[{hvec[[i,1,1]] dx, hvec[[i,1,2]]dy , fxy[[i]]},{i,1,Length[hvec]}];
Datatable = Partition[Flatten[Datatable],3];
G0 = ListContourPlot[Datatable,Contours->100,ColorFunction->"GrayTones"]

Putting in the conducting boundary, we have Figure 4.11.

G1 = Show[Graphics[{Line[{{0,.5},{.5-5dx,.5}}], Line[{{.5 + 5 dx,.5},{1,.5}}]}]];

Grout = Show[G0,G1]