03.3 – Randwertproblem: Diskretisierung und Systemmatrix

Es sei die Randwertaufgabe

$-{u^{\prime \prime }}+b{u^\prime }+cu = f\quad in\quad \Omega = \left( {0,1} \right)$

$u\left( 0 \right) = {g_0}$

$u\left( 1 \right) = {g_1}$

gegeben. $f$ , ${g_0}$ und ${g_1}$ seien so gewählt, das $u\left( x \right) = x\left( {1-x} \right){e^x}$ die exakte Lösung ist.

a) Sei zunächst $b = 0$ . Geben Sie das Gleichungssystem an, das bei der Diskretisierung der Randwertaufgabe mittels finiter Differenzen entsteht. Geben Sie die Konditionszahl der Systemmatrix für unterschiedliche Werte von $c$ und unterschiedliche Schrittweiten an (Matlab Befehl cond). Welchen Zusammenhang stellen Sie fest?

b) Sei nun $b > 0$ . Geben Sie jeweils die Systemmatrix für eine Diskretisierung mittels finiter Differenzen an, wobei der Term erster Ordnung einmal durch die Rückwärtsdifferenz und einmal durch die zentrale Differenz approximiert wird. Bestimmen Sie für beide Varianten und verschiedene Schrittweiten eine approximative Lösung und bestimmen Sie jeweils den maximalen Fehler. Wie ändert sich der Fehler bei Halbierung der Schrittweite?

Lösung

a )

$-{u^{\prime \prime }}+b{u^\prime }+cu = f\quad in\quad \Omega = \left( {0,1} \right),\quad u\left( 0 \right) = {g_0},\quad u\left( 1 \right) = {g_1}$

Wir wollen die Gleichung nun diskretisieren im Punkt ${x_k}$ .

Approximation der Ableitung mit einem Differenzenquotienten:

$-\frac{{{u_{k+1}}-2{u_k}+{u_{k-1}}}}{{{h^2}}}+c{u_k} = {f_k}$

$u\left( {{x_k}} \right) \approx {u_k},\quad {f_k} = f\left( {{x_k}} \right)$

Nun müssen wir die Werte ${u_k}$ bestimmen. Dabei müssen wir alle Werte gleichzeitig bestimmen und können nicht wie beim Euler-Verfahren rekursiv vorgehen. Es ist also ein Gleichungssystem aufzustellen und zu lösen. Wir schreiben als Matrix:

$\left( {\begin{array}{*{20}{c}} 1 & 0 & 0 & \cdots & 0 \\ {-\frac{1} {{{h^2}}}} & {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}} & 0 & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & 0 & {-\frac{1} {{{h^2}}}} & {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}} \\ 0 & \cdots & 0 & 0 & 1 \\ \end{array} } \right)\left( {\begin{array}{*{20}{c}} {{u_0}} \\ \vdots \\ \vdots \\ \vdots \\ {{u_N}} \\ \end{array} } \right) = \left( {\begin{array}{*{20}{c}} {{g_0}} \\ {{f_1}} \\ \vdots \\ {{f_{N-1}}} \\ {{g_1}} \\ \end{array} } \right)$

Um das Gleichungssystem zu vereinfachen, können wir die erste und die letzte Zeile herauslösen. Das resultierende Gleichungssystem hat dann direkt die praktische Eigenschaft, dass es tridiagonal und symmetrisch ist:

$\left( {\begin{array}{*{20}{c}} {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}} & 0 & \cdots & 0 \\ {-\frac{1} {{{h^2}}}} & {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}} & 0 & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & 0 & {-\frac{1} {{{h^2}}}} & {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}} \\ 0 & \cdots & 0 & {-\frac{1} {{{h^2}}}} & {c+\frac{2} {{{h^2}}}} \\ \end{array} } \right)\left( {\begin{array}{*{20}{c}} {{u_1}} \\ \vdots \\ \vdots \\ \vdots \\ {{u_{N-1}}} \\ \end{array} } \right)$

$= \left( {\begin{array}{*{20}{c}} {{f_1}+\frac{1} {{{h^2}}}{g_0}} \\ {{f_2}} \\ \vdots \\ {{f_{N-2}}} \\ {{f_{N-1}}+\frac{1} {{{h^2}}}{g_1}} \\ \end{array} } \right)$

MatLab code zur Berechnung der Konditionszahl:

function co = finDiff1(c)

sw=[0.2, 0.1, 0.05, 0.025, 0.0125, 0.00625];

for i = 1 : 6
	h = sw(i);
	N = 1 / h;
	v = [(c+2 / (h * h)) * ones(1, N-1)]';
	v1 = [-1 / (h * h) * ones(1, N-1)]';
	A = spdiags([v1, v, v1], -1 : 1, N-1, N-1); %weil man bei spa die Diagonalen
                                                         %eintragen kann
	co(i) = cond(full(A)); 			    %Umwandlung in normale Matrix,
                                                         %weil cond nur da geht
end

Zusammenhang:
Die Konditionszahl steigt quadratisch mit sinkender Schrittweite: $cond\left( A \right) = O\left( {{h^{-2}}} \right)$ .
Wenn das c größer wird, ist das Problem besser konditioniert.

