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

COMPLEX ANALYSIS: for Mathematics and Engineering, 
4th Edition,  2001,  ISBN: 0-7637-1425-9
Jones & Bartlett Publishers, Inc.  
40 Tall Pine Drive, Sudbury, MA 01776  
Tele.  (800) 832-0034, FAX:  (508)  443-8000  
E-mail:  mkt@jbpub.com  
Internet:  http://www.jbpub.com/  
This free software is compliments of the authors.
John H. Mathews,   mathews@fullerton.edu 
Russell W. Howell,  howell@westmont.edu

CHAPTER 11  FOURIER SERIES and the LAPLACE TRANSFORM

a = A = b = B = c = d = s = S = t = T = u = U = v = V = w = W = x = X = y = Y = z = Z = 
                                        0; 
Clear[a,a1,a2,a3,A,A1,A2,b,b1,b2,B1,B2,B,c,d,EqA,EqB,f,f1,f2,f12,F,F1,F2,F12,F3,g,G,h,H,i,int,j,k,L,Ly,L1,L2,Lf,LDf,ℒf,ℒDf,n,PQ,PR,r,res1,res2,s,solset,s1,s2,s3,S,S9,t,table,τ,T,θ,u,U,U1,U2,U7,v,v1,v2,v3,v4,val,V,w,w1,w2,W,x,X,y,y0,y1,Y,Yi,z,z1,z2,Z]; 
Needs["Graphics`FilledPlot`"]; 
Needs["Graphics`Colors`"];

Section 11.1	Fourier Series

Section 11.2	The Dirichlet Problem for the Unit Disk

Section 11.3	Vibrations in Mechanical Systems

Section 11.4	The Fourier Transform

Section 11.5	The Laplace Transform

Section 11.6	Laplace Transforms of Derivatives and Integrals

Section 11.7	Shifting Theorems and the Step Function

Section 11.8	Multiplication and Division by  t

Section 11.9	Inverting the Laplace Transform

Section 11.10  Convolution

Let  F (s)  and  G (s)  denote the transforms of  f (t)  and  g (t),  respectively.  Then the inverse of the product  F (s) G (s)  is given by the function  h (t) = (f*g) (t)  and is called the convolution of  f (t)  and  g (t)  and can be regarded as s generalized product of   f (t)  and  g (t).  Convolution will assist us in solving integral equations.

Theorem 10.24 (Convolution Theorem), Page 543. Let F (s) and G (s) denote the Laplace transforms of f (t) and g (t), respectively. Then the product  H (s) =F (s) G (s)  is the Laplace transform of the convolution of  f (t)  and  g (t),  and is denoted by  h (t) = (f*g) (t),  and has the integral representation 

		h[
          t] = (f*g)[t] = (∫_0^t) 
            f[τ] g[t - τ] dτ
or
h[
          t] = (g*f)[t] = (∫_0^t) 
            g[τ] 
              f[t - τ] dτ

Proof of Theorem 10.24, see text Page 543.

Example 11.29,  Page 545.  Use convolution to find the inverse Laplace transform of  H[s] = ((2  s)/((((s^2)   +   1))^2)).

Solution 11.29.

The transform H[s] is the product of  F[s] = 
        (1/((s^2)  +  1)),   and  G[
            s] = ((2  s)/((s^2)  +  1)), which are known to be the Laplace transforms of f[t] = Sin[t],   and  g[t] = 2 Cos[t],  respectively.  Enter H[s].

Remove[s, t]; 
Clear[f, F, g, G, h, H, int, τ, val]; 
F[s_]  =  (1/((s^2)  +  1)); 
G[s_]  =  ((2  s)/((s^2)  +  1)); 
H[s_]  =  ((2  s)/((((s^2)  +  1))^2)); 
Print["    F[s] = ", F[s]]; 
Print["    G[s] = ", G[s]]; 
Print["F[s]G[s] = ", F[s] G[s]]; 
Print[""]; 
Print["    H[s] = ", H[s]];

