微分方程数值解作业 2

Problem 1

To solve (2.4)

We can always obtain the equation at non-boundary positions,

The above FDE approximation has an error of

Question (a)

The first condition

and in Eqn. becomes

For the second condition

substitute with an one-sided approximation

we have

And the final equation becomes

The whole FDE system in matrix form

Question (b)

Therefore

FDE system in matrix form

Problem 2

Question (a)

Use central difference to approximate with

Substitute back

Question (b)

Given the differerntial equation

Simply susbtitute with in the above equation, and discard the higher order terms, we have

Question (c)

When

we will be able to seperate from , and the difference equation becomes

Collect terms on the left hand side

The above equation is linear, and could be written in matrix form. Denote

Assume the boundary conditions are and , we can obtain Equation at and

Finally, the FDE system is expressed in explicit matrix form

Question (d)

Therefore, when satisfies , Equation (2.4) will not contain term, and can be solve by the method in Question (b). That requires to be a constant.

Problem 3

At , use a third-ordered forward difference to approximiate and

And use a forth-ordered forward difference to approximiate

Substitute back to the equation we have one equation

with 5 variables to .

At to , use a forth-ordered central difference to approximiate and

and a forth-ordered central difference to approximiate

substitute them back we have equations

with variables from to .

Also, the boundry conditions provide 3 equations

Put Equations together, we could obtain a second-order non-linear FDE system with equations and variables to .

Problem 4

Question (a)

1
2
3
4

eqn = \[Epsilon] D[y[x], {x, 2}] - D[y[x], x] == -1 && y[0] == 1 && 
   y[1] == 3;
sol = (1 + # + (E^(#/\[Epsilon]) - 1)/(E^(1/\[Epsilon]) - 1)) &;
eqn /. {y -> sol}

1

True

Question (b)

Collect terms

At boundries

FDE system in matrix form

Since we have used first-ordered backward difference to approximate , the truncation error is .

Question (c)

It is obvious that the matrix in Equation meets the second condition of Theorem 2.1, therefore the matrix is invertible.

If we use Equation (2.16) to build the FDE system, Theorem 2.2 guarantees the existence of a unique solution if . That is, .

Question (d)

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

tridiagonalSolve[b_, a_, c_, z_] := Module[{y, v, w, n},
   n = Length[a];
   w = a[[1]];
   y = ConstantArray[0, n];
   y[[1]] = z[[1]]/w;
   v = ConstantArray[0, n];
   Do[
    v[[i]] = c[[i - 1]]/w;
    w = a[[i]] - b[[i]] v[[i]];
    y[[i]] = (z[[i]] - b[[i]] y[[i - 1]])/w;
    , {i, 2, n}];
   Do[
    y[[i]] = y[[i]] - v[[i + 1]] y[[i + 1]];
    , {i, n - 1, 1, -1}];
   Return[y];
   ];

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

solve1[\[Epsilon]_, n_] := Module[{c, a, b, h, x, y, z},
   h = 1/(n + 1);
   x = Range[h, 1 - h, h];
   b = ConstantArray[h + \[Epsilon], n];
   a = ConstantArray[-(h + 2 \[Epsilon]), n];
   c = ConstantArray[\[Epsilon], n];
   z = ConstantArray[-h^2, n];
   z[[1]] = -h^2 - h - \[Epsilon];
   z[[n]] = -h^2 - 3 \[Epsilon];
   y = tridiagonalSolve[b, a, c, z];
   Return[
    Interpolation[Join[{{0, 1}, {1, 3}}, Transpose[{x, y}]], 
     InterpolationOrder -> 1]];
   ];
solve2[\[Epsilon]_, n_] := Module[{c, a, b, h, x, y, z},
   h = 1/(n + 1);
   x = Range[h, 1 - h, h];
   b = ConstantArray[\[Epsilon] + h/2, n];
   a = ConstantArray[-2 \[Epsilon], n];
   c = ConstantArray[\[Epsilon] - h/2, n];
   z = ConstantArray[-h^2, n];
   z[[1]] = -h^2 - b[[1]];
   z[[n]] = -h^2 - 3 c[[n]];
   y = tridiagonalSolve[b, a, c, z];
   Return[
    Interpolation[Join[{{0, 1}, {1, 3}}, Transpose[{x, y}]], 
     InterpolationOrder -> 1]];
   ];

1
2
3
4
5
6
7
8
9
10
11

solveAndPlot[eps_, n_] := Module[{s1, s2},
  s1 = solve1[eps, n];
  s2 = solve2[eps, n];
  Plot[{symSol /. \[Epsilon] -> eps, s1[x], s2[x]}, {x, 0, 1}, 
   PlotLegends -> 
    Placed[{"Symbolic", "Method from (b)", "Equation (2.16)"}, Right],
    PlotLabel -> 
    StringForm["Solution for \[Epsilon]=``, N=``", eps, n]]
  ];
GraphicsColumn[{solveAndPlot[0.1, 10], solveAndPlot[0.1, 20], 
  solveAndPlot[0.1, 40]}]