COMPLEX ANALYSIS: Mathematica 3.0 Notebooks
(c) John H. Mathews, and
Russell W. Howell, 1998
Complimentary software to accompany our textbook:

COMPLEX ANALYSIS: for Mathematics and Engineering, 
3rd Edition,  1997,  ISBN: 0-7637-0270-6
Jones & Bartlett Publishers, Inc.
40 Tall Pine Drive, Sudbury, MA 01776
Tele.  (800) 832-0034, FAX:  (508)  443-8000
E-mail:  mkt@jbpub.com,     http://www.jbpub.com/
This free software is compliments of the authors.
mathews@fullerton.edu,  howell@westmont.edu

CHAPTER 10  APPLICATIONS of HARMONIC FUNCTIONS
Section 10.1  Preliminaries
Section 10.2  Invariance of Laplace's Equation and the Dirichlet Problem
Section 10.3  Poisson's Integral for the Upper Half-Plane
Section 10.4  Two-Dimensional Mathematical Models
Section 10.5  Steady State Temperatures
Section 10.6  Two Dimensional Electrostatics
Section 10.7  Two-Dimensional Fluid Flow
Section 10.8  The Joukowski Airfoil
Section 10.9  Schwarz-Christoffel Transformation
Section 10.10  Images of a Fluid Flow
Section 10.11  Sources and Sinks
GoTo Chapter 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11.
GoTo Contents

CHAPTER 10  APPLICATIONS of HARMONIC FUNCTIONS

\(SetOptions[Graphics, DefaultColor -> RGBColor[0, 0, 1]]; 

    Get["Graphics`ComplexMap`"]; 

    Clear[a, a0, a1, a2, a3, A, b, B, CR, Δx, Δy, 
      eq1, eq2, eqns, f, f1, F, F1, F2, F3, F4, g, G, lines, n, P, 
      Ψ, ϕ, Φ, r, r0, solset, θ, t, T, 
      T1, T2, Q, u, u0, U, U1, v, v0, V, V1, w, w0, W, W1, x, x1, x2, x3, y, 
      z, z1, z2, z3, z4]; \)

Section 10.1  Preliminaries

In most applications involving harmonic functions it is required to find a harmonic function that takes on prescribed values along certain contours.  We will assume the reader is familiar with the material in Sections 2.5, 3.3, 5.1, and 5.2.

Example 10.1,  Page 310.  A specific example of the general case. Find a function
u(x,y) that is harmonic in the vertical strip 1 < x < 2  and takes on the boundary values
u(1,y) = 6   for all y,
u(2,y) = 1   for all y.

\(Clear[a, b, U, U1, U2, x, y]; 
a\  = \ 1; \ 
b\  = \ 2; \ 

    U1\  = \ 6; \ 
U2\  = \ 1; 

    U[x_, y_]\  = \ U1\  + \ \(\(U2\  - \ U1\)\/\(b - a\)\) \((x - a)\); 

    ContourPlot[U[x, y], {x, 1, 2}, {y, \(-1\), 1}, 
      Contours\  -> \ {1, 2, 3, 4, 5, 6}, PlotPoints -> 10, 
      AspectRatio -> 2, ColorFunction -> Hue]; 

    Print["Contour plot for U[x,y] = ", U[x, y]]; \)

Example 10.2,  Page 310.  A specific example of the general case. Find a function
u(x,y) that is harmonic in the sector  0 < Arg(z) < ((2 π)/3)  and takes on the boundary values
u(x,0) = 1   for x > 0,
u(x,y) = 9   for all points on the ray  r > 0,  θ = α.

\(Clear[a, b, C1, C2, U, x, y]; 
a\  = \ \(2  π\)\/3; \ 

    C1\  = \ 1; \ 
C2\  = \ 9; 

    U[x_, y_]\  = \ C1\  + \ \(C2\  - \ C1\)\/a\ ArcTan[x, y]; 

    ContourPlot[U[x, y], {x, \(-2\), 2}, {y, 0.0000001, 2.0000001}, 
      Contours\  -> \ {2, 3, 4, 5, 6, 7, 8, 9}, PlotPoints -> 35, 
      AspectRatio -> 1\/2, ColorFunction -> Hue]; 

    Print["Contour plot for U[x,y] = ", U[x, y]]; \)

Example 10.3,  Page 311.  A specific example of the general case. Find a function
u(x,y) that is harmonic in the annulus   1 < |z| < R   and takes on the boundary values
u(x,y) = 1   when  |z| = 1,
u(x,y) = 6   when  |z| =4.

\(Clear[x, y, R, U, K1, K2]; 
R\  = \ 4; \ 
K1\  = \ 1; \ 
K2\  = \ 6; 
    
U[x_, \ y_]\  = \ 
      K1\  + \ \((K2\  - \ K1)\)*Log[\@\(x\^2\  + \ y\^2\)]\/Log[R]; 

    ContourPlot[U[x, y], {x, \(-4.2\), 4.2}, {y, \(-4.2\), 4.2}, 
      Contours\  -> \ {1, 2, 3, 4, 5, 6}, PlotPoints -> 25, 
      AspectRatio -> 1, ColorFunction -> Hue]; 

    Print["Contour plot for U[x,y] = ", U[x, y]]; \)

GoTo Chapter 10

Section 10.2  Invariance of Laplace's Equation and the Dirichlet Problem

Let  D  be a domain whose boundary is made up of piecewise smooth contours joined end to end.  The Dirichelet problem is to find a function  ϕ  that is harmonic in  D  such that  ϕ  takes on prescribed values at points on the boundary.  Let us first study the problem in the upper half plane.