Use the convolution formula to construct h[
          t] = (f*g)[t] = (∫_0^t) 
            f[τ] g[t - τ] dτ.

f[t_]  =  Sin[t]; 
g[t_]  =  2 Cos[t]; 
int  =  ∫f[τ] 
            g[t - τ] \[DifferentialD]τ; 
h[
          t_]  =  (∫_0^t) f[τ] 
            g[t - τ] \[DifferentialD]τ; 
Print["      f[τ] = ", 
        f[τ]]; 
Print["      g[t] = ", g[t]]; 
Print["    g[t-τ] = ", g[t - τ]]; 
Print["f[τ]g[t-τ] = ", 
        f[τ] g[t - τ]]; 
Print[""]; 
Print["∫f[τ]g[t-τ]\[DifferentialD]τ = ", 
        int]; 
Print["∫f[τ]g[t-τ]\[DifferentialD]τ = ", 
        Expand[int]]; 
Print[""]; 
\(\(Print[\*""h[t] = (∫_0^t)f[τ]g[t-τ]\[DifferentialD]τ""];\)\ \)
Print["h[t] = ", h[t]]; 
h[t_]  =  Simplify[h[t]]; 
Print[""]; 
Print["h[t] = ", h[t]];

Aside. We can use the other convolution formula to construct h[
          t] = (g*f)[t] = (∫_0^t) 
            g[τ] f[t - τ] dτ  .

f[t_]  =  Sin[t]; 
g[t_]  =  2 Cos[t]; 
int  =  ∫g[τ] 
            f[t - τ] \[DifferentialD]τ; 
h[
          t_]  =  (∫_0^t) g[τ] 
            f[t - τ] \[DifferentialD]τ; 
Print["      f[t] = ", f[t]]; 
Print["    f[t-τ] = ", f[t - τ]]; 
Print["      g[τ] = ", g[τ]]; 
Print["g[τ]f[t-τ] = ", 
        g[τ] f[t - τ]]; 
Print[""]; 
Print["∫g[τ]f[t-τ]\[DifferentialD]τ = ", 
        int]; 
Print["∫g[τ]f[t-τ]\[DifferentialD]τ = ", 
        Expand[int]]; 
Print[""]; 
\(\(Print[\*""h[t] = (∫_0^t)g[τ]f[t-τ]\[DifferentialD]τ""];\)\ \)
Print["h[t] = ", h[t]]; 
h[t_]  =  Simplify[h[t]]; 
Print[""]; 
Print["h[t] = ", h[t]];

Aside. We can check this with Mathematica's  result using the LaplaceTransform package.

h[t_]  =  
        InverseLaplaceTransform[H[s],  s,  t]; 
Print["H[s] = ", H[s]]; 
Print[""]; 
Print["h[t] = ", h[t]];

Example 11.30,  Page 545.  Use the convolution theorem to solve the integral equation  f[t]  =  
        2  Cos[t]  -  (∫_0^t) 
            g[t - τ] f[τ] dτ.

Solution 11.30.

Use the facts that  ℒ(t)=(1/(s^2))  and  ℒ(Cos[t])=(s/((s^2) +  1)).  
Set up the integral equation and Laplace transform each term.

Remove[s, t]; 
Clear[eqn, f, F, solset]; 
eqn  =  
        F[s]  ==  ((2  s)/((s^2) +  1))  -   (1/(s^2)) F[
                s]; 
solset  =  Solve[eqn,  F[s]]; 
\(\(Print[\*""f[t] = 2 Cos[t] - (∫_0^t) g[t-τ]f[τ]\[DifferentialD]τ""];\)\ \)
\(\(Print[\*"" F[s] = \!\(\(2  s\)\/\(s\^2 + \\ 1\)\) -  (1/(s^2)) F[s]""];\)\ \)
Print[eqn]; 
Print[solset]; 
F[s_]  =  
        Apart[(Subscript[solset, [[1, 1, 2]]])]; 
