微分方程数值解作业 1
微分方程
Problem 1
Question (a)
The interpolating linear function:
Substitute into (1.56)
Question (b)
Solve the 2nd order FDE
Assume
There is always one solution , therefore the method is not A-stable.
Question (c)
The order of the error of the linear approximation for :
The order of the truncation error of FDE: Second ordered method.
Problem 2
- Euler method
- Trapezoidal method
- Backward Euler method
Euler method
A-Stability: , conditionally A-Stable.
L-Stabliity: , not L-Stable.
Trapezoidal method
A-Stability: is true for all on the left half of the complex plane. It is a-stable.
L-Stability: for all on the left half of the complex plane but not for the right. Not L-Stable.
Backward Euler method
A-Stability: is true for all on the left half of the complex plane. It is A-stable even on partial right complex plane.
L-Stability: . It is L-stable.
Problem 3
Question (a)
FDE for and
Question (b)
Hamiltonian. Denote ,
Therefore, when the energy is conserved, it is required that
which crossponds to the trapezoid method. In Euler method,
The energy decreases. In backward Euler method,
The energy increases.
Question (c)
Deduce FDE for RK2. Denote
Only when the matrix on the last line is orthogonal. Therefore RK2 does not conserve the energy.
Problem 4
Question (a)
Use the one-sided backward numerical differentiation formula
Question (b)
This method is implicit. Soluiton on the model problem
Assume
When , both roots are less than 1. BDF2 is zero-stable.
Question (c)
Denote ,
Therefore BDF2 is A-stable.
Problem 5
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
| forwardEulerIt[f_, k_, t_, y1_] := y1 + k f[t, y1];
backwardEulerIt[f_, k_, t_, y1_] :=
NSolveValues[y2 == y1 + k f[t + k, y2], y2, Reals][[2]];
leapfrogIt[f_, k_, t_, y1_, y2_] := y1 + 2 k f[t, y1];
order[leapfrogIt] = 2;
trapezoidalIt[f_, k_, t_, y1_] :=
NSolveValues[y2 == y1 + 1/2 (f[t, y1] + f[t + k, y2]), y2,
Reals][[2]];
rk4It[f_, k_, t_, y1_] := Module[{k1, k2, k3, k4},
k1 = k f[t, y1];
k2 = k f[t + k/2, y1 + k1/2];
k3 = k f[t + k/2, y1 + k2/2];
k4 = k f[t + k, y1 + k3];
Return[y1 + 1/6 (k1 + 2 k2 + 2 k3 + k4)];
];
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| ivpOrder1Solve[f_, y0_, ts_, tt_, m_, method_] := Module[{k, y, ti},
k = (tt - ts)/m;
ti = Range[ts, tt, k];
y = ConstantArray[0, Length[ti]];
y[[1]] = y0;
Do[y[[i]] = method[f, k, ti[[i - 1]], y[[i - 1]]], {i, 2,
Length[ti]}];
Return[
Interpolation[Transpose[{ti, y}], InterpolationOrder -> 1][t]];
];
ivpOrder2Solve[f_, y0_, ts_, tt_, m_, method_] := Module[{k, y, ti},
k = (tt - ts)/m;
ti = Range[ts, tt, k];
y = ConstantArray[0, Length[ti]];
y[[1]] = y0;
y[[2]] = rk4It[f, k, ti[[1]], y[[1]]];
Do[y[[i]] = method[f, k, ti[[i - 2]], y[[i - 2]], y[[i - 1]]], {i,
3, Length[ti]}];
Return[
Interpolation[Transpose[{ti, y}], InterpolationOrder -> 1][t]];
];
|
1
2
3
4
5
| f[t_, y_] := 10 y (1 - y);
y0 = 1/10;
ts = 0;
tt = 1;
exact = DSolveValue[{y'[t] == f[t, y[t]], y[ts] == y0}, y[t], t]
|
1
2
| ndsolve =
NDSolveValue[{y'[t] == f[t, y[t]], y[ts] == y0}, y[t], {t, ts, tt}]
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| getResults[m_] := {
exact,
ndsolve,
ivpOrder1Solve[f, N[y0], ts, tt, m, forwardEulerIt],
ivpOrder1Solve[f, N[y0], ts, tt, m, backwardEulerIt],
ivpOrder2Solve[f, N[y0], ts, tt, m, leapfrogIt],
ivpOrder1Solve[f, N[y0], ts, tt, m, trapezoidalIt],
ivpOrder1Solve[f, N[y0], ts, tt, m, rk4It]
};
methodTitles = {"Exact", "NDSolve", "Forward Euler", "Backward Euler",
"Leapfrog", "Trapezoidal", "RK4"};
plotResults[m_] := Module[{results},
results = getResults[m];
Plot[results, {t, 0, 1}, PlotLegends -> methodTitles,
PlotLabel -> StringForm["M = ``", m]]
]
GraphicsColumn[plotResults /@ {10, 20, 40}]
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
| e1 = exact /. t -> 1;
getError[method_, m_] :=
Abs[(ivpOrder1Solve[f, N[y0], ts, tt, Round[m], method] /. t -> 1) -
e1];
errorplot =
LogLogPlot[{getError[forwardEulerIt, m],
getError[backwardEulerIt, m], getError[trapezoidalIt, m],
getError[rk4It, m]}, {m, 10, 1000},
PlotLegends -> {"Forward Euler", "Backward Euler", "Trapezoidal",
"RK4"}]
ticklines =
LogLogPlot[{1/m, 1/m^2, 1/m^4}, {m, 10, 1000}, PlotStyle -> Dashed,
PlotLegends -> "Expressions"]
Show[ticklines, errorplot]
|