Theorem 10.1,  Page 312.  (Invariance of Laplace's Equation)  Let  Φ (u, v) be harmonic in a domain  G  in the  w plane.  Then  Φ (u, v)  satisfies Laplace's equation  (∂_(u, u))Φ (u, v) + 
          (∂_(v, v))Φ (u, v) = 0  at each point  w = u + i v  in  G.   If  w = f (z) = u (x, y) + i v (x, y)  is a conformal mapping from a domain  D  in the z plane onto  G, then the composition  ϕ (x, y) =   
        Φ (u (x, y), u (x, y))  is harmonic in  D,  and  ϕ (x, y)  satisfies Laplace's equation  (∂_(x, x))ϕ (u, y) + 
          (∂_(y, y))ϕ (x, y) = 0  at each point  z = x + i y  in  D.

Theorem 10.2,  Page 314.  (N-Value Dirichlet Problem for the Upper Half Plane)  
Let  (Subscript[u, 1]) <  (Subscript[u, 2]) < ⋯ <  (Subscript[u, N - 1])  denote  N - 1  real constants.   The function
Φ (u, v)=(Subscript[a, N - 1])+(1/π)(∑_(k = 1)^(N - 1))((Subscript[a, k - 1]) - (Subscript[a, k])) Arg (w - (Subscript[u, k]))
		            =(Subscript[a, N - 1])+(1/π)(∑_(k = 1)^(N - 1))((Subscript[a, k - 1]) - (Subscript[a, k]))Arctan((v/(u - (Subscript[u, k]))))
is harmonic in the upper half plane  Im (w) > 0  and takes on the boundary values
\(Φ \((u, 0)\) = a\_0\ \ \ \ \ \ for\ \ \ u < u\_1, 

      Φ \((u, 0)\) = 
        \(a\_k\ \ \ \ \ \ for\ \ \ u\_k < u < \ u\_\(k + 1\)\ \ \ for\ \ k = 1
          \), 2, ⋯, N - 2, 

      Φ \((u, 0)\) = 
        a\_\(N - 1\)\ \ \ for\ \ \ u > \(u\_\(N - 1\).\)\)

Load Mathematica's  Graphics3D package.  Make sure this is done only ONCE during a Mathematica  session.

Get["Graphics`ParametricPlot3D`"]

Example 10.4,  Page 313.  Show that  ϕ(x,y)  =  arctan(((2 x)/((x^2)  +  (y^2)  -  1)))  
is harmonic in the disk  |z| < 1.

\(Clear[x, y, ϕ]; 

    ϕ[x_, y_]\  = \ ArcTan[\(2  x\)\/\(x\^2\  + \ y\^2\  - \ 1\)]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; \)

Verify that Laplace's equation holds.

\(ϕxx\  = \ Together[∂\_\(x, x\)ϕ[x, y]]; 

    ϕyy\  = \ Together[∂\_\(y, y\)ϕ[x, y]]; 

    Print[\*""\!\(ϕ\_xx\)(x,y) = "", \ ϕxx]; 

    Print[\*""\!\(ϕ\_yy\)(x,y) = "", \ ϕyy]; 

    Print[\*""\!\(ϕ\_xx\)(x,y) + \!\(ϕ\_yy\)(x,y) = "", 
      ϕxx\  + \ ϕyy]; \)

Show that  ϕ(x,y)  is the imaginary part of a certain analytic function.

\(Clear[f, g, u, U, U1, v, V, V1, x, y, z, w, W, W1]; 
x /: Im[x] = 0; 

    y /: Im[y] = 0; 

    g[z_]\  = \ \(I\  + \ z\)\/\(I\  - \ z\); 

    U\  = \ Together[
        ComplexExpand[Re[g[x + I\ y]], 
\ \ 
          TargetFunctions -> {Re, Im}]]; 

    V\  = \ Together[
        ComplexExpand[Im[g[x + I\ y]], 
\ \ 
          TargetFunctions -> {Re, Im}]]\ ; 

    W1\  = \ Log[U\  + \ I\ V]; 

    W2\  = \ ComplexExpand[W1, TargetFunctions -> {Re, Im}]; 

    W3\  = \ Expand[MapAll[Together, W2]]; 

    U1\  = \ MapAll[Together, 
        ComplexExpand[Re[W1], TargetFunctions -> {Re, Im}]]; 

    V1\  = \ ComplexExpand[Im[W1], TargetFunctions -> {Re, Im}]; 

    Q\  = \ ArcTan[
        Simplify[
          V1[[2]]\/V1
              [[1]]]]; 

    Print["   g[z] = ", g[z]]; 

    Print["g[x+Iy] = ", g[x + I\ y]]; 

    Print["g[x+Iy] = ", 
      ComplexExpand[g[x + I\ y], 
\ \ 
        TargetFunctions -> {Re, Im}]]; 

    Print["g[x+Iy] = ", U\  + \ I\ V]; 

    Print["  "]; 
Print["Log[g[x+Iy]] = ", W1]; 

    Print["Log[g[x+Iy]] = ", W3]; 
Print["  "]; 

    Print["Re[Log[g[x+Iy]]] = ", U1]; 

    Print["Im[Log[g[x+Iy]]] = ", V1]; 

    Print["Im[Log[g[x+Iy]]] = ", Q]; \)

Remark.  Thus,  ϕ(x,y)  =  arctan(((2 x)/((x^2)  +  (y^2)  -  1)))  =  Im(Log(((I  +  z)/(I  -  z)))).

Example 10.5,  Page 313.  Show that the function
Φ(u,v)  =  (1/π)ArcTan[(v/(u  -  (Subscript[u, 0])))]   =   (1/π)Arg[w - (Subscript[u, 0])]   is
harmonic in the upper half plane  Im(w) > 0  and has boundary values:
Φ(u,0)  =  0    for   u > (Subscript[u, 0]),
Φ(u,0)  =  1    for   u < (Subscript[u, 0]).

\(Clear[g, G, P, Ψ, Φ, Q, u, u0, U, v, V, w, w0]; 
    
g[w_, u0_]\ \  = \ 1\/π\ Log[w\  - \ u0]; 

    G[w_, u0_]\  = \ ComplexExpand[g[w, u0]]; 

    P[w_, u0_]\  = \ ComplexExpand[Re[g[w, u0]]]; 

    Q[w_, u0_]\  = \ ComplexExpand[Im[g[w, u0]]]; 

    Ψ[u_, v_, u0_]\  = \ 
      ReplaceAll[P[w, u0], w -> u + I\ v]; 

    Φ[u_, v_, u0_]\  = \ 
      ReplaceAll[Q[w, u0], w -> u + I\ v]; 

    Print["  g[w,u0] = ", g[w, u0]]; 

    Print["  g[w,u0] = ", G[w, u0]]; 

    Print["Ψ[u,v,u0] = ", P[w, u0], " = ", 
      Ψ[u, v, u0]]; 

    Print["Φ[u,v,u0] = ", P[w, u0], " = ", 
      Φ[u, v, u0]]; \)

For illustration, set  (Subscript[u, 0])  =1.

\(u = 2; \ v = 0; \ u0 = 1; 
Print["u = ", u, "    u0 = ", u0]; 
    
Print[u, " > ", u0, " is ", u > u0, 
      ",  Φ[u,v,u0] = ", Φ[u, v, u0]]; 
u = 0; \ 
    v = 0; \ u0 = 1; 
Print["u = ", u, "    u0 = ", u0]; 

    Print[u, " < ", u0, " is ", u < u0, 
      ",  Φ[u,v,u0] = ", Φ[u, v, u0]]; \)

\(Clear[u, v]; 

    Plot3D[Φ[u, v, 1], {u, 0, 2}, {v, 0, 1}, 
\ \ \ \ \ \ \ 
      AxesLabel -> {"u", "v", " "}, 
      ColorFunction -> Hue]; 

    Print["      Φ[u,v,1] = ", Φ[u, v, 1]]; \)

Example 10.6,  Page 315.  Find a function  ϕ(x,y) that is harmonic in the 
upper half-plane  Re z > 0,  which takes on the boundary values:
ϕ(x,0) = 4,   when   x < -1,
ϕ(x,0) = 1,   when   -1 < x < 0,
ϕ(x,0) = 3,   when   0 < x < 1,
ϕ(x,0) = 2,   when   1 < x.

\(Clear[a0, a1, a2, a3, ϕ, x, x1, x2, x3, y, z]; 
x1\  = \ \(-1\); 
    
x2\  = \ \ 0; 
x3\  = \ \ 1; 
a0\  = \ \ 4; 
a1\  = \ \ 1; 

    a2\  = \ \ 3; 
a3\  = \ \ 2; 
z\  = \ x\  + \ I\ y; 

    ϕ[x_, y_] = 
      a3 + \(1\/π\) \((a0 - a1)\) Arg[z - x1] + 
        \(1\/π\) \((a1 - a2)\) Arg[z - x2]\  + \ 
        \(1\/π\) \((a2 - a3)\) Arg[z - x3]; 

    Print["ϕ[x,y] = "]; \ Print[ϕ[x, y]]; \)

\(ContourPlot[ϕ[x, y], {x, \(-2.0\), 2.0}, {y, 0.0001, 1.5001}, 
      Contours -> 15, PlotPoints -> 35, AspectRatio -> 1.5\/4, 
      ColorFunction -> Hue]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; \)

Example 10.7,  Page 316.  Find a function  ϕ(x,y) that is harmonic in the 
upper half-plane  Re z > 0,  which takes on the boundary values:
ϕ(x,0) = 1    when    |x| < 1,
ϕ(x,0) = 0    when    |x| < 1.

\(Clear[a0, a1, a2, ϕ, x, x1, x2, y, v, z]; 
x1\  = \ \(-1\); 

    x2\  = \ \ 1; 
a0\  = \ \ 0; 
a1\  = \ \ 1; 
a2\  = \ \ 0; 

    z\  = \ x\  + \ I\ y; 

    ϕ[x_, y_] = 
      a2 + \(1\/π\) \((a0 - a1)\) Arg[z - x1] + 
        \(1\/π\) \((a1 - a2)\) Arg[z - x2]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; \)

\(Plot3D[ϕ[x, y], {x, \(-2\), 2}, {y, 0.00001, 1}, 
      AxesLabel -> {"x", "y", "v"}, 
      ColorFunction -> Hue]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; \)

Example 10.8,  Page 317.  Find a function  ϕ(x,y) that is harmonic in 
the unit disk  |z| < 1,  which takes on the boundary values:
ϕ(x,y) = 0    when   z = exp(iθ)  and  0 < θ < π,
ϕ(x,y) = 1    when   z = exp(iθ)  and  π < θ < 2π.

\(Clear[f, ϕ, Q, u, U, V, w, W, x, y, z]; 

    f[z_]\  = \ \(I \((1\  - \ z)\)\)\/\(1\  + \ z\); 

    W\  = \ ComplexExpand[f[x + \ I\ y], \ 
        TargetFunctions -> {Re, Im}]; 

    U\  = \ ComplexExpand[Re[W], \ TargetFunctions -> {Re, Im}]; 

    V\  = \ ComplexExpand[Im[W], \ TargetFunctions -> {Re, Im}]; 

    V\  = \ Together[V]; 
Q\  = \ 1\/π\ ArcTan[V\/U]; 

    ϕ[x_, y_]\  = \ 1\/π\ ArcTan[Cancel[V\/U]]; 

    Print["     f[z] = ", f[z]]; 

    Print["  f[x+Iy] = ", U\  + \ I\ V]; 

    Print[" "]; 
Print["Log[f[z]] = ", Log[f[z]]]; 

    Print["Log[f[z]] = ", Log[U\  + \ I\ V]]; 

    Print[" "]; 
Print["   ϕ[x,y] = ", Arg[f[z]]]; 

    Print["   ϕ[x,y] = ", Arg[U\  + \ I\ V]]; 

    Print[" "]; 
Print["   ϕ[x,y] = ", Q]; 

    Print["   ϕ[x,y] = ", ϕ[x, y]]; \)

Example 10.9,  Page 318.  Find a function  ϕ(x,y) that is harmonic in the 
upper half disk  H: y > 0,  |z| < 1,  which takes on the boundary values:
ϕ(x,y) = 0    when   z = exp(iθ)  and  0 < θ < π,
ϕ(x,0) = 1    for   -1 < x < 1.

\(Clear[F, U, V, w, W, ϕ, x, y, z]; 

    F[z_]\  = \ \(I \((1\  - \ z)\)\)\/\(1\  + \ z\); 

    W\  = \ ComplexExpand[F[x + \ I\ y], \ 
        TargetFunctions -> {Re, Im}]; 

    U\  = \ Together[ComplexExpand[Re[W]]]; 

    V\  = \ Together[ComplexExpand[Im[W]]]; 

    ϕ[x_, y_]\  = \ 2\/π\ ArcTan[V\/U]; 

    Print["w = ", F[z], " = ", U + I\ V]; 

    Print[" "]; 

    Print[\*""\!\(2\/π\)Log[w] = "", 
      \(2\/π\) Log[U + I\ V]]; 
Print[" "]; 

    Print["   ϕ[x,y] = ", ϕ[x, y]]; \)

\(Clear[f, r, θ]; 

    f[r_, θ_]\  = \ 
      Simplify[ReplaceAll[
          ϕ[x, y], {x -> r\ Cos[θ], 
            y -> r\ Sin[θ]}]]; 

    CylindricalPlot3D[f[r, θ], {r, 0.001, 1}, {θ, 0.001, 3.13}, 
      ViewPoint -> {1, 1, 1}, 
      AxesLabel -> {"x", "y", "u"}]; 

    Print["f[r,θ] = ", f[r, θ]]; \)

Example 10.10,  Page 320.  Find a function  ϕ(x,y) that is harmonic in the 
upper half disk  G:  x > 0, y > 0,  |z| < 1,  which takes on the boundary values:
ϕ(x,y) = 0    when   z = exp(iθ)  and  0 < θ < (π/2),
ϕ(x,0) = 1    for   0 ≤ x < 1,
ϕ(0,y) = 1    for   0 ≤ y < 1.

\(Clear[F, U, V, w, W, ϕ, x, y, z]; 

    F[z_]\  = \ \(I \((1\  - \ z\^2)\)\)\/\(1\  + \ z\^2\); 

    W\  = \ ComplexExpand[F[x + \ I\ y], \ 
        TargetFunctions -> {Re, Im}]; 

    U\  = \ Together[ComplexExpand[Re[W]]]; 

    V\  = \ Together[ComplexExpand[Im[W]]]; 

    ϕ[x_, y_]\  = \ 2\/π\ ArcTan[V\/U]; 

    Print["w = ", F[z], " = ", U + I\ V]; \ 

    Print[" "]\)\(Print[\*""\!\(2\/π\)Log[w] = "", 
      \(2\/π\) Log[U + I\ V]]; 
Print[" "]; 

    Print["   ϕ[x,y] = ", ϕ[x, y]]; \)

\(Clear[f, r, θ]; 

    f[r_, θ_]\  = \ 
      Simplify[ReplaceAll[
          ϕ[x, y], {x -> r\ Cos[θ], 
            y -> r\ Sin[θ]}]]; 

    CylindricalPlot3D[f[r, θ], {r, 0.001, 1}, {θ, 0.001, 1.56}, 
      ViewPoint -> {2, \(-1\), 1}, 
      AxesLabel -> {"x", "y", "u"}]; 

    Print["f[r,θ] = ", f[r, θ]]\)

GoTo Chapter 10

Section 10.3  Poisson's Integral for the Upper Half-Plane

Theorem 10.3,  Page 323.  (Poisson's Integral Formula)  Let  U(t)  be a real-valued 
function that is piecewise continuous and bounded for all real t.  The function
ϕ(x,y)  =  (y/π)(∫_(-∞)^∞)((U (t))/((((x - t))^2)  +  (y^2)))dt
is  harmonic in the upper half plane  Im(z) > 0 and has the boundary values
ϕ(x,0)  =  U(t)  wherever  U  is continuous.

Example 10.11,  Page 324.  Find the function  ϕ(x,y)  that is harmonic in 
the upper half-plane  Im(z) > 0  and has the boundary values.
     ϕ(x,0)  =  1     for    -1 < x < 1,
     ϕ(x,0)  =  0     for    |x| > 1.

\(Clear[t, ϕ, U, v, x, y, z]; 
U[t_]\  = \ 1; 

    ϕ[x_, y_]\  = \ 
      Apart[Together[
          y\/π\ 
            \(∫\_\(-1\)\%1\( U[t]\/\(\((x - t)\)\^2 + y\^2\)\) 
                \[DifferentialD]t\)]]; 

    v[x_, y_]\  = \ 
      ReplaceAll[ϕ[x, y], ArcTan[z_] -> \(-ArcTan[1\/z]\)]; 

    Print["  U[t]  =  ", U[t]]; \ \ 
Print[" "]; 

    Print[\*""ϕ[x,y]  =  (y/π) (∫_(-1)^1)U[t]/((((x - t))^2) + (y^2))\[DifferentialD]t""]; \ 
Print[" "]; 

    Print["ϕ[x,y]  =  ", ϕ[x, y]]; 
Print[" "]; 

    Print["ϕ[x,y]  =  ", v[x, y]]; \)

Using the trigonometric identity  arctan(-t)  =  - arctan(t), the above result can be written as
 ϕ(x,y)  =  (1/π) arctan((y/(x  -  1)))  -  (1/π) arctan((y/(x  +  1))).

\(Print["     ϕ[x,y] = ", ϕ[x, y]]; 

    Print[\*""\!\(lim\+\(y -> 0\)\) ϕ[-2,y] = "", 
      Limit[ϕ[\(-2\), y], \ y -> 0]]; 

    Print[\*"" \!\(lim\+\(y -> 0\)\) ϕ[0,y] = "", 
      Limit[ϕ[0, y], \ y -> 0]]; 

    Print[\*"" \!\(lim\+\(y -> 0\)\) ϕ[2,y] = "", 
      Limit[ϕ[2, y], \ y -> 0]]; \)

Example 10.12,  Page 324.  Find the funciton  ϕ(x,y)  that is harmonic in 
the upper half-plane  Im(z) > 0  and has the boundary values.
     ϕ(x,0)  =  x     for    -1 < x < 1,
     ϕ(x,0)  =  0     for    |x| > 1.

\(\(Clear[t, ϕ, U, v, x, y, z]; 
U[t_]\  = \ t; 

    ϕ[x_, y_]\  = \ 
      Expand[y\/π\ 
          \(∫\_\(-1\)\%1\( U[t]\/\(\((x - t)\)\^2 + y\^2\)\) 
              \[DifferentialD]t\)]; 

    v[x_, y_]\  = 
      ReplaceAll[ϕ[x, y], ArcTan[z_] -> \(-ArcTan[1\/z]\)]; 

    Print["  U[t]  =  ", U[t]]; \ \ 
Print[" "]; 

    Print[\*""ϕ[x,y]  =  (y/π) (∫_(-1)^1)U[t]/((((x - t))^2) + (y^2))\[DifferentialD]t""]; \ 
Print[" "]; 

    Print["ϕ[x,y]  =  ", ϕ[x, y]]; 
Print[" "]; 

    Print["ϕ[x,y]  =  ", v[x, y]]; 
\)\)

Using the identities   arctan(-t)  =  - arctan(t),   and    ln((a/b))  =  ln(a) - ln(b),  the
above result can be written as: ϕ(x,y)  =  (1/π) arctan((y/(x  -  1)))  -  (1/π) arctan((y/(x  +  1))).
 ϕ(x,y)  =  (x/π)arctan((y/(x  -  1))) - (x/π)arctan((y/(x  +  1))) + (x/(2 π))ln((((((x - 1))^2)  +   (y^2))/((((x + 1))^2)  +   (y^2)))).

\(Print["ϕ[x,y] = ", ϕ[x, y]]; 

    Print[\*"" (lim_(y -> 0)) ϕ[-2,y] = "", 
      Limit[ϕ[\(-2\), y], \ y -> 0]]; 

    Print[\*""(lim_(y -> 0)) ϕ[((-1)/2),y] = "", Limit[ϕ[\(-1\)\/2, y], \ y -> 0]]; 

    Print[\*""  (lim_(y -> 0)) ϕ[0,y] = "", 
      Limit[ϕ[0, y], \ y -> 0]]; 

    Print[\*"" (lim_(y -> 0)) ϕ[(1/2),y] = "", 
      Limit[ϕ[1\/2, y], \ y -> 0]]; 

    Print[\*""  (lim_(y -> 0)) ϕ[2,y] = "", 
      Limit[ϕ[2, y], \ y -> 0]]; \)

\(Plot3D[ϕ[x, y], {x, \(-1.5\), 1.5}, {y, 10\^\(-8\), 0.5}, 
      PlotPoints -> {18, 11}, 
      PlotRange -> {{\(-1.5\), 1.5}, {0, 0.5}, {\(-1.05\), 1.05}}, 
      AxesLabel -> {"x", "y", "u  "}, 
      Ticks -> {Range[\(-5\), 5, 1], Range[\(-5\), 5, 0.5], 
          Range[\(-1\), 1, 0.5]}, BoxRatios -> {1, 1, 0.8}, 
      ViewPoint -> {2.5, \(-3\), 2}, ColorFunction -> Hue]; 

    Print["ϕ[x,y]  =  ", ϕ[x, y]]; \)

Example 10.13,  Page 325.  Use Poisson's Integral formula to find the 
harmonic funciton  ϕ(x,y)  in the upper half-plane that takes on the boundary values.
     ϕ(t,0)  =  -1     for   t < -1,
     ϕ(t,0)  =   t      for   -1 < t < 1,
     ϕ(t,0)  =   1     for   1 < t.

Solution.  Using techniques from Section 10.2, the function
     v(x,y)  =  1 - (1/π) arctan((y/(x  +  1))) - (1/π) arctan((y/(x  -  1)))
is harmonic in the upper half plane and takes on the boundary values
     v(x,0)  =  -1     for     x < -1,
     v(x,0)  =   0     for    -1 < x < 1,
     v(x,0)  =   1     for     1 < x.
Thus, we should add it to the solution  ϕ(x,y)  in Example 10.12
to obtain the desired result.  However, with Mathematica  we need
to use   ArcTan[x+1,y]   instead of   ArcTan[(y/(x + 1))].

\(Clear[t, ϕ, U, v, V, x, y, z]; 
U[t_]\  = \ t; 

    ϕ[x_, y_]\  = \ 
      Expand[y\/π\ 
          \(∫\_\(-1\)\%1\( U[t]\/\(\((x - t)\)\^2 + y\^2\)\) 
              \[DifferentialD]t\)]; 

    v[x_, y_]\  = \ 
      ReplaceAll[ϕ[x, y], ArcTan[z_] -> \(-ArcTan[1\/z]\)]; 

    V[x_, y_]\  = \ 
      1\  - \ \(1\/π\) ArcTan[x + 1, y]\  - 
        \(1\/π\) ArcTan[x - 1, y]; 

    ϕ[x_, y_]\  = \ ϕ[x, y]\  + \ V[x, y]; 

    v[x_, y_]\  = \ v[x, y]\  + \ V[x, y]; 
Print["  U[t] = ", U[t]]; 
    \ \ 
Print[" "]; 
Print["V[x,y] = ", V[x, y]]; 

    Print[" "]; 

    Print[\*""ϕ[x,y] = (y/π)(∫_(-1)^1)U[t]/((((x - t))^2) + (y^2))\[DifferentialD]t + V[x,y]""]; \ 
Print[" "]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; \ 
Print[" "]; 

    Print["ϕ[x,y] = ", v[x, y]]; \)

\(Plot3D[ϕ[x, y], {x, \(-2\), 2}, {y, 10\^\(-8\), 0.5}, 
      PlotPoints -> {24, 11}, 
      PlotRange -> {{\(-2\), 2}, {0, 0.5}, {\(-1.05\), 1.05}}, 
      AxesLabel -> {"x", "y", "u  "}, 
      Ticks -> {Range[\(-5\), 5, 1], Range[\(-5\), 5, 0.5], 
          Range[\(-1\), 1, 0.5]}, BoxRatios -> {1, 1, 0.8}, 
      ViewPoint -> {\(-2.5\), \(-3\), 2}, ColorFunction -> Hue]; 

    Print["ϕ[x,y]  =  ", ϕ[x, y]]; \)

GoTo Chapter 10

Section 10.4  Two-Dimensional Mathematical Models

We now turn our attention to problems involving steady state heat flow, electrostatics, and ideal fluid flow that can be solved by conformal mapping techniques. The method uses conformal mapping to transform a region in which the problem is posed to one in which the solution is easy to obtain.

Theorem 10.4,  Page 328.  (Orthogonal Families of Level Curves)  
Let  f (z)  =  u (z, y)  +  i v (x, y)  be a harmonic function.  Then the following two families of curves are orthogonal:  {u (z, y) = (Subscript[K, 1]):  (Subscript[K, 1]) is a real constant}  and  {v (z, y) = (Subscript[K, 2]):  (Subscript[K, 2]) is a real constant}.
Proof.  Use the Cauchy Riemann equations  (∂_x) u (x, y) = (∂_y) v (x, y)  and (∂_y) u (x, y) = -(∂_x) v (x, y),  which are implemented by using the Mathematica  substitutions:

\(Clear[CR, x, y, u, v]; 

    CR\  = \ {∂\_y\ v[x, y] -> ∂\_x\ u[x, y], 
        ∂\_y\ u[x, y] -> \ \(-∂\_x\ v[x, y]\)}\)

The normal vectors to the curves   u(x,y) = (Subscript[K, 1])   and    v(x,y) = (Subscript[K, 2])   are:

\(N1\  = \ {∂\_x\ u[x, y], ∂\_y\ u[x, y]}; 

    N2\  = \ {∂\_x\ v[x, y], ∂\_y\ v[x, y]}; 

    Print[\*""\!\(N\_1\) = "", N1]; 

    Print[\*""\!\(N\_2\) = "", N1]; \)

We will show that  (Subscript[N, 1])  is orthogonal to  (Subscript[N, 2])  by forming the dot product:

\(Print[\*""\!\(N\_1\).\!\(N\_1\) = "", N1, ".", N2]; 

    Print[\*""\!\(N\_1\).\!\(N\_1\) = "", N1.N2]; \)

Using the Cauchy Riemann equations, this is simplified to get:

\(\(Print[\*""\!\(N\_1\).\!\(N\_1\) = "", ReplaceAll[N1.N2, CR]]; 
    \)\)

Therefore  (Subscript[N, 1])  is orthogonal to  (Subscript[N, 2]) .

GoTo Chapter 10

Section 10.5  Steady State Temperatures

In the theory of heat conduction the assumption is made that heat flows in the direction of decreasing temperature.  We also assume that the time rate at which heat flows across a surface area is proportional to the component of the temperature gradient in the direction perpendicular to the surface area.    If the temperature  T (x, y)  does not depend on time, then  T (x, y)  can be shown to be a harmonic function.

Example 10.14,  Page 331.  Find the temperature  T(x,y)  at each point in the 
horizontal strip  a < Im(z) < b  if the temperature at points on the boundary satisfies
     T(x,a)  =  (Subscript[T, 1])     
     T(x,b)  =  (Subscript[T, 2])

\(Clear[a, b, lines, T, T1, T2, x, y]; 

    T[x_, y_]\  = \ T1\  + \ \(T2\  - \ T1\)\/\(b\  - \ a\)\ \((y\  - \ a)\); 
    
Print["T[x,y] = ", T[x, y]]; \)

For illustration we can choose  a = 1,  b = 3,  (Subscript[T, 1]) = 50  and  (Subscript[T, 2]) = 10.

\(a\  = \ 1; 
b\  = \ 3; 
T1\  = \ 50; 
T2\  = \ 10; 

    Print["T[x,y] = ", T[x, y]]; 
Print["T[x,1] = ", T[x, 1]]; 

    Print["T[x,2] = ", T[x, 2]]; 
Print["T[x,3] = ", T[x, 3]]; \)

\(lines\  = \ {10, 20, 30, 40, 50}; 

    ContourPlot[T[x, y], {x, \(-3\), 3}, {y, 0, 3}, 
      Contours\  -> \ lines, PlotPoints -> 10, 
      AspectRatio -> 1\/2, ColorFunction -> Hue]; 

    Print["T[x,y] = ", T[x, y]]; 
Print["Contours = ", lines]; \)

Example 10.15,  Page 331.  Find the temperature  T(x,y)  at each point in the 
upper half plane  Im(z) > 0,  if the temperature at points on the boundary satisfies
     T(x,0)  =  (Subscript[T, 1])     for  x > 0,
     T(x,0)  =  (Subscript[T, 2])     for  x < 0.

Clear[a, b, T, T1, T2, x, y]\(T[x_, y_]\  = \ 
      T1\  + \ \(\(T2\  - \ T1\)\/π\) Arg[x\  + \ I\ y]; 

    Print["T[x,y] = ", T[x, y]]; \)

For illustration we can choose  a = 1,  b = 3,  (Subscript[T, 1]) = 10  and  (Subscript[T, 2]) = 50.

\(T1\  = \ 50; 
T2\  = \ 10; 
Print["       T[x,y] = ", T[x, y]]; 

    Print["       T[1,0] = ", T[1, 0]]; 

    Print[\*"" T[\!\(\@2\/2\),\!\(\@2\/2\)] = "", T[\@2\/2, \@2\/2]]; 

    Print["       T[0,1] = ", T[0, 1]]; 

    Print[\*""T[\!\(\(-\@2\)\/2\),\!\(\@2\/2\)] = "", 
      T[\(-\@2\)\/2, \@2\/2]]; 
Print["      T[-1,0] = ", T[\(-1\), 0]]; 
    \)

\(lines\  = \ Table[10 + 5  k, {k, 0, 8}]; 

    ContourPlot[T[x, y], {x, \(-3\), 3}, {y, 0, 3}, 
      Contours\  -> \ lines, PlotPoints -> 35, 
      AspectRatio -> 1\/2, ColorFunction -> Hue]; 

    Print["T[x,y] = ", T[x, y]]; 
Print["Contours = ", lines]; \)

Example 10.16,  Page 332.  Find the temperature  T(x,y)  at each point in the upper 
half-disk  H: Im z > 0,  |z| < 1  if the temperature at points on the boundary satisfies
     T(x,y)  =  100     for   (x,y)  on the upper semi-circle,
     T(x,0)  =  50       for   -1 < x < 1,  on the diameter.

\(Clear[f, lines, T, U, V, w, W, x, y]; 

    f[z_]\  = \(\ I \((1\  - \ z)\)\)\/\(1\  + \ z\); 

    W\  = \ ComplexExpand[f[x\  + \ I\ y]]; 

    U\  = \ Together[ComplexExpand[Re[W], TargetFunctions -> {Re, Im}]]; 
    
\(V\  = \ 
      Together[ComplexExpand[Im[W], TargetFunctions -> {Re, Im}]]; \); 

    T[x_, y_]\  = \ 100\  - \ 100\/π\ ArcTan[V\/U]; 

    Print["w = ", f[z], " = ", U\  + \ I\ V]; \ 

    Print[" "]; 

    Print[\*""T[x,y] = 100 - \!\(100\/π\) ArcTan[\!\(v\/u\)]""]; 

    Print["T[x,y] = ", T[x, y]]; 
Print[" "]; 

    Print[\*"" \!\(lim\+\(y -> 0\)\) T[\!\(\(-1\)\/2\),y] = "", 
      Limit[T[\(-1\)\/2, y], \ y -> 0]]; 

    Print[\*""   \!\(lim\+\(y -> 0\)\) T[0,y] = "", 
      Limit[T[0, y], \ y -> 0]]; 

    Print[\*""  \!\(lim\+\(y -> 0\)\) T[\!\(1\/2\),y] = "", 
      Limit[T[1\/2, y], \ y -> 0]]; 

    Print[\*"" T[\!\(\@2\/2\),\!\(\@2\/2\)] = "", T[\@2\/2, \@2\/2]]; 

    Print["       T[0,1] = ", T[0, 1]]; 

    Print[\*""T[\!\(\(-\@2\)\/2\),\!\(\@2\/2\)] = "", 
      T[\(-\@2\)\/2, \@2\/2]]; \)

\(lines\  = \ Table[60 + 10  k, {k, 0, 4}]; 

    ContourPlot[T[x, y], {x, \(-1\), 1}, {y, 0.01, 1}, 
      Contours\  -> \ lines, PlotPoints -> 35, 
      AspectRatio -> 1\/2, ColorFunction -> Hue]; 

    Print["T[x,y] = ", T[x, y]]; 
Print["Contours = ", lines]; \)

Example 10.17,  Page 334.  Find the temperature  T(x,y)  at each point in the 
upper half-plane where the temperature at points on the boundary satisfies
     T(x,0)  =  1      for  x > 1,
     T(x,0)  =  -1     for  x < -1,
     ((∂T)/(∂N))  =  (Subscript[T, y])(x,y)  =  0    for  -1 < x < 1.

\(Clear[lines, T, Ty, x, y]; 

    T[x_, y_]\  = \(\ 2\)\/π\ Re[ArcSin[x\  + \ I\ y]]; 

    Print["T[x,y] = ", T[x, y]]; 

    Print["T[x,y] = ", ComplexExpand[T[x, y]]]; \)

\(T[x_, y_]\  = \ 
      2\/π\ ArcSin[
          \(\@\(\((x + 1)\)\^2\  + \ y\^2\)\  - 
              \@\(\((x - 1)\)\^2\  + \ y\^2\)\)\/2]; 

    Ty[x_, y_]\  = \ ∂\_y\ T[x, y]; 

    Print["  T[x,y] = ", T[x, y]]; 

    Print[" T[-2,0] = ", T[\(-2\), 0]]; 

    Print[" T[-1,0] = ", T[\(-1\), 0]]; 

    Print[\*""T[\!\(\(-1\)\/2\),0] = "", T[\(-1\)\/2, 0]]; 

    Print["  T[0,0] = ", T[0, 0]]; 

    Print[\*"" T[\!\(1\/2\),0] = "", T[1\/2, 0]]; 

    Print["  T[1,0] = ", T[1, 0]]; 
Print["  T[2,0] = ", T[2, 0]]; 

    Print[" "]; 
Print[\*""\!\(T\_y\)[x,y] = "", Ty[x, y]]; 

    Print[\*""\!\(T\_y\)[\!\(\(-1\)\/2\),0] = "", Ty[\(-1\)\/2, 0]]; 

    Print[\*""  \!\(T\_y\)[0,0] = "", Ty[0, 0]]; 

    Print[\*"" \!\(T\_y\)[\!\(1\/2\),0] = "", Ty[1\/2, 0]]; \)

\(lines\  = \ Table[\(-1\)\  + \ 0.2  k, {k, 0, 10}]; 

    ContourPlot[T[x, y], {x, \(-4\), 4}, {y, 0.000001, 4}, 
      Contours\  -> \ 11, PlotPoints -> 25, 
      AspectRatio -> 1\/2, ColorFunction -> Hue]; 

    Print["T[x,y] = ", T[x, y]]; 
Print["Contours = "]; \ 

    Print[lines]; \)

GoTo Chapter 10

Section 10.6  Two Dimensional Electrostatics

A two-dimensional electrostatic field is produced by a system of charged wires, plates, and cylindrical conductors that are perpendicular to the z plane.  The wires, plates, and cylinders are assumed to be so long that the effects at the ends can be neglected.  This sets up an electric field  E (x, y),  the real and imaginary parts of  E (x, y)  are harmonic functions.

Example 10.18,  Page 343. Consider two parallel conducting planes that pass perpendicular to the z plane through the lines  x = a  and  y = b,  which are kept at the potentials  (Subscript[U, 1]) and (Subscript[U, 2]), respectively.   Then according to the result of Example 10.1, the electrical potential is
          ϕ(x,y)   = (Subscript[U, 1]) + (((Subscript[U, 2])  -  (Subscript[U, 1]))/(b  -  a))(x - a).

Example 10.19,  Page 343.  Find the electrical potential  ϕ(x,y)  in the region between 
the two infinite coaxial cylinders r = a  and  r = b, which are kept at potentials (Subscript[U, 1]) and (Subscript[U, 2]),
respectively.   For an illustration we will use the boundary values:
 ϕ(x,y) = 0      when  |z| = 1,
 ϕ(x,y) = 40    when  |z| = 2.

\(Clear[a, b, lines, ϕ, U1, U2]; 
a\ \  = \ 1; 
b\ \  = \ 2; 

    U1\  = \ 0; 
U2\  = \ 40; 

    ϕ[x_, y_]\  = \ 
      U1\  + \ \(\(U2\  - \ U1\)\/\(Log[b] - Log[a]\)\) 
          \((Log[\@\(x\^2 + y\^2\)]\  - \ Log[a])\); 

    Print["ϕ[x,y] = ", ϕ[x, y]]; 

    Print["ϕ[1,0] = ", ϕ[1, 0]]; 

    Print["ϕ[0,1] = ", ϕ[0, 1]]; 

    Print["ϕ[2,0] = ", ϕ[2, 0]]; 

    Print["ϕ[0,2] = ", ϕ[0, 2]]; \)

\(\(lines\  = \ Table[0 + 10  k, {k, 0, 4}]; 

    ContourPlot[ϕ[x, y], {x, \(-2\), 2}, {y, \(-2\), 2}, 
      Contours\  -> \ lines, PlotPoints -> 25, 
      AspectRatio -> 1, ColorFunction -> Hue]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; 

    Print["Contours = ", lines]; \ \)\)

Example 10.20,  Page 343.   Find the electrical potential  ϕ(x,y)  produced by 
two charged half-planes that are perpendicular to the z plane and pass through the 
rays  x < -1, y = 0  and  x > 1, y = 0  where the planes are kept at the fixed potentials:
 ϕ(x,0) = -300,  for  x < -1
 ϕ(x,0) = 300,   for  x > 1.

\(Clear[solset, x, y, ϕ, u]; 

    u[x_, y_]\  = \ 600\/π\ Re[ArcSin[x\  + \ I\ y]]; 

    Print["u[x,y] = ", u[x, y]]; 

    Print["u[x,y] = ", ComplexExpand[u[x, y]]]; \)

\(ϕ[x_, y_]\  = \ 
      600\/π\ 
        ArcSin[\(\@\(\((x + 1)\)\^2\  + \ y\^2\)\  - 
              \@\(\((x - 1)\)\^2\  + \ y\^2\)\)\/2]; 

    Print[" ϕ[x,y] = ", ϕ[x, y]]; 

    Print["ϕ[-1,0] = ", ϕ[\(-1\), 0]]; 

    Print[" ϕ[0,1] = ", ϕ[0, 1]]; 

    Print[" ϕ[1,0] = ", ϕ[1, 0]]; \)

\(\(lines\  = \ Table[\(-300\)\  + \ 60  k, {k, 0, 10}]; 

    ContourPlot[u[x, y], {x, \(-4\), 4}, {y, \(-4\), 4}, 
      Contours\  -> \ lines, PlotPoints -> 25, 
      AspectRatio -> 1, ColorFunction -> Hue]; 

    Print["u[x,y] = ", u[x, y]]; 
Print["Contours = "]; \ 

    Print[lines]; \ \)\)

The following commands will graph the contours as parametric curves.

\(Clear[F, X, Y, z, Z, w]; 
F[z_]\  = \ 600\/π\ ArcSin[z]; 

    solset\  = \ Solve[w\  == \ F[z], z]; 

    Z[u_, v_]\  = \ 
      ReplaceAll[solset[[1, 1, 2]], 
        w -> u\  + \ I\ v]; 

    X[u_, v_]\  = \ ComplexExpand[Re[Z[u, v]]]; 

    Y[u_, v_]\  = \ ComplexExpand[Im[Z[u, v]]]; 
Print["w = ", F[z]]; 

    Print["Z = ", Z[u, v]]; 
Print["Z = ", ComplexExpand[Z[u, v]]]; 
    
Print["X = ", X[u, v]]; 
Print["Y = ", Y[u, v]]; \)

\(\(eqns\  = \ N[Table[{X[k\ 60, v], Y[k\ 60, v]}, {k, \(-5\), 5}]]; 

    ParametricPlot[Evaluate[eqns], {v, \(-500\), 500}, 
      PlotRange -> {{\(-4\), 4}, {\(-4\), 4}}, AspectRatio -> 1]; 

    Print["u[x,y] = ", u[x, y]]; 
Print["Contours = "]; \ 

    Print[lines]; \ \)\)

Example 10.21,  Page 343.  Find the electrical potential  ϕ(x,y)  in the disk  
D: |z| < 1   that satisfies the boundary values:
 ϕ(x,y) = 80,   for  z = eiθ  and     0 < θ < (π/2),
 ϕ(x,y) =   0,   for  z = eiθ  and  (π/2) < θ < 2π.

\(Clear[f, F, ϕ, u, U, v, V, w, W, x, X, y, Y, z, Z]; 

    f[z_]\  = \ 
      \(\((1\  - \ I)\) \((z\  - \ I)\)\)\/\(z\  - \ 1
          \); 
W\  = \ \ ComplexExpand[f[x\  + \ I\ y]]; 

    U\  = \ Together[ComplexExpand[Re[W], TargetFunctions -> {Re, Im}]]; 
    
\(V\  = \ 
      Together[ComplexExpand[Im[W], TargetFunctions -> {Re, Im}]]; \); 

    ϕ[x_, y_]\  = \ 80\/π\ ArcTan[V\/U]; 
Print["w = ", f[z]]; 
    
Print["w = u + I v = ", U\  + \ I\ V]; \ 

    Print[" "]; 

    Print[\*""ϕ[x,y] = \!\(80\/π\) ArcTan[\!\(v\/u\)]""]; 

    Print["ϕ[x,y] = ", ϕ[x, y]]; \)

\(Clear[F, u, v, w, X, Y, z, Z]; 

    F[z_]\  = \ 
      \(\((1\  - \ I)\) \((z\  - \ I)\)\)\/\(z\  - \ 1
          \); 
solset\  = \ Solve[w\  == \ F[z], z]; 

    Z[u_, v_]\  = \ 
      ReplaceAll[solset[[1, 1, 2]], 
        w -> u\  + \ I\ v]; 

    X[u_, v_]\  = \ 
      Together[ComplexExpand[Re[Z[u, v]], 
\ \ \ \ \ \ \ \ \ \ \ 
          TargetFunctions -> {Re, \ Im}]]; 

    Y[u_, v_]\  = \ 
      Together[ComplexExpand[Im[Z[u, v]], 
\ \ \ \ \ \ \ \ \ \ \ 
          TargetFunctions -> {Re, \ Im}]]; 
Print["w = ", F[z]]; 

    Print["z = ", 
      solset[[1, 1, 2]]]; 

    Print["Z = ", Z[u, v]]; 
Print["X = ", X[u, v]]; 

    Print["Y = ", Y[u, v]]; \)

\(n\  = \ 8; 

    x[r_]\  = \ 
      Simplify[X[r\ Cos[\(k\ π\)\/n], r\ Sin[\(k\ π\)\/n]]]; 

    y[r_]\  = \ 
      Simplify[Y[r\ Cos[\(k\ π\)\/n], r\ Sin[\(k\ π\)\/n]]]; 

    eqns\  = \ N[Table[{x[r], y[r]}, {k, 0, n}]]; 

    ParametricPlot[Evaluate[eqns], {r, 0.05, 40}, 
      PlotRange -> {{\(-1.1\), 1.1}, {\(-1.1\), 1.1}}, 
      AspectRatio -> 1, DefaultColor -> RGBColor[0, 0, 1]]; \)

GoTo Chapter 10

Section 10.7  Two-Dimensional Fluid Flow

Suppose that a fluid flows over the complex plane and that the velocity at the point 
z = x + iy  is given by the vector  V(x,y)  =  p(x,y) + i q(x,y) .   We also require that the 
velocity does not depend on time and that the components  p(x,y)  and  q(x,y)  have 
continuous partial derivatives.  It is shown that p(x,y)  and  q(x,y)  are harmonic functions.

Theorem 10.5,  Page 352.  (Invariance of Flow)  Let  (Subscript[F, 1]) (w) = 
        Φ (u, v) + i Ψ (u, v)  denote the complex potential for a fluid flow in a domain  G  in the  w plane where the velocity is  (Subscript[V, 1])(u, v)= (OverBar[((Subscript[F, 1]))'(w)]).  If the function    w = S (z) = u (x, y) + i v (x, y)  is a one-to-one conformal mapping from a domain  D  in the  z plane onto  G,  then the composite function
(Subscript[F, 2]) (z) = 
        (Subscript[F, 1]) (S (z)) = 
          Φ (u (x, y), v (x, y)) + 
            i Ψ (u (x, y), v (x, y))  
is the complex potential for a fluid flow in  D  where the velocity is  (Subscript[V, 2])(x, y)= (OverBar[((Subscript[F, 2]))'(z)]).

Load Mathematica's  ComplexMap package.  Make sure this is done only ONCE during a Mathematica  session.

Get["Graphics`ComplexMap`"]

Example 10.22,  Page 353.  The complex potential  F (z) = (a + ib) z  has the velocity potential and stream function given by  ϕ (x, y) = a x - b y  and   φ (x, y) = b x + a y,  respectively.
Solution.  The solution is given in the textbook.

Example 10.23,  Page 354.  Consider the complex potential  F(z)  =  (A/2)(z^2)   where  
A  is a positive real number.  The velocity potential and stream function are given by:
u(x,y)  =  (A/2) ((x^2) - (y^2)),
v(x,y)  =  A x y.

\(Clear[A, f, F, x, y, u, v, z]; 
F[z_]\  = \ A\/2\ z\^2; 

    u[x_, y_]\  = \ 
      ComplexExpand[A\/2\ Re[\((x\  + \ I\ y)\)\^2]]; 

    v[x_, y_]\  = \ 
      ComplexExpand[A\/2\ Im[\((x\  + \ I\ y)\)\^2]]; 

    Print["The  complex potential is   F[z] = ", F[z]]; 

    Print["The velocity potential is u[x,y] = ", u[x, y]]; 

    Print["The   stream function  is v[x,y] = ", v[x, y]]; \)

For illustration purposes, choose A = 1.

\(\(\ f[z_]\  = \ \@z; 

    CartesianMap[f, {\(-4\), 4, 0.4}, {0, 4, 0.4}, AspectRatio -> 1, 
      Ticks -> {Range[\(-3\), 3, 1], Range[0, 3, 1]}, 
      PlotRange -> {{0, 2}, {0, 2}}]; 

    Print["        f[z] = ", f[z]]; \)\)

Example 10.24,  Page 354.  Find the complex potential for an ideal fluid flowing 
from left to right across the complex plane and around the unit circle  |z| = 1.

Solution:  Consider the complex potential  f(z)  =  A ( z + (1/z))   where  A  is a positive 
real number.  The velocity potential and stream function are given by:
u(x,y)  =  A ( x  + (x/((x^2)  +  (y^2)))),
v(x,y)  =  A ( y  - (y/((x^2)  +  (y^2)))).

\(Clear[A, f, F, g, solset, x, y, u, v, w, z]; 

    F[z_]\  = \ z\  + \ 1\/z; 

    u[x_, y_]\  = \ 
      ComplexExpand[Re[F[x + I\ y]], 
        TargetFunctions -> {Re, Im}]; 

    v[x_, y_]\  = 
      ComplexExpand[Im[F[x + I\ y]], 
        TargetFunctions -> {Re, Im}]; 

    Print["The  complex potential is   F[z] = ", F[z]]; 

    Print["The velocity potential is u[x,y] = ", u[x, y]]; 

    Print["The   stream function  is v[x,y] = ", v[x, y]]; \)

\(lines\  = \ Table[0 + 0.2  k, {k, 0, 7}]; 

    ContourPlot[v[x, y], {x, \(-2\), 2}, {y, 0, 2}, 
      Contours\  -> \ lines, PlotPoints -> 25, 
      AspectRatio -> 1\/2, ColorFunction -> Hue]; 

    Print["v[x,y] = ", v[x, y]]; 
Print["Contours = ", lines]; \)

The inverse function of  f(z)  = z + (1/z)  is  g(z)  =  ((z  +  (Sqrt[(z^2) - 4]) )/2).

\(Clear[f, g, w, z, Z0, Z1, Z2]; 
f[z_]\  = \ z\  + \ 1\/z; 

    solset\  = \ Solve[w == f[z], z]; 

    g[w_]\  = \ solset[[2, 1, 2]]; 

    Z0\  = \ PowerExpand[ExpandAll[g[f[z]]]]; 

    Z1\  = ReplaceAll[Z0\ , 
        \@\(1 - 2\ z\^2 + z\^4\)\  -> \ \(-\((1 - z\^2)\)\)]; 

    Z2\  = \ Together[Z1]; 
Print["f[z] = ", f[z]]; 

    Print["g[w] = ", g[w]]; 
Print[" "]; 

    Print["g[f[z]] = ", Z0]; 
Print["g[f[z]] = ", Z1]; 

    Print["g[f[z]] = ", Z2]; 
Print[" "]; 

    Print["f[g[w]] = ", f[g[w]]]; 

    Print["f[g[w]] = ", Together[f[g[w]]]]; \)

We can use it to find the image of the flow.

\(Clear[g1, g2, z]; 
g1[z_]\  = \ \ \(z\  + \ \@\(z\^2\  - \ 4\)\)\/2; 

    g2[z_]\  = \ \ \(z\  - \ \@\(z\^2\  - \ 4\)\)\/2; 

    graph1\  = \ 
      CartesianMap[g1, {0.00001, 3.00001, 0.2}, {0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; \ 

    graph2\  = \ 
      CartesianMap[
        g2, {\(-3.00001\), \(-0.00001\), 0.2}, {\ 0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; \ 

    Show[graph1, graph2, 
      Ticks -> {Range[\(-2\), 2, 0.5], Range[0, 2, 0.5]}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["        g(z) = ", g1[z]]; \)

Example 10.25,  Page 355.  Consider the complex potential  f(z) = A(Sqrt[(z^2)  +  1])
for an ideal fluid flowing from left to right across the complex plane slit along the 
segment from  -i  to  i.

\(Clear[A, f, F, lines, x, y, u, v, z]; 
F[z_]\  = \ \@\(z\^2\  + \ 1\); 
    
u[x_, y_]\  = \ 
      ComplexExpand[Re[F[x + I\ y]], 
        TargetFunctions -> {Re, \ Im, \ }]; 

    v[x_, y_]\  = \ 
      ComplexExpand[Im[F[x + I\ y]], 
        TargetFunctions -> {Re, \ Im, \ }]; 

    Print["The  complex potential is   F[z] = ", F[z]]\)\(Print["The velocity potential is u[x,y] = "]; 
Print[u[x, y]]; 

    Print["The   stream function  is v[x,y] = "]; 
Print[v[x, y]]; \)

\(lines\  = \ Table[0 + 0.2  k, {k, 0, 7}]; 

    ContourPlot[v[x, y], {x, 0.000001, 2}, {y, 0, 2}, 
      Contours\  -> \ lines, PlotPoints -> 25, 
      AspectRatio -> 1, ColorFunction -> Hue]; 

    Print["v[x,y] = ", v[x, y]]; \ 
Print["Contours = ", lines]; \)

The inverse function of  f(z)  =  (Sqrt[(z^2) + 1])  is  g(z)  =  (Sqrt[(z^2) - 1]).
We can use it to find the image of the flow.

Clear[g1, g2, z]\(g1[z_]\  = \ \ \@\(z\^2\  - \ 1\); 

    g2[z_]\  = \ \(-\@\(z\^2\  - \ 1\)\); 

    graph1\  = \ 
      CartesianMap[g1, {0.00001, 2.00001, 0.1}, {0.00001, 1.00001, 0.1}, 
        DisplayFunction -> Identity]; \ 

    graph2\  = \ 
      CartesianMap[
        g2, {\(-2.00001\), \(-0.00001\), 0.1}, {\ 0.00001, 1.00001, 0.1}, 
        DisplayFunction -> Identity]; \ 

    Show[graph1, graph2, 
      Ticks -> {Range[\(-2\), 2, 0.5], Range[0, 2, 0.5]}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["          g(z) = ", g1[z]]; \)

GoTo Chapter 10

Section 10.8  The Joukowski Airfoil

The function  J(z) = z + (1/z)  was studied by the Russian scientist N. E. Joukowski. 
It will be shown that the image of a circle passing through  (Subscript[z, 1]) = 1  and containing the point 
(Subscript[z, 2]) = -1  is mapped onto a curve that is shaped like the cross section of an airplane wing. 
We call this curve the Joukowski airfoil.

Load Mathematica's  ComplexMap package.
Make sure this is done only ONCE during a Mathematica  session.

Get["Graphics`ComplexMap`"]

The Joukowski Airfoil, Page 361.  Find the image of a fluid flow past the 
Joukowski airfoil.

\(Clear[F, F1, F2, F3, F4, dx, dy, r0, w, z, z1, z2, z3, z4]; 

    Δx\  = \ \(-\(1\/10\)\); 
Δy\  = \ \ 1\/5; 

    r0\  = \ \@\(\((1 - Δx)\)\^2 + 
          \((0 - Δy)\)\^2\); 

    z1\  = \ \ \(z\  + \ \@\(z\^2\  - \ 4\)\)\/2; 

    z2\  = \ \ \(z\  - \ \@\(z\^2\  - \ 4\)\)\/2; 

    z3\  = \ \(\(-z\)\  - \ \@\(z\^2\  - \ 4\)\)\/2; 

    z4\  = \ \(\(-z\)\  + \@\(\ z\^2\  - \ 4\)\)\/2; 

    w[z_]\  = \ \ 
      r0\ z\  + \ Δx\  + \ I\ Δy; 

    F[z_]\  = \ \ z\  + \ 1\/z; 
F1[z_]\  = \ \ F[w[z1]]; 

    F2[z_]\  = \ \ F[w[z2]]; 
F3[z_]\  = \ \ F[w[z3]]; 

    F4[z_]\  = \ \ F[w[z4]]; 
Print[\*""\!\(F\_1\)[z] = "", F1[z]]; 

    Print[\*""\!\(F\_2\)[z] = "", F2[z]]; 

    Print[\*""\!\(F\_3\)[z] = "", F3[z]]; 

    Print[\*""\!\(F\_4\)[z] = "", F4[z]]; \)

\(graph1\  = \ 
      CartesianMap[F1, {0.00001, 3.00001, 0.2}, {0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; \ 

    graph2\  = \ 
      CartesianMap[
        F2, {\(-3.00001\), 0.00001, 0.2}, {\ 0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; \ 

    graph3\  = \ 
      CartesianMap[F3, {0.00001, 3.00001, 0.2}, {0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; \ 

    graph4\  = \ 
      CartesianMap[
        F4, {\(-3.00001\), 0.00001, 0.2}, {\ 0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; \ 

    Show[graph1, graph2, graph3, graph4, 
\ \ \ \ \ 
      Ticks -> {Range[\(-2\), 2, 0.5], Range[\(-2\), 2, 0.5]}, 

      \ \ \ \ \ DisplayFunction -> $DisplayFunction]; 

    Print["f[z] = ", F1[z]]; \)

GoTo Chapter 10

Section 10.9  Schwarz-Christoffel Transformation

Theorem 10.6,  Page 370.  (Schwarz-Christoffel Formula)  Let  P  be a polygon in the w plane with vertices  (Subscript[w, 1]),(Subscript[w, 2]),⋯,(Subscript[w, n])  and exterior angles  (Subscript[α, 1]),(Subscript[α, 2]),⋯,(Subscript[α, n]), where -π<(Subscript[α, k])<π,  for  k = 1, 2, ⋯, n.  There exists a one-to-one conformal mapping  w = f (z)  from the upper half plane Im (z) > 0  onto  G  that satisfies the boundary conditions
	(Subscript[w, k]) = f ((Subscript[x, k]))   for   k = 1, 2, , n - 1    and(Subscript[w, n])=f(∞)     where   (Subscript[x, 1])<(Subscript[x, 2])<⋯<(Subscript[x, n - 1])<∞.
The derivative is   f'(z)= A (((z - (Subscript[x, 1])))^(-(Subscript[α, 1])/π)) 
          (((z - (Subscript[x, 2])))^(-(Subscript[α, 2])/π)) ⋯ 
          (((z - (Subscript[x, n - 1])))^(-(Subscript[α, n - 1])/π))
and the function  f(z)  can be expressed as an indefinite integral,
f (z)=B+A∫(((z - (Subscript[x, 1])))^(-(Subscript[α, 1])/π)) 
                (((z - (Subscript[x, 2])))^(-(Subscript[α, 2])/π)) 
                ⋯ 
                (((z - (Subscript[x, n - 1])))^(-(Subscript[α, n - 1])/
                      π)) dz
where  A  and  B  are suitably chosen constants.  Two of the points  {(Subscript[x, k])}  may be chosen arbitrarily, and the constants  A  and  B  determine the size and position of  P.

Load Mathematica's  ComplexMap package.  Make sure this is done only ONCE during a Mathematica  session.

Get["Graphics`ComplexMap`"]

Example 10.26,  Page 373.  Use the Schwarz Christoffel formula to verify that  
w  =  Arcsin(z)  maps the upper half plane  Im(z) > 0  onto the semi infinite strip  
((-π)/2) < u < (π/2),  v >0.

\(Clear[A, B, f1, F, z]; 
f1[z_]\  = \(\ A\)\/\@\(z\^2\  - \ 1\); 

    F[z_]\  = \ 
      B\  + \ ∫\(\(\ A\)\/\@\(z\^2\  - \ 1\)\) \[DifferentialD]z; 

    Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(\(\\ A\)\/\@\(z\^2\\  - \\ 1\)\)\[DifferentialD]z""]; 
Print["F[z] = ", F[z]]; \)

Remark.   No matter how tempting the integrand  B + A Log(z + (Sqrt[(z^2)  -  1]))
may seem, the Logarithm does not have the proper branch cut for this problem, and 
we must use the familiar integrand  B + A i ArcSin(z).  Switching the order of terms 
in the square root will solve this problem.

\(Clear[A, B, eq1, eq2, f, f1, F, z]; 

    f1[z_]\  = \ A\/\@\(1\  - \ z\^2\); 

    F[z_]\  = \ 
      B\  + \ ∫\(A\/\@\(1\  - \ z\^2\)\) \[DifferentialD]z; 

    Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(A\/\@\(1\\  - \\ z\^2\)\)\[DifferentialD]z""]; 
Print["F[z] = ", F[z]]; \)

Now solve for the coefficients  A  and  B.

\(Clear[eq1, eq2, f, z]; 
eq1\  = \ F[\(-1\)]\  == \(\ \(-π\)\)\/2; 
    
eq2\  = \ \ F[1]\  == \ \ π\/2; 

    solset\  = \ Solve[{eq1, eq2}, {A, B}]; 

    f[z_]\  = \ 
      ReplaceAll[F[z], solset[[1]]]; 

    Print[eq1]; 
Print[eq2]; 

    Print[solset[[1]]]; 

    Print["f[z] = B + A ArcSin[z]"]; 
Print["f[z] = ", f[z]]; \)

\(CartesianMap[f, {\(-2\), 2, 0.5}, {0.00001, 4.00001, 0.5}, 
      Ticks -> {Range[\(-1\), 1, 1], Range[0, 2, 1]}, 
      AxesLabel -> {"u", "v"}]; 

    Print["The mapping  w = f[z] = ", f[z]]; \)

Example 10.27,  Page 373.  Use the Schwarz Christoffel formula to verify that   
w  =  (Sqrt[(z^2) - 1])   maps the upper half plane  Im(z) > 0  onto the upper half plane
Im(w) > 0  slit along the segment from  0  to  i.

\(Clear[A, B, f1, F, z]; 

    f1[z_]\  = \(\ A\ z\)\/\(\@\(z\  + \ 1\)\ \@\(z\  - \ 1\)\); 

    F[z_]\  = \ 
      B\  + \ ∫\(\(\ A\ z\)\/\(\@\(z\  + \ 1\)\ \@\(z\  - \ 1\)\)\) 
            \[DifferentialD]z; 
Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(\(\\ A\\ z\)\/\(\@\(z\\  + \\ 1\)\\ \@\(z\\  - \\ 1\)\)\)\[DifferentialD]z""]; 

    Print["F[z] = ", F[z]]; \)

\(Clear[A, B, f1, F, z]; 
f1[z_]\  = \ \(A\ z\)\/\@\(z\^2\  - \ 1\); 

    F[z_]\  = \ 
      B\  + \ ∫\(\(A\ z\)\/\@\(z\^2\  - \ 1\)\) \[DifferentialD]z; 
    
Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(\(A\\ z\)\/\@\(z\^2\\  - \\ 1\)\)\[DifferentialD]z""]; 
Print["F[z] = ", F[z]]; \)

Now solve for the coefficients  A  and  B.

\(eq1\  = \ F[1] == 0; 
eq2\  = \ F[0] == I; 

    solset\  = \ Solve[{eq1, eq2}, {A, B}]\)\(f[z_]\  = \ 
      ReplaceAll[F[z], solset[[1]]]; 

    Print[eq1]; 
Print[eq2]; 

    Print[solset[[1]]]; 

    Print[\*""f[z] = B + A\!\(\@\(\(-1\) + z\^2\)\)""]; 

    Print["f[z] = ", f[z]]; \)

\(Clear[f1, f2, z]; 
f1[z_]\  = \ \ \ \@\(z\^2\  - \ 1\); 

    f2[z_]\  = \ I\ \@\(1\  - \ z\^2\); 

    graph1\  = \ 
      CartesianMap[f1, {0, 2, 0.2}, {0.00001, 1.00001, 0.1}, 
        PlotRange -> {{\(-2\), 2}, {0, 1.5}}, 
        DisplayFunction -> Identity]; \ 

    graph2\  = \ 
      CartesianMap[f2, {\(-2\), 0, 0.2}, {0.00001, 1.00001, 0.1}, 
        PlotRange -> {{\(-2\), 2}, {0, 1.5}}, 
        DisplayFunction -> Identity]; \ 

    Show[graph1, graph2, PlotRange -> {{\(-2\), 2}, {0, 1.5}}, 
      Ticks -> {Range[\(-2\), 2, 1], Range[0, 2, 1]}, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The mapping  w = f[z] = ", f[z]]; \)

Example 10.28,  Page 374.  Show that    
     w  =  (1/π)Arcsin(z)  + (i/π)Arcsin((1/z))  +  ((1  +  i)/2)
maps the upper half plane  Im(z) > 0  onto the channel in the right half plane
bounded by the coordinate axes and rays  x ≥ 1, y = 1  and  y ≥ 1, x =1.

\(Clear[A, B, f1, F]; 

    f1[z_]\  = \ \(A\ \((z\  - \ 1)\)\)\/\(z\ \@\((z\^2\  - \ 1)\)\); 

    F[z_]\  = \ 
      B\  + \ ∫\(\(A\ \((z\  - \ 1)\)\)\/\(z\ 
                  \@\((z\^2\  - \ 1)\)\)\) \[DifferentialD]z; 

    Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(\(A\\ \((z\\  - \\ 1)\)\)\/\(z\\ \@\((z\^2\\  - \\ 1)\)\)\)\[DifferentialD]z""]; 

    Print["F[z] = ", F[z]]; \)

This is one, formula for the integral.   However, we will use the following form 
of the integral to continue the computations.

\(F[z_]\  = \ 
      A\ \((ArcSin[z]\  + \ I\ ArcSin[1\/z])\)\  + \ B; 

    Print["F[z] = ", F[z]]; \)

Now solve for the coefficients  A  and  B.

\(eq1\  = \ F[\(-1\)]\  == \ 0; 

    eq2\  = \ F[1]\ \  == \ 1 + I; 

    solset\  = \ Solve[{F[\(-1\)] == 0, F[1] == 1 + I}, {A, B}]; 
    
f[z_]\  = \ 
      ReplaceAll[F[z], solset[[1]]]; 

    Print[eq1]; 
Print[eq2]; 

    Print[solset[[1]]]; 

    Print[\*""f[z] = B + A (ArcSin[z] + I ArcSin[(1/z)])""]; 
Print["f[z] = ", f[z]]; \)

\(CartesianMap[f, {\(-4\), 4, 0.5}, {0.00001, 2.00001, 0.2}, 
      Ticks -> {Range[0, 2, 0.5], Range[0, 2, 0.5]}]; 

    Print["The map  w = f[z] = ", f[z]]; \)

GoTo Chapter 10

Section 10.10  Images of a Fluid Flow

We have already examined several two-dimensional fluid flows and have discovered 
that the image of a flow under a conformal transformation is a flow.   The mapping  
w = u(x,y) + i v(x,y) ,  which is obtained by using the Schwarz-Christoffel formula, 
will allow us to find the streamlines for flows in domains in the w plane that are 
bounded by straight line segments.

Load Mathematica's  ComplexMap package.
Make sure this is done only ONCE during a Mathematica  session.

Get["Graphics`ComplexMap`"]

Example 10.29,  Page 380.  Show that the mapping  
     w  =  (1/π) ( (Sqrt[(z^2) - 1]) +  Log(z + (Sqrt[(z^2) - 1])) )
maps the upper half plane  Im(z) > 0  onto the domain in the upper half plane
that lies above the rays  u ≤ 0, v = 0  and  u ≥ 0, v = 0.

\(Clear[A, B, f, f1, F, z]; 

    f1[z_]\  = \ \(A\ \@\(z\  + \ 1\)\)\/\@\(z\  - \ 1\); 

    F[z_]\  = \ 
      B\  + \ ∫\(\(A\ \@\(z\  + \ 1\)\)\/\@\(z\  - \ 1\)\) 
            \[DifferentialD]z; 
Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(\(A\\ \@\(z\\  + \\ 1\)\)\/\@\(z\\  - \\ 1\)\)\[DifferentialD]z""]; 

    Print["F[z] = ", F[z]]; \)

\(f1[z_]\  = \ \(A\ \((z\  + \ 1)\)\)\/\@\(z\^2\  - \ 1\); 

    F[z_]\  = \ 
      B\  + \ ∫\(\(A\ \((z\  + \ 1)\)\)\/\@\(z\^2\  - \ 1\)\) 
            \[DifferentialD]z; 
Print["f[z] = ", f1[z]]; 

    Print[\*""F[z] = B + ∫\!\(\(A\\ \((z\\  + \\ 1)\)\)\/\@\(z\^2\\  - \\ 1\)\)\[DifferentialD]z""]; 

    Print["F[z] = ", F[z]]; \)

Now solve for the coefficients  A  and  B.

\(\(eq1\  = \ F[\(-1\)]\  == \ \ 0; 

    eq2\  = \ \ F[1]\  == \ \(-I\); 

    solset\  = \ Solve[{eq1, eq2}, {A, B}]; 

    f[z_]\  = \ 
      ReplaceAll[F[z], solset[[1]]]; 

    Print[eq1]; 
Print[eq2]; 

    Print[solset[[1]]]; 

    Print[\*""f[z] = B + A (Sqrt[-1 + (z^2)]) + A Log[z+(Sqrt[-1 + (z^2)])]""]; 
Print["f[z] = ", f[z]]; 
\)\)

Unfortunately, the branch cuts don't work out right and we will have to do this one by hand.

\(f[z_]\  = \ ReplaceAll[F[z], {A -> 1\/π, B -> 0}]; 

    Print["f[z] = ", f[z]]; \)

Using this formula in the book we form the two branches of the function.
But for the second branch, we cannot evaluate at  t = 0.

\(Clear[f1, f2, z]; 

    f1[z_]\  = \ \(\ \@\((z\^2 - 1)\) + Log[z + \@\(z\^2 - 1\)]\)\/π; 

    f2[z_]\  = \ \(\(-\@\((z\^2 - 1)\)\) + Log[z - \@\(z\^2 - 1\)]\)\/π; 
    
Print[\*""\!\(f\_1\)[z] = "", f1[z]]; 

    Print[\*""\!\(f\_2\)[z] = "", f2[z]]; \)

\(g1\  = \ 
      CartesianMap[f1, {0, 2, 0.5}, {0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; 

    g2\  = \ CartesianMap[
        f2, {\(-6.00001\), \(-0.00001\), 0.5}, {\ 0.00001, 1.00001, 0.2}, 
        DisplayFunction -> Identity]; 

    Show[g1, g2, DisplayFunction -> $DisplayFunction, 
      Ticks -> {Range[\(-1\), 1, 0.5], Range[0, 1, 0.5]}]; 

    Print["The mapping w = ", f[z]]\)

Here is the function that produces the desired result.
The graph is drawn parametrically in this case.

\(f[z_]\  = \ 
      \(I \@\( 1\  - \ z\^2\)\  + \ \(I\ π\)\/2\ 
           - \ I\ ArcSin[z]\)\/π; 

    eqns\  = \ 
      Table[{Re[f[x + I\ 0.2  k]], 
          Im[f[x + I\ 0.2  k]]}, {k, 1, 10}]; 

    ParametricPlot[Evaluate[eqns], {x, \(-9\), 4}, 
      PlotRange -> {{\(-1\), 2}, {\(-0.1\), 1.6}}, 
      Ticks -> {None, None}, AspectRatio -> 1.8\/3, 
      DefaultColor -> RGBColor[0, 0, 1]]; 

    Print["    f[z] = ", f[z]]; \)

GoTo Chapter 10

Section 10.11  Sources and Sinks

If the two-dimensional motion of an ideal fluid consists of an outward radial flow from a point and is symmetric in all directions, then the point is called a simple source.  For fluid flows, a sink is a negative source.

Load Mathematica's  ComplexMap package.
Make sure this is done only ONCE during a Mathematica  session.

Get["Graphics`ComplexMap`"]

Example 10.30,  Page 387.  Show that the complex potential for a fluid flowing 
from the source  at   z = +1  to the sink at   z = -1  is    F(z)  =  log(((z  -  1)/(z  +  1))).

\(Clear[f, F, W, x, y, z, Z]; 

    F[Z_]\  = \ Log[\(Z\  - \ 1\)\/\(Z\  + \ 1\)]; 

    solset\  = \ Solve[F[Z]\  == \ W, Z]; 

    f[x_, y_]\  = \ \ 
      ReplaceAll[solset[[1, 1, 2]], 
        W -> x\  + \ I\ y]; 

    Print["W = F[Z] = ", F[Z]]; 
Print["The inverse function is:"]; 
    
Print["f[x,y] = ", f[x, y]]; \)

\(n\  = \ 8; 

    eqns\  = \ 
      N[Table[{Re[f[x, \(\ k\ π\)\/n]], Im[f[x, \ \(k\ π\)\/n]]}, {k, 
            1, 2  n}]]; 

    ParametricPlot[Evaluate[eqns], {x, \(-4\), 4}, AspectRatio -> 1, 
      Ticks -> None, Axes -> None, 
      PlotRange -> {{\(-2.5\), 2.5}, {\(-2.5\), 2.5}}, 
      DefaultColor -> RGBColor[0, 0, 1]]; 

    Print["          f[x,y] = ", f[x, y]]; \)

Example 10.31,  Page 388.  Let two sources of unit strength be located at the points
z = 1 and z = -1. The resulting complex potential for a fluid flow is  F(z)  =  log((z^2) - 1).

\(Clear[f1, f2, F, W, x, y, z, Z]; 
F[Z_]\  = \ Log[Z\^2\  - \ 1]; 

    solset\  = \ Solve[F[Z]\  == \ W, Z]; 

    f1[x_, y_]\  = \ \ 
      ReplaceAll[solset[[1, 1, 2]], 
        W -> x\  + \ I\ y]; 

    f2[x_, y_]\  = \ \ 
      ReplaceAll[solset[[2, 1, 2]], 
        W -> x\  + \ I\ y]; 

    Print["W = F[Z] = ", F[Z]]; 
Print["The inverse function is:"]; 
    
Print["f[x,y] = ", f2[x, y]]; \)

\(n\  = \ 12; 

    eqns1\  = \ 
      N[Table[{Re[f1[x, \(k\ π\)\/n]], Im[f1[x, \(k\ π\)\/n]]}, {k, 
            1, n - 1}]]; 

    eqns2\  = \ 
      N[Table[{Re[f1[x, \(k\ π\)\/n]], Im[f1[x, \(k\ π\)\/n]]}, {k, 
            12 + 1, 12 + n - 1}]]; 

    eqns3\  = \ 
      N[Table[{Re[f2[x, \(k\ π\)\/n]], Im[f2[x, \(k\ π\)\/n]]}, {k, 
            1, n - 1}]]; 

    eqns4\  = \ 
      N[Table[{Re[f2[x, \(k\ π\)\/n]], Im[f2[x, \(k\ π\)\/n]]}, {k, 
            12 + 1, 12 + n - 1}]]; 

    eqns\ \  = \ Union[eqns1, eqns2, eqns3, eqns4]; 

    ParametricPlot[Evaluate[eqns], {x, \(-3\), 3}, AspectRatio -> 1, 
      Ticks -> None, PlotRange -> {{\(-2\), 2}, {\(-2\), 2}}, 
      DefaultColor -> RGBColor[0, 0, 1]]; 

    Print["             f[x,y] = ", f2[x, y]]; \)

Example 10.32,  Page 389.  Show that the complex potential for a fluid flow in the 
infinite strip ((-π)/2) < Re(z) < (π/2) with a single source at the origin is  F(z) = log(sin(z)).

\(Clear[F, W, x, y, z, Z]; 
F[Z_]\  = \ Log[Sin[Z]]; 

    solset\  = \ Solve[F[Z]\  == \ W, Z]; 

    f[x_, y_]\  = \ \ 
      ReplaceAll[solset[[1, 1, 2]], 
        W -> x\  + \ I\ y]; 

    Print["W = F[Z] = ", F[Z]]; 
Print["The inverse function is:"]; 
    
Print["f[x,y] = ", f2[x, y]]; \)

\(n\  = \ 12; 

    eqns\  = \ 
      Table[{Re[f[x, \(k\ π\)\/n]], Im[f[x, \(k\ π\)\/n]]}, {k, 1, 
          2  n}]; 

    ParametricPlot[Evaluate[eqns], {x, \(-3\), 3}, 
      AspectRatio -> 2\/1.6, Ticks -> None, Axes -> None, 
      PlotRange -> {{\(-1.6\), 1.6}, {\(-2\), 2}}, 
      DefaultColor -> RGBColor[0, 0, 1]]; 

    Print["         f[x,y] = ", f2[x, y]]; \)

Example 10.33,  Page 389.  Show that the complex potential for a fluid flow in the
infinite strip  ((-π)/2) < Re(z) < (π/2)  with a two sources at  z = ((-π)/2)  and  z = (π/2)  is   
     F(z)  =  log(((sin (z)  -  1)/(sin (z)  +  1))).

\(Clear[F, W, x, y, z, Z]; 

    F[Z_]\  = \ Log[\(Sin[Z] - 1\)\/\(Sin[Z] + 1\)]; 

    solset\  = \ Solve[F[Z]\  == \ W, Z]; 

    f[x_, y_]\  = \ \ 
      ReplaceAll[solset[[1, 1, 2]], 
        W -> x\  + \ I\ y]; 

    Print["W = F[Z] = ", F[Z]]; 
Print["The inverse function is:"]; 
    
Print["f[x,y] = ", f[x, y]]; \)

\(n\  = \ 8; 

    eqns\  = \ 
      Table[{Re[f[x, \(k\ π\)\/n]], Im[f[x, \(k\ π\)\/n]]}, 

        \ \ \ \ \ \ \ \ \ \ \ \ \ {k, 1, 2  n}]; 

    ParametricPlot[Evaluate[eqns], {x, \(-6\), 6}, 
      AspectRatio -> 2\/1.6, Ticks -> None, Axes -> None, 
      PlotRange -> {{\(-1.6\), 1.6}, {\(-2.4\), 2.4}}, 
      DefaultColor -> RGBColor[0, 0, 1]]; 

    Print["     f[x,y] = ", f[x, y]]; \)

GoTo Chapter 10