Print[""]; 
Print["F[s] = ", F[s]];

Use known Laplace transforms to find the solution.

f1[t_] = t Sin[t]; 
f2[t_] = Cos[t]; 
F1[s_] = LaplaceTransform[f1[t], t,  s]; 
F2[s_] = 
        LaplaceTransform[f2[t], t, s,  ]; 
\(\(Print[\*""    \!\(f\_1\)[t] = "", 
        f1[t]];\)\ \)
\(\(Print[\*""ℒ(\!\(f\_1\)[t]) = "", 
        F1[s]];\)\ \)
Print[""]; 
\(\(Print[\*""    \!\(f\_2\)[t] = "", 
        f2[t]];\)\ \)
\(\(Print[\*""ℒ(\!\(f\_2\)[t]) = "", 
        F2[s]];\)\ \)
Print[""]; 
Print["F[s] = ", F[s]]; 
Print[""]; 
\(\(Print[\*""- ℒ((Subscript[f, 1])[t]) + 2 ℒ((Subscript[f, 1])[t]) = "", \(-\ F1[s]\) + 
          2  F2[s]];\)\ \)
Print[""]; 
\(\(Print[\*""f[t] = - \!\(f\_1\)[t] + 2 \!\(f\_2\)[t]   = "", \(-\ 
            f1[t]\) + 2  f2[t]];\)\ \)

Aside. We can check this with Mathematica's  result using the LaplaceTransform package.

f[t_]  =  
        InverseLaplaceTransform[F[s],  s,  t]; 
Print["F[s] = ", F[s]]; 
Print[""]; 
Print["f[t] = ", f[t]];

Engineers and physicists sometimes consider forces that produce large effects that are applied over a very short time interval. The force acting at the time an earthquake starts is an example. This leads to the idea of a unit impulse function δ(t).  Consider the small positive constant a, then the function (Subscript[δ, a])(t)  is defined by 

			(Subscript[δ, a])(t)={({{(1/a), for 0<t<a}, {0, otherwise        }})}.     

The unit impulse function is obtained by letting the interval width go to zero, i.e.,

			δ(t)=(lim_(a  ->  0)) (Subscript[δ, a])(t).  

Although δ(t) is called the Dirac  Delta Function, it is not an ordinary function. To be precise it is a distribution, and the theory of distributions permits manipulations of δ(t) as though it were a function. For our work, we will treat δ(t) as a function and investigate its properties.

Example 11.31,  Page 546. Find the Laplace transform of the Dirac delta function  δ[t].

Solution 11.31.

For any  a > 0  we define  (Subscript[δ, a])[t]  on the interval  [0, a]  as follows.

Remove[a, c, d, s, t]
Clear[f, F]; 
δ[a_, t_]  =  (1/a); 
F[a_, 
          s_]   =   (∫_0^a) δ[a, 
              t] (E^(-s t)) \[DifferentialD]t; 
Print["δ[a,t] = ", δ[a, 
          t], "  for  0 ≤ t ≤ a"]; 
Print[""]; 
\(\(Print[\*""F[a,s] = (∫_0^a)δ[a,t](E^(-st))\[DifferentialD]t""];\)\ \)
\(\(Print[\*""F[a,s] = (∫_0^a)"", δ[a, 
          t], \*""(E^(-st))\[DifferentialD]t""];\)\ \)
Print["F[a,s] = ", F[a, s]]; 
Print[""]; 
F[a_, s_]   = Together[F[a, s]]; 
Print["F[a,s] = ", F[a, s]];

Plot  (Subscript[δ, a])[t]  for  a = (1/10), (1/20), (1/50), (1/100).

A = {(1/10), (1/20), (1/50), 
          (1/100)}; 
