微分方程数值解作业 4

Problem 1

Numerical domain of dependence for each method

The CFL condition require that the numerical domain of dependence should contain the physical domain of dependence, which is . Therefore, for these methods,

1. Can be used. Satifies the CFL condition when
2. Cannot be used. It never satisfies the CFL condition.
3. Can be used. Satifies the CFL condition when
4. Cannot be used. It never satisfies the CFL condition.

Problem 2

Question (a)

Question (b)

The numerical domain of dependence for the Fromm method is

The CFL condition requires that the numerical domain of dependence should contain the physical domain of dependence, which is a single point . To satisfy the CFL condition, we need

that is

To discuss the stability of the Fromm method, we use the assumption

Substitute this into the Fromm method, we have

Divide by , we have

The amplification factor is

To be stable, we need , that is . Therefore, we have

Therefore,

Finally we have that the Fromm method is stable when .

To discuss the monotonicity of the Fromm method, it can be noticed that the four coefficients of the Fromm method can never be all positive at the same . That is, the Fromm method is not monotonic.

Question (c)

The Fromm method

Determine the truncation error by using Taylor expansion ( subscripts are omitted for simplicity)

Notice that

Therefore, we have

Problem 3

Question (a)

Use the Taylor expansion

Question (b)

The differential equation to solve is

On grid points , where , we use forward difference for and central difference for ,

Extract from the equation, we have

Substitute with formula in Question (a), and use to denote , we have

Question (c)

To discuss the stability, use the assumption

Substitute this into FDE, we have

The amplification factor is

To be stable, we need , that is

Therefore,

That is, the stable condition is .

Problem 4

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

makeTridiagonal[n_, {a_, b_, c_}] := 
  DiagonalMatrix[ConstantArray[c, n - 1], -1] + 
   DiagonalMatrix[ConstantArray[b, n]] + 
   DiagonalMatrix[ConstantArray[a, n - 1], 1];
upwind[n_, \[Lambda]_] := 
  makeTridiagonal[n, {0, 1 - \[Lambda], \[Lambda]}];
laxFriedrichs[n_, \[Lambda]_] := Module[{mat},
   mat = 
    makeTridiagonal[
     n, {1/2 (1 - \[Lambda]), 0, 1/2 (1 + \[Lambda])}];
   mat[[n, n - 1]] = \[Lambda];
   mat[[n, n]] = 1 - \[Lambda];
   Return[mat];
   ];
laxWendroff[n_, \[Lambda]_] := Module[{mat},
   mat = 
    makeTridiagonal[
     n, {-\[Lambda]/2 (1 - \[Lambda]), 
      1 - \[Lambda]^2, \[Lambda]/2 (1 + \[Lambda])}];
   mat[[n, n - 1]] = \[Lambda];
   mat[[n, n]] = 1 - \[Lambda];
   Return[mat];
   ];

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

advectionSolve[a_, g_, n_, l_, r_, m_, td_, method_] := 
  Module[{x, t, k, h, \[Lambda], mat, u},
   h = (r - l)/(n + 1);
   x = Range[l + h, r - h, h];
   k = td/m;
   t = Range[0, td, k];
   \[Lambda] = (a k)/h;
   mat = method[n, \[Lambda]];
   u = ConstantArray[0, {n, m + 1}];
   u[[All, 1]] = g /@ x;
   Do[u[[All, j]] = mat . u[[All, j - 1]], {j, 2, m + 1}];
   Return[
    Interpolation[
     Join[Flatten[MapThread[List, {Outer[List, x, t], u}, 2], 1], 
      Map[{{l, #}, 0} &, t], Map[{{r, #}, 0} &, t]], 
     InterpolationOrder -> 1]];
   ];

Test the code with the following example

1
2
3
4
5
6
7
8

a = 1;
g = Function[x, If[0 <= x <= 1, 1., 0.]];
l = -10;
r = 10;
td = 7;
sol = advectionSolve[a, g, 201, l, r, 72, td, laxFriedrichs];
Manipulate[
 Plot[sol[x, t], {x, -10, 10}, PlotRange -> {-1, 2}], {t, 0, 7}]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

exact = DSolveValue[{D[u[x, t], x] + D[u[x, t], t] == 0, 
    u[x, 0] == If[0 <= x <= 1, 1, 0]}, u[x, t], {x, t}];
makeFigure411[a_, g_, n_, l_, r_, ms_, td_, method_, exact_] := 
  GraphicsColumn[Function[m, Module[{xs, sol},
      xs = Subdivide[l, r, n + 1];
      sol = advectionSolve[a, g, n, l, r, m, td, method];
      ListLinePlot[{sol[#, td] & /@ xs, 
        exact /. {x -> #, t -> td} & /@ xs},
       PlotRange -> {-1, 1.2},
       DataRange -> {l, r},
       PlotLegends -> 
        Placed[{StringForm["M = ``", m], "Exact"}, {0.2, 0.1}], 
       PlotStyle -> {Thick, Directive[Black, Dashed]},
       AxesLabel -> {"x", "solution"},
       AspectRatio -> 1/3,
       ImageSize -> Medium
       ]
      ]] /@ ms, ImageSize -> Large];

1
2

makeFigure411[a, g, 201, l, r, {70, 72}, td, #, exact] & /@ {upwind, 
  laxWendroff, laxFriedrichs}

Solution at t=7 for upwind method

Solution at t=7 for Lax-Wendroff method

Solution at t=7 for Lax-Friedrichs method

For all cases, we have

That indicates all these methods are unstable. As a result, it can be observed that all three solutions oscillate at the discontinuity. When , we have

That indicates all these methods are stable. It can be seen that the oscillation no longer exists in the solution. Among these three methods, Lax-Wendroff method is not monotonic, while the other two methods are monotonic. This is reflected in the solution, where the Lax-Wendroff method has overshoots at the discontinuity, while other two methods do not.

许可协议

本文采用 署名-相同方式共享 4.0 国际 许可协议，转载请注明出处。

分享文章