b )

$-{u^{\prime \prime }}+b{u^\prime }+cu = f\quad in\quad \Omega = \left( {0,1} \right),\quad u\left( 0 \right) = {g_0},\quad u\left( 1 \right) = {g_1}$

Zentrale Differenz:

$-\frac{{{u_{k+1}}-2{u_k}+{u_{k-1}}}}{{{h^2}}}+b\frac{{{u_{k+1}}-{u_{k-1}}}}{{2h}}+c{u_k} = {f_k}$

Rückwärtsdifferenz:

$-\frac{{{u_{k+1}}-2{u_k}+{u_{k-1}}}}{{{h^2}}}+b\frac{{{u_k}-{u_{k-1}}}}{h}+c{u_k} = {f_k}$

Die Matrix sieht für die zentrale Differenz folgendermaßen aus:

$\left( {\begin{array}{*{20}{c}} {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}+\frac{b} {{2h}}} & 0 & \cdots & 0 \\ {-\frac{1} {{{h^2}}}-\frac{b} {{2h}}} & {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}+\frac{b} {{2h}}} & 0 & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & 0 & {-\frac{1} {{{h^2}}}-\frac{b} {{2h}}} & {c+\frac{2} {{{h^2}}}} & {-\frac{1} {{{h^2}}}+\frac{b} {{2h}}} \\ 0 & \cdots & 0 & {-\frac{1} {{{h^2}}}-\frac{b} {{2h}}} & {c+\frac{2} {{{h^2}}}} \\ \end{array} } \right)\left( {\begin{array}{*{20}{c}} {{u_1}} \\ \vdots \\ \vdots \\ \vdots \\ {{u_{N-1}}} \\ \end{array} } \right)$

$= \left( {\begin{array}{*{20}{c}} {{f_1}+\frac{1} {{{h^2}}}{g_0}} \\ {{f_2}} \\ \vdots \\ {{f_{N-2}}} \\ {{f_{N-1}}+\frac{1} {{{h^2}}}{g_1}} \\ \end{array} } \right)$

Die Matrix für die Rückwärtsdifferenz sieht folgendermaßen aus:

$\left( {\begin{array}{*{20}{c}} {c+\frac{2} {{{h^2}}}+\frac{b} {h}} & {-\frac{1} {{{h^2}}}} & 0 & \cdots & 0 \\ {-\frac{1} {{{h^2}}}-\frac{b} {{2h}}} & {c+\frac{2} {{{h^2}}}+\frac{b} {h}} & {-\frac{1} {{{h^2}}}} & 0 & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & 0 & {-\frac{1} {{{h^2}}}-\frac{b} {{2h}}} & {c+\frac{2} {{{h^2}}}+\frac{b} {h}} & {-\frac{1} {{{h^2}}}} \\ 0 & \cdots & 0 & {-\frac{1} {{{h^2}}}-\frac{b} {{2h}}} & {c+\frac{2} {{{h^2}}}+\frac{b} {h}} \\ \end{array} } \right)\left( {\begin{array}{*{20}{c}} {{u_1}} \\ \vdots \\ \vdots \\ \vdots \\ {{u_{N-1}}} \\ \end{array} } \right)$

$= \left( {\begin{array}{*{20}{c}} {{f_1}+\frac{1} {{{h^2}}}{g_0}} \\ {{f_2}} \\ \vdots \\ {{f_{N-2}}} \\ {{f_{N-1}}+\frac{1} {{{h^2}}}{g_1}} \\ \end{array} } \right)$

Es ergibt sich quadratische Konvergenz für die zentrale Differenz, normale (lineare) Konvergenz für die Rückwärtsdifferenz. Das liegt daran, dass man die erste Ableitung approximieren muss. Dieses geschieht bei der zentralen Differenz mit der Konvergenzordnung 2, bei der Rückwärtsdifferenz mit der Konvergenzordnung 1. Dieser Unterschied zieht sich dann durch den Rest der Rechnung durch.

Interpretation: $\left\| {u-{u_h}} \right\| \leq c{h^\alpha }$ mit

$\alpha = 1$ für die Rückwärtsdifferenz
$\alpha = 2$ für die zentrale Differenz

$\left. {\begin{array}{*{20}{c}} {{e_{h1}} = \left\| {u-{u_{h1}}} \right\| \leq ch_1^\alpha } \\ {{e_{h2}} = \left\| {u-{u_{h2}}} \right\| \leq ch_2^\alpha } \\ \end{array} } \right\}\frac{{{e_{h1}}}} {{{e_{{h_2}}}}} = {\left( {\frac{{{h_1}}} {{{h_2}}}} \right)^\alpha }\quad \Rightarrow \quad \alpha = \frac{{\ln \left( {\frac{{{e_{h1}}}} {{{e_{h2}}}}} \right)}} {{\ln \left( {\frac{{h1}} {{h2}}} \right)}}$

Mathematical Engineering

03.3 – Randwertproblem: Diskretisierung und Systemmatrix

Lösung

a )

b )

Ähnliche Artikel

Kommentar verfassen