\(\(Do[
        c = A\_\(\([[\)\(i\)\(]]\)\); \ 
FilledPlot[δ[c, t], {t, 0, c}, 
          PlotRange -> {{0, 0.1}, {0, 
                100}}, 
PlotStyle -> Magenta, 
          Fills -> Magenta, AspectRatio -> 1, 
          Ticks -> {{1\/10, 1\/20, 1\/50, 1\/100}, Range[0, 100, 20]}, 
          AxesLabel -> {"t", "y"}]; \ 
Print["δ[", 
          c, ",t] = ", δ[c, 
            t], "  for  0 ≤ t ≤ ", c];\ , {i, 1, 
          4}];\)\ \)

Find  (lim_(a -> 0)) ℒ((Subscript[δ, a])[t])  to obtain the Laplace transform of  δ[t].

δ[t_]  =  
        DiracDelta[t]; 
Print["   δ[a,t] = ", δ[a, 
          t], "  for  0 ≤ t ≤ a"]; 
Print["ℒ(δ[a,t]) = ", 
        F[a, s]]; 
Print[""]; 
Print["     δ[t] = ", δ[
          t]]; 
\(\(Print[\*""     δ[t] = (lim_(a -> 0)) "", δ[a, 
          t], "  for  0 ≤ t ≤ a"];\)\ \)
Print[""]; 
\(\(Print[\*""  ℒ(δ[t]) = (lim_(a -> 0)) "", F[a, s]];\)\ \)
Print["  ℒ(δ[t]) = ", 
        Limit[F[a, s], a -> 0]];

Aside. We can check this with Mathematica's  result using the LaplaceTransform package.

δ[t_]  =  
        DiracDelta[t]; 
F[s_]  =  LaplaceTransform[δ[t],  t,  s]; 
Print["δ[t] = ", δ[
          t]]; 
Print[""]; 
Print["F[s] = ", F[s]];

Example 11.32,  Page 547.  Solve the initial value problem y''[t]  + 4 y'[t]  + 13 y[t] =3 δ[t] with y[0]=0 and y'[0]=0.

Solution 11.32.

Taking transforms results in the following equation.

Remove[s, t, y, Y]; 
Clear[de, f, F, eqn, solset]; 
f[t_]  =  3  DiracDelta[t]; 
F[s_]  =  
        LaplaceTransform[f[t],  t,  s]; 
de  =  y''[t]  + 4 y'[t]  + 13 y[t]  ==  
          f[t]; 
y0 = 0; 
y1 = 0; 
eqn = 
        ReplaceAll[
          de, {f[t] -> F[s], 
            y[t] -> Y[s], y'[t] -> 
              s Y[s] - y0, y''[t] -> 
              (s^2) Y[s] - s y0 - y1}]; 
solset  =  Solve[eqn,  Y[s]]; 
Print[de]; 
Print[" y[0] = ", y0]; 
Print["y'[0] = ", y1]; 
Print[""]; 
Print[eqn]; 
Print[""]; 
Print[solset]; 
Y[s_]  =  
        (Subscript[solset, [[1, 1, 2]]]); 
Print["Y[s] = ", Y[s]];

From a table of Laplace transforms we get  ℒ((e^(-2 t)) 
              Sin[3  t])=(3/((((s + 2))^2) + (3^2))) = (3/(13 + 4 s + (s^2))), and the solution is  y[t]  =  (e^(-2 t)) Sin[3  t].

Aside. We can check this with Mathematica's  result using the LaplaceTransform package.

y[t_]  =  
        InverseLaplaceTransform[Y[s],  s,  t]; 
Print["Y[s] = ", Y[s]]; 
Print[""]; 
Print["y[t] = ", y[t]]; 
y[t_] = ComplexExpand[y[t]]; 
Print["y[t] = ", y[t]]; 
y[t_] = TrigReduce[y[t]]; 
Print[""]; 
Print["y[t] = ", y[t]];

