Numerical Computation of Internal and External Flows.

By C.Hirsh

Chapter 4

Page 175 Hirsch

Taylor series in equation 4.2  .30 is obtained from,  
    see note with this page:

Series[ 
      2 (((1 + ((δ^(2 ))/4)))^(((-1)/2))) 
        ArcSinh[(δ/2)], {δ, 0, 7}]

(36/(7!))

Page 176

My result:

((-n)/((6 24 64))) (2 (n - 2) (((5  n + 22))/90) + 
            2 (2  n - 1) (3/10) - (45/7)) // Simplify

Books result:

((-n)/(4^5)) ((5/7) + (((n - 1))/5) + (n - 
                  1) (((n - 2))/(3^5)))//Simplify

Page 177-178

Equation 4.2.41

Series[
(((u[x + 2  h] - 2  u[x + h] + u[x]))/(h^2)), 
{h, 0, 1}]

Equation 4.2.42

Series[
(((u[x] - 2  u[x - h] + u[x - 2  h]))/(h^2)), 
{h, 0, 1}]

Equation 4.2.43

Series[
(((u[x + h] - 2  u[x] + u[x - h]))/(h^2)), 
{h, 0, 2}]

Equation 4.2.44

Series[
(((u[x + (3/2) h] - u[x + (h/2)] + -u[x - (h/2)] + u[x - 3  (h/2)]))/(2  (h^2))), 
{h, 0, 
        2}]

Equation 4.2.45

Series[
(((2  u[x] - 5  u[x + h] + 4  u[x + 2  h] - u[x + 3  h]))/(h^2)), 
{h, 0, 2}]

Equation 4.2.46

Series[
(((2  u[x] - 5  u[x - h] + 4  u[x - 2  h] - u[x - 3  h]))/(h^2)), 
{h, 0, 2}]

Equation 4.2.47

Series[
(((-u[x + 2  h] + 16  u[x + h] - 30  u[x] + 16  u[x - h] - u[x - 2  h]))/(12  (h^2))), 
{h, 0, 
        4}]

Equation 4.2.48

Series[(1/(4 12  (h^2))) (-5 u[x + (5/2) h] + 
            39  u[x + (3/2) h] - 34  u[x + (h/2)] - 34  u[x - (h/2)] + 
            39  u[x - (3/2) h] - 
            5  u[x - (5/2) h]), 
{h, 0, 4}]

Equation 4.2.50

Series[(((k[x + (h/2)] (u[x + h] - u[x]) - k[x - (h/2)] (u[x] - u[x - h])))/(h^2)), {h, 0, 2}]

Page 179

Equation 4.2.52

Series[((u[x + 3  h] - 3  u[x + 2  h] + 3  u[x + h] - u[x])/(h^3)), {h, 
        0, 2}]

2  Function[x, (x^3) - (3/2) (x^4)][e - 1] // Expand

Equation 4.2.53

Series[((-3 u[x + 4  h] + 14  u[x + 3  h] - 24  u[x + 2  h] + 18  u[x + h] - 5  u[x])/(2  (h^3))), {h, 0, 2}]

Page 184

The coefficients in 4.3.15 become.

\(Solve[{\(a\_+\) + a\_0 + \(a\_-\) == 
          0, 
Δx \((\(a\_+\) - \(a\_-\))\) + \(b\_+\) + b\_0 + \(b\_-\) == 
          0, 
\(Δx\^2\/2\) \((\(a\_+\) + \(a\_-\))\) + Δx \((\(b\_+\) - \(b\_-\))\) + \(c\_+\) + 
            c\_0 + \(c\_-\) == 
          0, 
\(Δx\^3\/6\) \((\(a\_+\) - \(a\_-\))\) + \(Δx\^2\/2\) \((\(b\_+\) + \(b\_-\))\) + Δx \((\(c\_+\) - \(c\_-\))\) == 
          0}, 
{\(a\_+\), a\_0, \(a\_-\), b\_0}]\)

The truncation error for 4.3.10 will then become.

lhs=Series[
(Subscript[a, 0]) u[0]+\(\(a\_-\)\ u[\(-Δx\)]\)+\(\(a\_+\)\ u[Δx]\)+(Subscript[b, 0]) u'[0]+\(b\_-\) u'[-Δx]+\(b\_+\) u'[Δx]+(Subscript[c, 0]) u''[0]+\(c\_-\) u''[-Δx]+\(c\_+\) u''[Δx],
{Δx, 0, 7}]

And upon simplifying to make the first four terms zero one gets:

lhs /. Out[1] // Simplify

