微分方程数值解作业 3

Problem 1

Question (a)

Assume (3.30)

and substitute it into (3.51)

With setting , we have

Collecting terms with and on each side, we have

Use the identity , we have

Then use , we can get

Finally, we have

which contains the desired .

Question (b)

The stability condition is

Breakdown the absolute value, we have

Since , the right inequality is always true. For the left inequality, we have

which holds irrespective of the value of under the condition

Problem 2

Question (a)

The stencil is shown below.

Limits on and are

Question (b)

FDE for the Dufort-Frankel scheme is

Extract and use Taylor expansion ( subscript is omitted), we have

Use and drop higher order terms, we have

This is of order . Due to the presence of term, the method is not consistent when , where when . Therefore the method is only conditionally consistent.

Question (c)

The method is explicit. can be directly calculated from , and without solving a linear system or inverse function.

Question (d)

Use notation, the FDE can be written as

Extract , we have

Substitute with the generic solution , we have

Assume , and divide both sides by , we have

which yields the following 2nd order equation

The roots of this equation are

For the method to be stable, we need . For real roots where , we have

For complex roots where , we have

and

Note that guarantees that , and the method is unconditionally stable.

Question (e)

From , we have

This indicates that though deacys exponentially, it does so slowly when is large. Therefore, the method is stable but not L-stable.

Problem 3

Question (a)

From equation (3.22)

we can replace terms with Taylor expansion at , and obtain

Cancel duplicated terms, we have

where has an order of . Note that the fifth spatial derivative terms in also cancel out, so the remain in is instead of .

Question (b) (c)

Taking , we have

When applied to the heat equation

we have

assuming is smooth enough, we have

Therefore, we can see that the terms in cancel out, and the error is reduced to .

Question (d)

Apply Taylor expansion to , we have

omitted. Cancel out terms to get

And extract

It is possible but tedious to show that in can be finally reduced to in . is at least of order. In heat equation, we have , which indicates that when , the two terms in cancel out, and the error is reduced to .

Question (e)

In the above sub-figure, is fixed for each curve and is varing along the horizontal axis. This does not match situation in the question.

This sub-figure has fixed for each curve, and is decreasing along the horizontal axis. For a given , curve has larger and smaller . It can be observed that curve has larger error than curve at the same , indicating that decreasing while keeping increases the error. From this perspective, there is no contradiction.

Problem 4

Question (a)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

makeAMatrix[n_, \[Lambda]_] := 
  DiagonalMatrix[ConstantArray[1 - 2 \[Lambda], n]] + 
   DiagonalMatrix[ConstantArray[\[Lambda], n - 1], 1] + 
   DiagonalMatrix[ConstantArray[\[Lambda], n - 1], -1];
makeBMatrix[n_, \[Lambda]_] := 
  DiagonalMatrix[ConstantArray[1 + 2 \[Lambda], n]] + 
   DiagonalMatrix[ConstantArray[-\[Lambda], n - 1], 1] + 
   DiagonalMatrix[ConstantArray[-\[Lambda], n - 1], -1];
explicitIt[u1_, f1_, f2_, k_, a_, b_] := a . u1 - k f1;
implicitIt[u1_, f1_, f2_, k_, a_, b_] := Inverse[b] . (u1 - k f2);
crankNicolsonIt[u1_, f1_, f2_, k_, a_, b_] := 
  Inverse[b + IdentityMatrix[Length[b]]] . (a . u1 + u1 - k (f1 + f2));
heatSolve[f_, g_, n_, m_, td_, it_] := 
  Module[{x, t, h, k, \[Lambda], a, b, fval, u},
   h = 1/(n + 1);
   x = Range[h, 1 - h, h];
   k = td/m;
   t = Range[0, td, k];
   \[Lambda] = k/h^2;
   a = makeAMatrix[n, \[Lambda]];
   b = makeBMatrix[n, \[Lambda]];
   fval = Outer[f, x, t];
   u = ConstantArray[0, {n, m + 1}];
   u[[All, 1]] = g /@ x;
   Do[u[[All, j]] = 
     it[u[[All, j - 1]], fval[[All, j - 1]], fval[[All, j]], k, a, 
      b], {j, 2, m + 1}];
   Return[
    Interpolation[
     Join[Flatten[MapThread[List, {Outer[List, x, t], u}, 2], 1], 
      Map[{{0, #}, 0} &, t], Map[{{1, #}, 0} &, t]]]
    ];
   ];

Check the code with the following test case

1
2
3
4
5
6
7

f = Function[{x, t}, 0];
g = Function[{x}, N[Sin[2 \[Pi] x]]];
n = 20;
m = 5;
td = 0.1;
sol = heatSolve[f, g, n, m, td, explicitIt];
Plot3D[sol[x, t], {x, 0.`, 1.`}, {t, 0.`, 0.1`}, PlotRange -> Full]

Find the analytical solution

1
2
3
4

exactsol = 
 DSolveValue[{D[u[x, t], t] == D[u[x, t], {x, 2}], 
   u[x, 0] == Sin[2 \[Pi] x]}, u[x, t], {x, t}]
Plot3D[exactsol, {x, 0.`, 1.`}, {t, 0.`, 0.1`}, PlotRange -> Full]