Aside. We use Mathematica/s DSolve package to solve the D.E. and plot the solution.
Notice.  The initial condition for Mathematica must be changed to y'[0] = 3.  This is required because  \(\(y'\)[\(0\^-\)] = \(0\ \ and\ \ \(y'\)[\(0\^+\)] = 3\)\).

Clear[t, y]; 
de  =  y''[t]  + 4 y'[t]  + 13 y[t]  == 
          3  DiracDelta[t]; 
ics = {y[0] == 0,  y'[0] == 3}; 
solset  =  
        DSolve[Flatten[{de,  ics}],  y[t],  
          t]; 
Print[de]; 
Print[ics]; 
Print[solset]; 
y[t_]  =  
        (Subscript[solset, [[1, 1, 2]]]); 
\(\(Plot[y[t], {t, 0, 2.5}, 

        PlotRange -> {{0, 2.5}, {\(-0.06\), 0.45}}, 
        PlotStyle -> Magenta, AxesLabel -> {"t", "y"}, 
        Ticks -> {Range[0, 2.5, 0.5], 
            Range[0, 0.4, 0.1]}];\)\ \)
Print["y[t] = ", y[t]];

Theorem 10.25 (Initial value Problem - IVP convolution method), Page 548. The unique solution to the initial value problem  

	a y''[t]  +  b y'[t]  +  c y[t] = g[t]  with  y[0] = (Subscript[y, 0]),   and  y'[0] = (Subscript[y, 1])  

is given by  

	y[t]  =  u[t]  +  (h*g)[t],  

where  u[t]  is the solution to the homogeneous equation  a u''[t]  +  b u'[t]  +  c u[t]  =  0  with  u[0] = (Subscript[y, 0]),   and  u'[0] = (Subscript[y, 1]),  and  h[t]  has the Laplace transform given by  H[s] = (1/(a (s^2)  +  b s  +  c)).

Proof of Theorem 10.25, see text Page 548.

Example 11.33, Page 549. Use convolution to solve the initial value problem y''[t]  +  y[t]  = Tan[t] with y[0]=1 and y'[0]=2.

Solution 11.33.

The solution is accomplished in two steps. (a) First solve  u''[t]  +  u[t]  = 0  with  u[0] = 1  and  u'[0] = 2.  
(b)  Second solve for  v''[t]  +  v[t]  = Tan[t]  with  v[0] = 0  and  v'[0] = 0.  Then the desired solution is  y[t]  =  u[t]  +  v[t].  

(a) First solve  u''[t]  +  u[t]  = 0  with  u[0] = 1  and  u'[0] = 2.

Remove[s, t, u, U, v, y]; 
Clear[de, eqn, f, g, H, int, solset, τ, v1, v2, v3, v4, 
        val]; 
f[t_]  =  Tan[t]; 
de  =  u''[t]  +  u[t]  ==  
          0; 
u0 = 1; 
u1 = 2; 
eqn = 
        ReplaceAll[
          de, {u[t] -> U[s], u'[t] -> 
              s U[s] - u0, u''[t] -> 
              (s^2) U[s] - s u0 - u1}]; 
solset  =  Solve[eqn,  U[s]]; 
Print[de]; 
Print[" u[0] = ", u0]; 
Print["u'[0] = ", u1]; 
Print[""]; 
Print[eqn]; 
Print[""]; 
Print[solset]; 
U[s_]  =  
        (Subscript[solset, [[1, 1, 2]]]); 
u[t_]  =  InverseLaplaceTransform[U[s], s, t]; 
Print["U[s] = ", U[s]]; 
U[s_]  = 
        Map[Simplify, ExpandAll[U[s]]]; 
Print["U[s] = ", U[s]]; 
Print[" "]; 
Print[" u[t] = ", u[t]]; 
Print[" u[0] = ", u[0]]; 
Print["u'[0] = ", u'[0]];

(b)  Second solve for  v''[t]  +  v[t]  = Tan[t]  with  v[0] = 0  and  v'[0] = 0. 
The convolution theorem is used to make the construction of  v[t] = (h*g)[t]  where  H[s] = (1/((s^2)  +  1))  and  g[t]  = Tan[t].

H[s_]  =  (1/((s^2)  +  1)); 
h[t_]  =  InverseLaplaceTransform[H[s], s, t]; 
g[t_]  =  Tan[t]; 
int  =  ∫h[t - τ] 
            g[τ] \[DifferentialD]τ; 
v[
          t_]  =  (∫_0^t) h[t - τ] 
            g[τ] \[DifferentialD]τ; 
Print["      h[t] = ", h[t]]; 
Print["      g[t] = ", g[t]]; 
Print["h[t-τ]g[τ] = ", 
        h[t - τ] g[τ]]; 
Print["∫h[t-τ]g[τ]\[DifferentialD]τ = ", 
        int]; 
Print[" "]; 
\(\(Print[\*"" v[t] = (∫_0^t)h[t-τ]g[τ]\[DifferentialD]τ""];\)\ \)
Print[" v[t] = ", v[t]]; 
Print[" v[0] = ", v[0]]; 
Print["v'[0] = ", v'[0]];

Now guide Mathematica using  trigonometric substitutions to simplify the expression for  v[t].
Notice.  This cannot be done automatically.

v1  =  
        Cos[t] (Log[Cos[(t/2)] - Sin[(t/2)]] - 
                Log[Cos[(t/2)] + Sin[(t/2)]]) + Sin[t]; 
v2  =  
        Cos[t] Log[((Cos[(t/2)] - Sin[(t/2)])/(Cos[(t/2)] + Sin[(t/2)]))] + 
          Sin[t]; 
v3  =  
        Cos[t] Log[((Cos[(t/2)]^2 - Sin[(t/2)]^2)/(((Cos[(t/2)] + Sin[(t/2)]))^2))] + Sin[t]; 
v4  =  
        Cos[t] Log[Cos[t]/(1 + 2 Cos[(t/2)] Sin[(t/2)])] + 
          Sin[t]; 
v[t_]  =  
        Cos[t] Log[Cos[t]/(1 + Sin[t])] + Sin[t]; 
y[t_]  =  u[t]  +  v[t]; 
Print[" u[t] = ", u[t]]; 
Print[" u[0] = ", u[0]]; 
Print["u'[0] = ", u'[0]]; 
Print[" v[t] = ", v[t]]; 
Print[" "]; 
Print[" y[t] = u[t] + v[t]"]; 
Print["y [t] = ", y[t]]; 
Print[" y[0] = ", y[0]]; 
Print["y'[0] = ", y'[0]];

Aside. We can substitute y[t] into the D.E.'s and verify that it is the solution.

\(
\(Print["y''[t]+y[t] = ", \(y''\)[t] + 
          y[t]];\)\ \)
Print[" "]; 
Print["y''[t]+y[t] = ", 
        Simplify[y''[t] + y[t]]];

Aside. We use Mathematica/s DSolve package to solve the D.E. and plot the solution.

Clear[t, Y]; 
de  =  Y''[t]  +  Y[t]  == 
          Tan[t]; 
ics = {Y[0] == 1,  Y'[0] == 2}; 
solset  =  
        DSolve[Flatten[{de,  ics}],  Y[t],  
          t]; 
Print[de]; 
Print[ics]; 
Print[solset];

Remark. This involves the version unsimplified of v[t].  So we will plot our previously simplified function y[t].

Plot[y[t], {t, 0, (π/2)}, 
        PlotRange -> {{0, (π/2)}, {0, 
              3}}, 
Ticks -> {Range[
              0, (π/2), (π/4)], Range[0, 3, 1]}, 
        PlotStyle -> Magenta, 
        AxesLabel -> {"t", "y"}]; 
Print["y[t] = ", y[t]];

Section 11.10  Exercises for Convolution
See textbook page 549.

GO TO  THE  Table of Contents