Coefficient[lhs /. Out[1], u''''[0]]

Coefficient[lhs /. Out[1], u''''''''[0]]

(((Δ((x^7)/(8!))))^(-1)) Out[10] // 
      Expand

Not the same!

\(Solve[{2 \((\(b\_+\) - \(b\_-\))\) + \(10\/Δx\) \((\(c\_+\) + \(c\_-\))\) - \(2\/Δx\) c\_0 == 
          0, 
2 \((\(b\_+\) + \(b\_-\))\) + \(14\/Δx\) \((\(c\_+\) - \(c\_-\))\) == 
          0, 
4 \((\(b\_+\) - \(b\_-\))\) + \(28\/Δx\) \((\(c\_+\) + \(c\_-\))\) - \(2\/Δx\) 
              c\_0 == 0}, 
{\(b\_+\), \(b\_-\), 
        c\_0}]\)

Page 198.  Problem 4.3

Solve[{
a + b + c + d + e == 
          0, 
2  a + b - d - 2  e == 
          1, 
4  a + b + d + 4  e == 
          0, 
8  a + b - d - 8  e == 
          0, 
16  a + b + d + 16  e == 0}, {a, b, c, 
        d, e}]

Residual of the method becomes:

32  a + b - d - 32  e /. Out[18] // Simplify

4/5!

Solve[{
a + b + c + d + e == 
          0, 
2  a + b - d - 2  e == 
          0, 
4  a + b + d + 4  e == 
          2, 
8  a + b - d - 8  e == 
          0, 
16  a + b + d + 16  e == 0}, {a, b, c, 
        d, e}]

64  a + b + d + 64  e /. Out[22] // Simplify

8/6!

Page 198.  Problem 4.4

Solve[{
a + b + c + d + e == 
          0, 
2  a + b - d - 2  e == 
          0, 
4  a + b + d + 4  e == 
          0, 
8  a + b - d - 
            8  e == 3!, 
16  a + b + d + 
            16  e == 0}, {a, b, c, d, e}]

32  a + b - d - 32  e /. Out[25]

30/5!

Solve[{
a + b + c + d + e == 
          0, 
2  a + b - d - 2  e == 
          0, 
4  a + b + d + 4  e == 
          0, 
8  a + b - d - 8  e == 
          0, 
16  a + b + d + 16  e == 4!}, {a, 
        b, c, d, e}]

(32  a + b - d - 32  e)/5! /. Out[28]

(64  a + b + d + 64  e)/6! /. Out[28]

Page 199.  Problem 4.5

Δx = 0.1;
W = 7;

e = 
        Table[(E^x), {x, -W Δx, 
            W Δx, Δx}];
s = 
        Table[Sin[x], {x, -W Δx, 
            W Δx, Δx}];
c = 
        Table[Cos[x], {x, -W Δx, 
            W Δx, Δx}];

First derivative.

der1Exact[e] = 1;
der1Exact[s] = 1;
der1Exact[c] = 0;

der1Back1 = (((# - RotateRight[#]))/Δ)x &;
der1Back2 = (((3  # - 4  RotateRight[#] + RotateRight[#, 2]))/(2  Δx)) &;
der1Fowd1 = (((RotateLeft[#] - #))/Δ)x &;
der1Fowd2 = (((-3 # + 4  RotateLeft[#] - RotateLeft[#, 2]))/(2  Δx)) &;
der1Cent2 = (((RotateLeft[#] - RotateRight[#]))/(2 Δx)) &;
der1Cent4 = (((-RotateLeft[#, 2] + 8  RotateLeft[#] - 8  RotateRight[#] + RotateRight[#, 2]))/(12 Δx)) &;

der1TrunkBack1 = Δx D[#, {x, 2}] /. {x -> 
                0. } &;
der1TrunkBack2 = Δ((x^2)/3) 
              D[#, {x, 3}] /. {x -> 0. } &;
der1TrunkFowd1 = -Δx D[#, {x, 2}]/2 /. {x -> 
                0. } &;
der1TrunkFowd2 = Δ((x^2)/3) 
              D[#, {x, 3}] /. {x -> 0. } &;
der1TrunkCent2 = ((-Δ(x^2))/6) 
              D[#, {x, 3}] /. {x -> 0. } &;
der1TrunkCent4 = Δ((x^4)/30) 
              D[#, {x, 5}] /. {x -> 0. } &;

der1Approx = {der1Back1[#][[W + 1]], der1Back2[#][[W + 
                  1]], 
der1Fowd1[#][[W + 
                  1]], der1Fowd2[#][[W + 
                  1]], 
der1Cent2[#][[W + 
                  1]], der1Cent4[#][[W + 1]]} &;

der1Trunk = {der1TrunkBack1[#], 
            der1TrunkBack2[#], 
der1TrunkFowd1[#], 
            der1TrunkFowd2[#], 
der1TrunkCent2[#], 
            der1TrunkCent4[#]} &;

Calculations for three different functions.

\(\(Print[\*""--Approximation/Exact for first der of (E^x) at 0.--""];\)\)
Print["Exact first derivative at 0."];
der1Exact[e]
Print["Approx first derv at 0. for several methods."];
der1Approx[e]
Print["Difference between exact and approx."];
der1Exact[e] - der1Approx[e]
Print["Truncation error for each method."];
der1Trunk[(E^x)]

\(\(Print["--Approximation/Exact for first der of Cos[x] at 0.--"];\)\)
Print["Exact first derivative at 0."];
der1Exact[c]
Print["Approx first derv at 0. for several methods."];
der1Approx[c]
Print["Difference between exact and approx."];
der1Exact[c] - der1Approx[c]
Print["Truncation error for each method."];
der1Trunk[Cos[x]]

\(\(Print["--Approximation/Exact for first der of Sin[x] at 0.--"];\)\)
Print["Exact first derivative at 0."];
der1Exact[s]
Print["Approx first derv at 0. for several methods."];
der1Approx[s]
Print["Difference between exact and approx."];
der1Exact[s] - der1Approx[s]
Print["Truncation error for each method."];
der1Trunk[Sin[x]]

Conclusion: The truncation error in most cases ( exect when due to the special functions considered ) is a very good approximation to the error.  It seems to be almost always an over estimate.  But the trend is that if you know the truncation error you know the error involved.

Second derivative.

der2Exact[e] = 1;
der2Exact[s] = 0;
der2Exact[c] = -1;

der2Back1 = (((# - 2  RotateRight[#] + RotateRight[#, 2]))/Δ)(x^2) &;
der2Back2 = (((2  # - 5  RotateRight[#] + 4  RotateRight[#, 2] - RotateRight[#, 3]))/Δ)(x^2) &;
der2Fowd1 = (((RotateLeft[#, 2] - 2  RotateLeft[#] + #))/Δ)(x^2) &;
der2Fowd2 = (((2  # - 5  RotateLeft[#] + 4  RotateLeft[#, 2] - RotateLeft[#, 3]))/Δ)(x^2) &;
der2Cent2 = (((RotateLeft[#] - 2  # + RotateRight[#]))/Δ)(x^2) &;
der2Cent4 = (((-RotateLeft[#, 2] + 16  RotateLeft[#] - 30  # + 16  RotateRight[#] - RotateRight[#, 2]))/(12 Δ(x^2))) &;

der2TrunkBack1 = Δx D[#, {x, 3}] /. {x -> 
                0. } &;
der2TrunkBack2 = ((-11)/12) Δ(x^2) 
              D[#, {x, 4}] /. {x -> 0. } &;
der2TrunkFowd1 = -Δx D[#, {x, 3}] /. {x -> 
                0. } &;
der2TrunkFowd2 = (11/12) Δ(x^2) 
              D[#, {x, 4}] /. {x -> 0. } &;
der2TrunkCent2 = ((-Δ(x^2))/12) 
              D[#, {x, 4}] /. {x -> 0. } &;
der2TrunkCent4 = Δ((x^4)/90) 
              D[#, {x, 6}] /. {x -> 0. } &;

der2Approx = {der2Back1[#][[W + 1]], der2Back2[#][[W + 
                  1]], 
der2Fowd1[#][[W + 
                  1]], der2Fowd2[#][[W + 
                  1]], 
der2Cent2[#][[W + 
                  1]], der2Cent4[#][[W + 1]]} &;

der2Trunk = {der2TrunkBack1[#], 
            der2TrunkBack2[#], 
der2TrunkFowd1[#], 
            der2TrunkFowd2[#], 
der2TrunkCent2[#], 
            der2TrunkCent4[#]} &;

Calculations for three different functions.

\(\(Print[\*""--Approximation/Exact for second der of (E^x) at 0.--""];\)\)
Print["Exact second derivative at 0."];
der2Exact[e]
Print["Approx second derv at 0. for several methods."];
der2Approx[e]
Print["Difference between exact and approx."];
der2Exact[e] - der2Approx[e]
Print["Truncation error for each method."];
der2Trunk[(E^x)]

Print["--Approximation/Exact for second der of Cos[x] at \
0.--"];
Print["Exact second derivative at 0."];
der2Exact[c]
Print["Approx second derv at 0. for several methods."];
der2Approx[c]
Print["Difference between exact and approx."];
der2Exact[c] - der2Approx[c]
Print["Truncation error for each method."];
der2Trunk[Cos[x]]

Print["--Approximation/Exact for second der of Sin[x] at \
0.--"];
Print["Exact second derivative at 0."];
der2Exact[s]
Print["Approx second derv at 0. for several methods."];
der2Approx[s]
Print["Difference between exact and approx."];
der2Exact[s] - der2Approx[s]
Print["Truncation error for each method."];
der2Trunk[Sin[x]]

Conclusion: Again the truncation error is the same order of magnitude as the actual error.  Look at the order of magnitude of the truncation error ( given a particular value of Δx ) and that is the same magnitude at the absolute error.

Page 199.  Problem 4.6

The coefficients in 4.3.15 become.

Clear[a, b, c, e, Δx];
\(Solve[{\(a\_+\) + a\_0 + \(a\_-\) == 
          0, 
Δx \((\(a\_+\) - \(a\_-\))\) + \(b\_+\) + b\_0 + \(b\_-\) == 
          0, 
\(Δx\^2\/2\) \((\(a\_+\) + \(a\_-\))\) + Δx \((\(b\_+\) - \(b\_-\))\) + \(c\_+\) + 
            c\_0 + \(c\_-\) == 
          0, 
\(Δx\^3\/6\) \((\(a\_+\) - \(a\_-\))\) + \(Δx\^2\/2\) \((\(b\_+\) + \(b\_-\))\) + Δx \((\(c\_+\) - \(c\_-\))\) == 
          0}, 
{\(a\_+\), \(b\_+\), \(c\_+\), 
        a\_0}]\)

The truncation error for 4.3.10 will then become.

lhs=Series[
(Subscript[a, 0]) u[0]+\(\(a\_-\)\ u[\(-Δx\)]\)+\(\(a\_+\)\ u[Δx]\)+(Subscript[b, 0]) u'[0]+\(b\_-\) u'[-Δx]+\(b\_+\) u'[Δx]+(Subscript[c, 0]) u''[0]+\(c\_-\) u''[-Δx]+\(c\_+\) u''[Δx],
{Δx, 0, 7}]

And upon simplifying to make the first four terms zero one gets:

lhs /. Out[50] // Simplify

This needs to be grouped in terms of the derivatives of the unknown function u.

Coefficient[Out[52], u''''[0]] // 
      Simplify
Coefficient[Out[52], u'''''[0]] // 
      Simplify
Coefficient[Out[52], u''''''[0]] // 
      Simplify
Coefficient[Out[52], u'''''''[0]] // 
      Simplify
Coefficient[Out[52], u''''''''[0]] // Simplify

\(Solve[{8\ Δx\^2\ \(a\_-\) - 
            20\ Δx\ \(b\_-\) + 
            36\ \(c\_-\) - Δx\ b\_0 + 6\ c\_0 == 
          0, 
16\ Δx\^2\ \(a\_-\) - 
            32\ Δx\ \(b\_-\) + 32\ \(c\_-\) - 
            3\ Δx\ b\_0 + 16\ c\_0 == 
          0, 
16\ Δx\^2\ \(a\_-\) - 
            36\ Δx\ \(b\_-\) + 60\ \(c\_-\) - 
            3\ Δx\ b\_0 + 15 == 
          0}, 
{\(a\_-\), \(b\_-\), \(c\_-\)}]\)

Put this result into the coefficient of the seventh derivative of u.

Out[65] /. Out[59] // 
      Simplify

Out[66] /. Out[59] // Simplify

Why does it not depend on a ratio?

Page 199.  Problem 4.7

\(Solve[{
2 \( b\_+\) - 2 \( b\_-\) == 
          2, 
\(b\_+\) + \(b\_-\) == 
          0}, 
{\(b\_+\), \(b\_-\)}]\)

Chapter 5

Page 219.

Inverse[{{1, 1, 1}, 
{(Subscript[x, j - 1]), (Subscript[x, j]), 
            (Subscript[x, j + 1])}, 
{((Subscript[x, j - 1])^2), ((Subscript[x, j])^2), 
            ((Subscript[x, j + 1])^2)}}] // Simplify

Out[
          305] . {(Subscript[x, j + 1]) - 
            (Subscript[x, j - 1]), (((((Subscript[x, j + 1])^2) - ((Subscript[x, j - 1])^2)))/2), (((((Subscript[x, j + 1])^3) - ((Subscript[x, j - 1])^3)))/3)} // Simplify

Page 232.

Problem 5.4

Integrate[(2/Δ)(x^2)(2 χ - ((Subscript[x, i - 1]) - Δx + (Subscript[x, i - 1]) - 
                  2  Δx))(2/Δ)(x^2)(2  χ - ((Subscript[x, i - 1]) + (Subscript[x, i - 1]) - Δx)),{χ, (Subscript[x, i - 1]) - 2  Δx, 
          (Subscript[x, i - 1])}]

Integrate[(2/Δ)(x^2)(2 χ - ((Subscript[x, i - 1]) - Δx + (Subscript[x, i - 1]) - 
                  2  Δx))((-4)/Δ)(x^2)(2  χ - ((Subscript[x, i - 1]) + (Subscript[x, i - 1]) - 2  Δx)),{χ, (Subscript[x, i - 1]) - 2  Δx, (Subscript[x, i - 1])}]

Integrate[(2/Δ)(x^2)(2 χ - ((Subscript[x, i - 1]) - Δx + (Subscript[x, i - 1]) - 
                  2  Δx))(2/Δ)(x^2)(2  χ - ((Subscript[x, i - 1]) - Δx + (Subscript[x, i - 1]) - 
                  2  Δx)),{χ, (Subscript[x, i - 1]) - 2  Δx, (Subscript[x, i - 1])}]

Page 234.

Problem 5.4

tentL[j_] := 
        Which[
j == -2, (2/Δ)(x^2) (χ - ((Subscript[x, i]) - Δ(x/2))) (χ - 
                (Subscript[x, i])), 
j == -1, ((-4)/Δ)(x^2) (χ - ((Subscript[x, i]) - Δx )) (χ - (Subscript[x, i])), 
j == 
            0, (2/Δ)(x^2) (χ - ((Subscript[x, i]) - Δx )) (χ - ((Subscript[x, i]) - Δ(x/2)))
];

tentR[j_] := 
        Which[
j == 
            2, (2/Δ)(x^2) (χ - 
                (Subscript[x, i])) (χ - ((Subscript[x, i]) + Δ(x/2))), 
j == 
            1, ((-4)/Δ)(x^2) (χ - 
                (Subscript[x, i])) (χ - ((Subscript[x, i]) + Δx )), 
j == 
            0, (2/Δ)(x^2) (χ - ((Subscript[x, i]) + Δ(x/2))) (χ - ((Subscript[x, i]) + Δx ))
];

With tent=tentL

Integrate[
      D[tent[-1], χ] D[ tent[-2], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[
      D[tent[-1], χ] D[ tent[-1], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[
      D[tent[-1], χ] D[ tent[0], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[
      tent[-1] tent[-2], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

Integrate[
      tent[-1] tent[-1], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

Integrate[
      tent[-1] tent[0], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

Integrate[
      D[tentL[0], χ] D[tentL[-2], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]
Integrate[
      D[tentL[0], χ] D[tentL[-1], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]
Integrate[
      D[tentL[0], χ] D[tentL[0], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[
      D[tentR[0], χ] D[tentR[0], χ], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      D[tentR[0], χ] D[tentR[1], χ], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      D[tentR[0], χ] D[tentR[2], χ], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]

Integrate[
      tentL[0] tentL[-2], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[0] tentL[-1], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[0] tentL[0], {χ, (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[
      tentR[0] tentR[0], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      tentR[0] tentR[1], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      tentR[0] tentR[2], {χ, (Subscript[x, i]), (Subscript[x, i]) + Δx}]

Series[(1/(4 (h^2))) (-u[x - 2  h] + 8  u[x - h] - 
            14  u[x] + 8  u[x + h] - u[x + 2  h]), {h, 0, 3}]

Problem 5.5

tentL[j_] := 
        Which[
j == -2, (2/Δ)(x^2) (χ - ((Subscript[x, i]) - Δ(x/2))) (χ - 
                (Subscript[x, i])), 
j == -1, ((-4)/Δ)(x^2) (χ - ((Subscript[x, i]) - Δx )) (χ - (Subscript[x, i])), 
j == 
            0, (2/Δ)(x^2) (χ - ((Subscript[x, i]) - Δx )) (χ - ((Subscript[x, i]) - Δ(x/2)))
];

tentR[j_] := 
        Which[
j == 
            0, (2/Δ)(x^2) (χ - ((Subscript[x, i]) + Δ(x/2))) (χ - ((Subscript[x, i]) + Δx )), 
j == 
            1, ((-4)/Δ)(x^2) (χ - 
                (Subscript[x, i])) (χ - ((Subscript[x, i]) + Δx )), 
j == 
            2, (2/Δ)(x^2) (χ - 
                (Subscript[x, i])) (χ - ((Subscript[x, i]) + Δ(x/2)))
];

Integrate[
      tentL[-1] D[tentL[-2], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]
Integrate[
      tentL[-1] D[tentL[-1], χ], {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]
Integrate[
      tentL[-1] D[tentL[0], χ], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

Integrate[
      tentL[-1] tentL[-2], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[-1] tentL[-1], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[-1] tentL[0], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

\(Series[\(1\/\(2\ h\)\) \((u[x - 2  h] - u[x])\) , {h, 0, 
        1}]\)

Integrate[
      tentL[0] D[tentL[-2], χ], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[0] D[tentL[-1], χ], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[0] D[tentL[0], χ], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

Integrate[
      tentR[0] D[tentR[0], χ], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      tentR[0] D[tentR[+1], χ], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      tentR[0] D[tentR[+2], χ], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]

Integrate[
      tentL[0] tentL[-2], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[0] tentL[-1], {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]
Integrate[
      tentL[0] tentL[0], {χ, (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[
      tentR[0] tentR[0], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      tentR[0] tentR[1], {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]
Integrate[
      tentR[0] tentR[2], {χ, (Subscript[x, i]), (Subscript[x, i]) + Δx}]

The expression given cannot be correct as a simple taylor expansion will show.

Series[(1/(6  h)) (u[x - 2  h] - 4  u[x - h] + 4  u[x + h] - 
            u[x + 2  h]), {h, 0, 2}]

Series[(1/(4  h)) (u[x - 2  h] - 4  u[x - h] + 4  u[x + h] - 
            u[x + 2  h]), {h, 0, 2}]

Problem 5.6

\(\(\( (*\ 
      i = element\ \(\(number\)\(.\)\)\ *) \)\(
\)\(tent[
          i_, j_, k_] := 
        Which[i == "L", 
Which[
j == 1 && k == 0, 
            1 - \(3\/Δx\^2\) \((x - xim1)\)\^2 + 
              2 \((x - xim1)\)\^3\/Δx\^3, 
j == 2 && k == 
                0, \(\((x - xim1)\)\^2\/Δx\^2\) \((3 - 
                  2 \((x - xim1)\)\/Δx)\), 
j == 1 && k == 
                1, \(\((1 - \((x - xim1)\)\/Δx)\)\^2\) \((x - xim1)\)\/Δx, 
j == 2 && 
              k == 
                1, \(\((x - xim1)\)\^2\/Δx\^2\) \((\((x - xim1)\)\/Δx - 
                  1)\)], 
i == "R", 
Which[
j == 1 && k == 0, 
            1 - \(3\/Δx\^2\) \((x - xi)\)\^2 + 
              2 \((x - xi)\)\^3\/Δx\^3, 
j == 2 && k == 
                0, \(\((x - xi)\)\^2\/Δx\^2\) \((3 - 
                  2 \((x - xi)\)\/Δx)\), 
j == 1 && k == 
                1, \(\((1 - \((x - xi)\)\/Δx)\)\^2\) \((x - xi)\)\/Δx, 
j == 2 && 
              k == 
                1, \(\((x - xi)\)\^2\/Δx\^2\) \((\((x - xi)\)\/Δx - 1)\)]
];\)\)\)

Integrations over left element.

Integrate[
      D[tent["L", 1, 0], x] tent["L", 2, 0], {x, xim1, 
        xim1 + Δx}]
Integrate[
      D[tent["L", 1, 1], x] tent["L", 2, 0], {x, xim1, 
        xim1 + Δx}]
Integrate[
      D[tent["L", 2, 0], x] tent["L", 2, 0], {x, xim1, 
        xim1 + Δx}]
Integrate[
      D[tent["L", 2, 1], x] tent["L", 2, 0], {x, xim1, 
        xim1 + Δx}]

Integrations over right element.

Integrate[
      D[tent["R", 1, 0], x] tent["R", 1, 0], {x, xi, 
        xi + Δx}]
Integrate[
      D[tent["R", 1, 1], x] tent["R", 1, 0], {x, xi, 
        xi + Δx}]
Integrate[
      D[tent["R", 2, 0], x] tent["R", 1, 0], {x, xi, 
        xi + Δx}]
Integrate[
      D[tent["R", 2, 1], x] tent["R", 1, 0], {x, xi, 
        xi + Δx}]

Series[((-(u'[x - h] - 2  u'[x] + u'[x + h]))/10) + (((u[x + h] - u[x - h]))/2), {h, 0, 3}]

I claim the value of 3 is incorrect for the following reason: Simply taylor expanding gives very incorrect results.

Series[((-(u'[x - h] - 3  u'[x] + u'[x + h]))/10) + (((u[x + h] - u[x - h]))/2), {h, 0, 3}]

Problem 5.8

Integrate[(y - 
              yjm1) (1 - (((y - yjm1))/Δ)x) + (x - 
              xim1) (1 - (((x - xim1))/Δ)x), {x, xim1, 
        xim1 + Δx}, {y, yjm1, 
        yjm1 + Δx}]

Integrate[(((y - 
                  yjm1))/Δ)x (1 - (((y - yjm1))/Δ)x) + ((((((x - xim1))/Δ)x))^2), {x, xim1, 
        xim1 + Δx}, {y, yjm1, 
        yjm1 + Δx}]

Integrate[((((((y - yjm1))/Δ)x))^2) + ((((((x - xim1))/Δ)x))^2), {x, xim1, xim1 + Δx}, {y, yjm1, 
        yjm1 + Δx}]

Problem 5.9

tent[i_] := 
        Which[
i == 
            1, (1/4) (1 - ξ) (1 - η), 
i == 
            2, (1/4) (1 + ξ) (1 - η), 
i == 
            3, (1/4) (1 + ξ) (1 + η), 
i == 
            4, (1/4) (1 - ξ) (1 + η)
];

Integrate[
      tent[3] tent[4], {ξ, -1, 1}, {η, -1, 1}]

c = {c1 -> (x2 - x1) (y4 - y1) - (x4 - x1) (y2 - 
                    y1), c2 -> (x2 - x1) (y3 - y2) - (x3 - 
                    x2) (y2 - 
                    y1), 
c3 -> (x3 - 
                    x4) (y3 - y2) - (x3 - x2) (y3 - y4), 
          c4 -> (x3 - x4) (y4 - y1) - (x4 - x1) (y3 - 
                    y4)};

4  c1 + 2  c2 + c3 + 2  c4 /. c // Simplify

2  c1 + 4  c2 + 2  c3 + c4 /. c // Simplify

Problem 5.14

Integrate[(((χ - 
                (Subscript[x, i]) + Δx))/Δ)x ((((Subscript[x, i]) - χ))/Δ)x, {χ, (Subscript[x, i]) - Δx, 
        (Subscript[x, i])}]

Integrate[((((χ - (Subscript[x, i]) + Δx))^2)/Δ)(x^2), {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]
Integrate[(((((Subscript[x, i]) + Δx - χ))^2)/Δ)(x^2), {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]

Integrate[((((Subscript[x, i]) + Δx - χ))/Δ)x (((χ - (Subscript[x, i])))/Δ)x, {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]

Integrate[(((χ - (Subscript[x, i]) + Δx))/Δ)x, {χ, (Subscript[x, i]) - Δx, (Subscript[x, i])}]

Integrate[χ - (Subscript[x, i]) + Δx, {χ, 
        (Subscript[x, i]) - Δx, (Subscript[x, i])}]
Integrate[
      (Subscript[x, i]) + Δx - χ, {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]

Integrate[
      (Subscript[x, i]) + Δx - χ, {χ, (Subscript[x, i]), 
        (Subscript[x, i]) + Δx}]

Problem 5.16

Integrate[(x - xi + Δx) (y - 
            yj + Δx) (1 - (((y - yj + Δx))/Δ)x), {x, xi - Δx, xi}, {y, 
        yj - Δx, yj}]

Problem 5.17

tent[i_] = 
        Which[i == 1, ((-1)/2) x (1 - x), i == 2, 
          1 - (x^2), i == 3, (1/2) x (1 + x)];

Integrate[tent[3] tent[1], {x, -1, 1}]

Integrate[tent[3] tent[2], {x, -1, 1}]

Integrate[tent[3]^2, {x, -1, 1}]

Integrate[tent[1]^2, {x, -1, 1}]

Integrate[tent[1] tent[2], {x, -1, 1}]

Integrate[tent[1] tent[3], {x, -1, 1}]

Integrate[
      tent[3] D[tent[1], x], {x, -1, 1}]
Integrate[
      tent[3] D[tent[2], x], {x, -1, 1}]
Integrate[tent[3] D[tent[3], x], {x, -1, 1}]

Integrate[
      tent[1] D[tent[1], x], {x, -1, 1}]
Integrate[
      tent[1] D[tent[2], x], {x, -1, 1}]
Integrate[tent[1] D[tent[3], x], {x, -1, 1}]

Integrate[tent[2] tent[1], {x, -1, 1}]
Integrate[tent[2] tent[2], {x, -1, 1}]
Integrate[tent[2] tent[3], {x, -1, 1}]

Integrate[
      tent[2] D[tent[1], x], {x, -1, 1}]
Integrate[
      tent[2] D[tent[2], x], {x, -1, 1}]
Integrate[tent[2] D[tent[3], x], {x, -1, 1}]

Page 255

fAB = (((2  u[x + h, y] + u[x + h, y - h] - 2  u[x, y] - u[x, y - h] + u[x + h, y + h] - u[x, y + h]))/(4 h));
fCD = (((2 (u[x, y] - u[x - h, y]) + u[x, y + h] - u[x - h, y + h] + u[x, y - h] - u[x - h, y - h]))/(4 h));

Series[fAB, {h, 0, 2}]
Series[fCD, {h, 0, 2}]

Series[fAB - fCD, {h, 0, 2}]

gBC = 
        (1/(4 h))(2 (u[x, y + h] - u[x, y]) + u[x + h, y + h] - 
            u[x + h, y] - u[x - h, y] + 
            u[x - h, y + h]);
gDA = (((2 (u[x, y] - u[x, y - h]) + u[x + h, y] - u[x + h, y - h] + u[x - h, y] - u[x - h, y - h]))/(4 h));

Series[gBC, {h, 0, 2}]
Series[gDA, {h, 0, 2}]

Series[gBC - gDA, {h, 0, 2}]

(fAB - fCD) h + (gBC - gDA) h // Simplify

Series[Out[14], {h, 0, 2}] // Simplify

Chapter 7

\(\(f[x_] := 
        Which[x ≤  .8, 
          0,  .8 ≤ x < 1, \((x -  .8)\)\/ .2, 
          1 ≤ x ≤ 1.2, \(-\((x - 1.2)\)\)\/ .2, 
          x > 1.2, 0];\)\)

initialCondition = 
        Plot[f[x], {x,  .5, 
            1.8}, 
PlotStyle -> {RGBColor[1, 0, 
                0]}];

n = 12;
Δx = (((1.7 - 0.5))/n);
x = Table[0.5 + Δx i, {i, 0, n}];

ic = Map[f, x];
ListPlot[ic];

Page 282.  Problem 7.1

Numerical scheme to use.

os[x_, σ_] := 
        Identity[
            x] - (σ/2) (RotateLeft[x] - 
                RotateRight[x]);

exercise[nts_, σ_] := 
        Module[{exact, ts, approx}, 
exact = 
            Plot[f[x - nts Δx σ], {x, 0.5, 
                1.8}, 
PlotStyle -> 
                RGBColor[0, 1, 0], 
DisplayFunction -> 
                Identity
]; 
ts = 
            Nest[os[#, σ] &, ic, nts]; 
approx = 
            ListPlot[
Thread[{x, 
                  ts}], 
PlotJoined -> True, 
              PlotRange -> 
                All, 
DisplayFunction -> 
                Identity]; 
Show[{initialCondition, exact, 
              approx}];
];

exercise[3, 
      1.0]

Table[exercise[i, 1.0], {i, 2, 5}];

Table[exercise[i, 0.5], {i, 2, 5}];

Table[exercise[i, 1.5], {i, 2, 5}];

If I make σ very small will this delay the istability?  No still get divergence from the expected wave form.

Table[exercise[i, 0.05], {i, 2, 15}];

Page 282.  Problem 7.2

Numerical scheme to use.

os[x_, σ_] := 
        Identity[
            x] - σ (Identity[x] - RotateRight[x]);

exercise[nts_, σ_] := 
        Module[{exact, ts, approx}, 
exact = 
            Plot[f[x - nts Δx σ], {x, 0.5, 
                1.8}, 
PlotStyle -> 
                RGBColor[0, 1, 0], 
DisplayFunction -> 
                Identity
]; 
ts = 
            Nest[os[#, σ] &, ic, nts]; 
approx = 
            ListPlot[
Thread[{x, 
                  ts}], 
PlotJoined -> True, 
              PlotRange -> 
                All, 
DisplayFunction -> 
                Identity]; 
Show[{initialCondition, exact, 
              approx}];
];

exercise[3, 
      1.0]

exercise[3, 0.5]

exercise[3, 1.5]

Table[exercise[i, 1.0], {i, 2, 5}];

Page 282.  Problem 7.3

Numerical scheme to use.

os[x_, σ_] := 
        Identity[
            x] - σ (Identity[x] - RotateRight[x]);

Get initial conditions.

a = 0.5;
b = 1.7;
getIC[ΔxIn_] := 
        Module[{n}, 
Δx = ΔxIn; 
n = (((b - a))/Δ)x; 
x = 
            Table[a + Δx i, {i, 0, 
                n}]; 
ic = 
            Map[f, x]
];

getIC[0.05]

ListPlot[
        ic];

exercise[nts_, σ_] := 
        Module[{exact, ts, approx}, 
exact = 
            Plot[f[x - nts Δx σ], {x, a, 
                b}, 
PlotStyle -> 
                RGBColor[0, 1, 0], 
DisplayFunction -> 
                Identity
]; 
ts = 
            Nest[os[#, σ] &, ic, nts]; 
approx = 
            ListPlot[
Thread[{x, 
                  ts}], 
PlotJoined -> True, 
              PlotRange -> 
                All, 
DisplayFunction -> 
                Identity]; 
Show[{initialCondition, exact, 
              approx}];
];

exercise[3, 
      0.5]

σExercise = 0.5;
exercise[4, σExercise]
exercise[8, σExercise]
exercise[16, σExercise]

getIC[0.025]

σExercise = 0.5;
exercise[4, σExercise]
exercise[8, σExercise]
exercise[16, σExercise]

Chapter 8

General Routines

discretize[a_, b_, nCells_] := 
      Module[{Δx}, 
        
Δx = (((b - a))/nCells); 
        
Return[
          Table[N[a + Δx i], {i, 0, nCells}]]; 
        
];

initialCondition[a_, b_, nCells_, iCFunction_] := 
      Module[{xIC, yIC}, 
xIC = discretize[a, b, nCells]; 
        
yIC = Map[iCFunction, xIC]; 
        
Return[Thread[{xIC, yIC}]]; 
];

Page 304

Series[((1 - 4  β Sin[(ϕ/2)]^2)/(E^(-β (ϕ^2)))), {ϕ, 0, 7}]

Page 306

G[σ_] := 
        Plot[(Sqrt[Cos[ϕ]^2 + (σ^2) Sin[ϕ]^2]), {ϕ, 
            0, π}, 
DisplayFunction -> 
            Identity];

Show[{G[0.25], G[0.5], G[0.75]}, 
        DisplayFunction -> $DisplayFunction];

These next plots don't look correct???

p[σ_] := 
        Plot[ArcTan[σ Tan[ϕ]]/(σ ϕ), {ϕ, 
            0, π}, 
DisplayFunction -> 
            Identity];
Show[{p[0.25], p[0.5], 
          p[0.75]}, 
DisplayFunction -> $DisplayFunction];

Page 307

<< Graphics`

p[σ_] := 
        PolarPlot[(Sqrt[1 - 4 σ (1 - σ) Sin[(ϕ/2)]^2]), {ϕ, 0, 
            2 π}, 
DisplayFunction -> 
            Identity];

Show[{p[0.25], p[0.5], 
          p[0.75]}, 
DisplayFunction -> $DisplayFunction];

\(temp[x_] := 4  x \((1 - x)\); 
temp[0.25]\)temp[0.75]

Why is this not correct???

p[σ_] := 
      PolarPlot[
        ArcTan[(((σ Sin[ϕ]))/(1 - σ + 
                σ Cos[ϕ]))]/(σ ϕ), {ϕ, 0, 
          2 π}, 
DisplayFunction -> Identity]; 
    
Show[{p[0.25], p[0.5], p[0.75]}, 
      
DisplayFunction -> $DisplayFunction];

Page 310

p[σ_] := 
        Plot[(Sqrt[1 - 4 (σ^2) (1 - (σ^2)) Sin[(ϕ/2)]^4]), {ϕ, 0, π}, 
DisplayFunction -> 
            Identity];

Show[{p[0.25], p[0.5], 
          p[0.75]}, 
DisplayFunction -> $DisplayFunction];

Why is this not correct???

\(p[σ_] := 
      Plot[ArcTan[
            \((σ\ Sin[ϕ])\)\/\(1 - 
                2  σ\^2\ Sin[ϕ\/2]\^2\)]\/\(σ\ ϕ\), {
          ϕ, 0, π}, 
        
DisplayFunction -> Identity]; 
    


    Show[{p[0.25], p[0.5], p[0.75]}, 
      
DisplayFunction -> $DisplayFunction]; \)

Page 314

\(pp[β_] := 
      Plot[\((2\ β\ Cos[ϕ] + 
              \@\(1 - 4\ \(β\^2\) Sin[ϕ]\^2\))\)\/\(1 + 
            2  β\), {ϕ, 0, π}, 
        
DisplayFunction -> Identity]; 
    


    Show[{pp[0.25], pp[0.5]}, 
      
DisplayFunction -> $DisplayFunction]; \)

Page 339

Problem 8.5

uExact[α_, k_, 
          t_] := 
Plot[(E^(-α (k^2) (π^2) t)) Sin[k π x], {x, 0, 
            1}, 
PlotStyle -> 
            RGBColor[0, 1, 0], 
DisplayFunction -> 
            Identity];

\(temp = Table[uExact[1\/3, 1, \(1\/10\) i], {i, 0, 10}]; 


    Show[temp, DisplayFunction -> $DisplayFunction]; \)

explicitCentral[v_List, β_] := 
        Module[{temp}, 
temp = (β RotateRight[
                    v] + (1 - 2 β) v + β RotateLeft[
                    v]); 
temp[[1]] = 
            0; 
temp[[Length[temp]]] = 
            0; 
Return[
            temp];
];

diffusionOneDimOneStep[a_, b_, iCFunction_, nCells_, scheme_, β_, 
        nTimeSteps_] := 

Module[{ic, 
          icx, icy}, 
ic = 
          initialCondition[a, b, nCells, 
            iCFunction]; 
icx = 
          Map[First, ic]; 
icy = 
          Map[Last, ic]; 
yAtnTimeSteps = 
          FixedPoint[scheme[#, β] &, icy, 
            nTimeSteps]; 
Return[
          Thread[{icx, yAtnTimeSteps}]];
]

Exer85[k_, nCells_, β_, nTimeSteps_] := 
        Module[{α = 1, Δx = (((1 - 0))/nCells), ue, au, 
            lpa}, 
ue = 
            uExact[1, 
              k, β  Δ((x^2)/α) 
                nTimeSteps]; 
ua = 
diffusionOneDimOneStep[0, 1, Sin[k π #1] &, nCells, 
              explicitCentral, β, nTimeSteps]; 
lpa = 
            ListPlot[ua, 
              DisplayFunction -> Identity]; 
Show[{ue, 
              lpa}, DisplayFunction -> $DisplayFunction];
];

Exer85[1, 30, (1/3), 5]
Exer85[1, 30, (1/3), 10]

Exer85[1, 30, (1/6), 5]
Exer85[1, 30, (1/6), 10]

Exer85[5, 30, (1/3), 5]
Exer85[5, 30, (1/3), 10]

Exer85[5, 30, (1/6), 5]
Exer85[5, 30, (1/6), 10]

Exer85[10, 30, (1/3), 5]
Exer85[10, 30, (1/3), 10]

Exer85[10, 30, (1/6), 5]
Exer85[10, 30, (1/6), 10]

Series[(((1 - 4 β Sin[(ϕ/2)]^2))/(E^(-β (ϕ^2)))), {ϕ, 0, 6}] // Simplify

Page 341

Problem 8.18.

Series[
        ArcSin[σ Sin[ϕ]]/(σ ϕ), {ϕ, 0, 4}] // 
      Simplify